R复刻VBA牛顿-拉夫逊迭代模型：浮点精度差异致计算异常排查

我正在从零开发R代码，复刻VBA中基于牛顿-拉夫逊（Newton Raphson）迭代的模型，该模型需多次迭代直至输出差异稳定。目前发现：

• 从Excel读取输入值E时，R中显示14位小数（如365318.42516561900266），而Excel仅显示9位（如365318.425165619）
• 后续alpha计算结果与VBA存在细微差异：R输出为-3.98568702613257741518、-3.98577596443294002171；VBA输出为-3.98568702613256、-3.98577596443293
• 怀疑是R与Excel的运算机制差异导致模型运行失败

R代码 #

# 创建alpha偏差的dobj
APCI_summary$dobj <- 2 * APCI_summary$weight * (APCI_summary$final_exposure * APCI_summary$m_xy - APCI_summary$calc_deaths)
dalphaobj_values <- aggregate(APCI_summary$dobj, by = list(APCI_summary$age), FUN = sum)$x
dalphaobj_m <- matrix(dalphaobj_values, nrow = length(dalphaobj_values), ncol = 1)

# 创建alpha偏差的d2obj
APCI_summary$d2obj <- 2 * APCI_summary$weight * APCI_summary$final_exposure * APCI_summary$m_xy
d2alphaobj_values <- aggregate(APCI_summary$d2obj, by = list(APCI_summary$age), FUN = sum)$x
d2alphaobj_m <- diag(d2alphaobj_values)

# 创建alpha惩罚项的dobj和d2obj
alpha_table <- data.frame(age = seq(age_min, age_max), alpha_x = APCI_summary$alpha_x[APCI_summary$year == year_max])
alpha_x_values <- alpha_table$alpha_x

# 将提取的列转换为矩阵
Alpha <- matrix(alpha_x_values, nrow = length(alpha_x_values), ncol = 1)
tAlpha <- t(Alpha)

dP_temp_alpha <- 2 * lambda_alpha * tDiff_3 %*% Diff_3 %*% Alpha

d2P_temp_alpha <- 2 * lambda_alpha * tDiff_3 %*% Diff_3

# 分别合并dobj和d2obj
dalpha_total_m <- dalphaobj_m + dP_temp_alpha
d2alpha_total_m <- d2alphaobj_m + d2P_temp_alpha

# 最终更新alpha值的步骤
d2Obj_inv <- solve(d2alpha_total_m) %*% dalpha_total_m
Alpha <- Alpha - d2Obj_inv
Alpha_data <- data.frame(age = seq(age_min, age_max), alpha_x = Alpha)

# 使用新Alpha值更新计算
APCI_summary <- APCI_summary %>%
  left_join(Alpha_data, by = "age") %>%
  mutate(alpha_x = coalesce(alpha_x.y, alpha_x.x)) %>%
  select(-alpha_x.x, -alpha_x.y)

APCI_summary$log_m_xy <- log_m_xy(APCI_summary$alpha_x, APCI_summary$beta_x, APCI_summary$year, APCI_summary$kappa_x, APCI_summary$gamma_x)
APCI_summary$m_xy <- exp(APCI_summary$log_m_xy)

APCI_summary$deviance <- APCI_summary$weight * (APCI_summary$calc_deaths * log(APCI_summary$calc_deaths)
                                                - APCI_summary$calc_deaths
                                                - APCI_summary$calc_deaths * log(APCI_summary$final_exposure * APCI_summary$m_xy)
                                                + APCI_summary$final_exposure * APCI_summary$m_xy)

VBA代码 #

Private Sub ImproveAlpha(Lambda)

Dim X, Y, W, C, dlogm, X1, X2
Dim dObj(), d2Obj() As Double    ' 目标函数的一阶和二阶微分
ReDim dObj(RP.MinAge To RP.MaxAge, 1 To 1)
ReDim d2Obj(RP.MinAge To RP.MaxAge, RP.MinAge To RP.MaxAge) As Double

' 偏差项的贡献
For X = RP.MinAge To RP.MaxAge
    For Y = RP.MinYear To RP.MaxYear
        W = RP.Weight(Y)
        C = Y - X
        If RP.MinCohort <= C And C <= RP.MaxCohort Then
            dlogm = 1
            dObj(X, 1) = dObj(X, 1) + 2 * W * (E(X, Y) * m(X, Y) - D(X, Y)) * dlogm
            d2Obj(X, X) = d2Obj(X, X) + 2 * W * (E(X, Y) * m(X, Y)) * dlogm ^ 2
        End If
    Next Y
Next X


' 正则化惩罚项的贡献
Dim Diff, P, dP, d2P, dP_temp

Diff = DifferenceMatrix(RP.NumAges, ORDERALPHA)
P = MMult(Transpose(Diff), Diff)

ReDim dP(RP.MinAge To RP.MaxAge, 1 To 1)
ReDim d2P(RP.MinAge To RP.MaxAge, RP.MinAge To RP.MaxAge)
dP_temp = MMult(P, VectorToMatrix(Alpha))


For X1 = RP.MinAge To RP.MaxAge
    dObj(X1, 1) = dObj(X1, 1) + 2 * Lambda * dP_temp(X1 + 1 - RP.MinAge, 1)
    For X2 = RP.MinAge To RP.MaxAge
        d2Obj(X1, X2) = d2Obj(X1, X2) + 2 * Lambda * P(X1 + 1 - RP.MinAge, X2 + 1 - RP.MinAge)
    Next X2
Next X1


' 更新参数
Dim Delta
Delta = MMult(MInverse(d2Obj), dObj)
For X = RP.MinAge To RP.MaxAge
    Alpha(X) = Alpha(X) - Delta(1 + X - RP.MinAge, 1)
Next X

End Sub

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