r - 最大似然的 Hessian 矩阵 - Gauss vs. R

Question

我正在努力解决以下问题。简而言之:两个不同的软件包(Aptech 和 R 的 Gauss)在最大似然过程中产生完全不同的 Hessian 矩阵。我使用相同的程序(BFGS)，完全相同的数据，相同的最大似然公式(它是一个非常简单的 logit 模型)，具有完全相同的起始值，令人困惑的是，我得到了相同的参数和对数结果 -可能性。只有 Hessian 矩阵在两个程序中不同，因此标准误差和统计推断的估计不同。

在这个特定的例子中，它没有出现太大的偏差，但是模型的每一个增加的复杂性都会增加差异，所以如果我尝试估计我的最终模型，两个程序都会产生完全不同的结果。

有谁知道，这两个程序在计算 Hessian 的方式以及获得相同结果的正确方法上有何不同？

编辑:在 R (Gauss) 代码中，向量 X (alt) 是自变量，由双列向量组成，第一列完全是一列，第二列是受试者的 react 。向量 y (itn) 是因变量，由一列包含受试者的 react 组成。示例(R 代码和数据集)取自 http://www.polsci.ucsb.edu/faculty/glasgow/ps206/ps206.html ，仅作为重现和隔离问题的示例。

我附上了代码(Gauss 和 R 语法)和输出。

任何帮助将不胜感激。谢谢 :)

高斯:

start={ 0.95568840 , -0.20459156 };

library maxlik,pgraph;
maxset;
_max_Algorithm = 2;
_max_Diagnostic = 1;
{betaa,f,g,cov,ret} = maxlik(XMAT,0,&ll,start);
call maxprt(betaa,f,g,cov,ret);
print _max_FinalHess;

proc ll(b,XMAT);
local exb, probo, logexb, yn, logexbn, yt, ynt, logl;

    exb = EXP(alt*b);
    //print exb;
    probo = exb./(1+exb);

    logexb = ln(probo);

    yn = 1 - itn;
    logexbn = ln(1 - probo);

    yt = itn';
    ynt = yn';

    logl = (yt*logexb + ynt*logexbn);

    print(logl);

retp(logl);
endp;

回复:

startv <- c(0.95568840,-0.20459156)

logit.lf <- function(beta) {

    exb <- exp(X%*%beta) 
    prob1 <- exb/(1+exb) 

    logexb <- log(prob1)

    y0 <- 1 - y
    logexb0 <- log(1 - prob1)

    yt <- t(y)
    y0t <- t(y0)

    logl <- -(yt%*%logexb + y0t%*%logexb0)

    return(logl)
}

logitmodel <- optim(startv, logit.lf, method="BFGS", control=list(trace=TRUE, REPORT=1), hessian=TRUE)
logitmodel$hessian

高斯输出:

return code =    0
normal convergence

Mean log-likelihood       -0.591820
Number of cases     1924

Covariance matrix of the parameters computed by the following method:
Inverse of computed Hessian

Parameters    Estimates     Std. err.  Est./s.e.  Prob.    Gradient
------------------------------------------------------------------
P01              2.1038        0.2857    7.363   0.0000      0.0000
P02             -0.9984        0.2365   -4.221   0.0000      0.0000

高斯黑森:

0.20133256       0.23932571 
0.23932571       0.29377761 

R 输出:

initial  value 1153.210839 
iter   2 value 1148.015749
iter   3 value 1141.420328
iter   4 value 1138.668174
iter   5 value 1138.662148
iter   5 value 1138.662137
iter   5 value 1138.662137
final  value 1138.662137 
converged

          Coeff.  Std. Err.          z       p value
[1,]  2.10379869 0.28570765  7.3634665 1.7919000e-13
[2,] -0.99837955 0.23651060 -4.2212889 2.4290942e-05

R黑森州:

          [,1]      [,2]
[1,] 387.34106 460.45379
[2,] 460.45379 565.24412

Answer 1

它们只是按比例不同。 GAUSS 数比 R 数小 1924 倍左右。
我认为 GAUSS 将数字保持在较小的范围内以保证数值稳定性。