R语言中函数内调用lprim与脚本直接调用结果不一致的排障求助

我在求解参数theta时遇到一个问题：在Newton函数内部调用lprim时，theta的计算结果几乎为0，函数随即终止返回该值；但在函数外部直接调用lprim，传入相同参数得到的结果却是44.55，不知道该从何处排查Bug。

相关代码如下：

set.seed(0)

lprim <- function(X, theta, sigma) {
  n = length(X)
  return((1 / sigma) * (n - 2 * sum(exp(-(X - theta) / sigma) / (1 + exp(-(X - theta)/sigma)))))
}
lbis = function(X, theta, sigma){
  return(-2 / sigma ^ 2 * sum( exp(-(X - theta) / sigma) / (1 + exp(-(X - theta)/sigma))^2))
}

X = rlogis(50, 1, 1)
Newton = function(start, koniec, X, sigma) {
  i = 0
  theta = start
  while (abs(lprim(X, theta, sigma)) > 0.000001 & i < 5000){
    a = lprim(X, theta, sigma)
    b = lbis(X, theta, sigma)
    theta = theta - a / b
    i = i + 1
  }
  x = abs(lprim(X, theta, sigma))
  return(c(theta, sigma))
}
debug(lprim)
z = Newton(1, 5000, X, 1)
abs(lprim(X, z[0], z[1]))

排查要点 #

• 索引错误：R语言向量索引从1开始，你最后调用abs(lprim(X, z[0], z[1]))用了z[0]，这会返回NA，实际应该用z[1]（Newton返回的第一个元素是theta）。这是外部调用结果异常的直接原因。
• 迭代过程监控：在Newton函数的循环内添加打印语句，比如print(c(i, theta, a, b))，查看每一步的theta、lprim值(a)和二阶导数值(b)，确认theta是否在迭代中骤变到0。
• 二阶导数符号验证：Newton迭代公式为theta = theta - f(theta)/f'(theta)，你的lbis返回负数，逻辑分布对数似然的二阶导数确实为负，公式逻辑正确，但可手动计算几个点验证lbis的计算是否准确。
• 数值稳定性优化：当(X - theta)/sigma绝对值较大时，exp(-(X - theta)/sigma)可能出现数值溢出/下溢。可将表达式简化为1/(1 + exp((X - theta)/sigma))替换原有的exp(-x)/(1+exp(-x))，提升计算稳定性。

内容的提问来源于stack exchange，提问作者Kat