Матрица Гессе в максимальном правдоподобии - Гаусс против R

Question

Я борюсь со следующей проблемой. В двух словах: два разных пакета программ (Gauss от Aptech и R) дают совершенно разные гессианские матрицы в процедуре максимального правдоподобия. Я использую ту же процедуру (BFGS), точно такие же данные, ту же формулу максимального правдоподобия (это очень простая модель логита) с точно такими же начальными значениями, и, как ни странно, я получаю одинаковые результаты для параметров и журнала. вероятность. Только гессенские матрицы различны по обеим программам, и поэтому оценка стандартных ошибок и статистического вывода различна.

В этом конкретном примере не наблюдается значительных отклонений, но каждое усложнение модели увеличивает разницу, поэтому, если я попытаюсь оценить свою окончательную модель, обе программы будут давать совершенно разные результаты.

Кто-нибудь знает, как обе программы отличаются тем, как они вычисляют гессенский и, возможно, правильный способ получения одинаковых результатов?

РЕДАКТИРОВАТЬ: В коде R (Гаусса), вектор X (alt) является независимой переменной, состоящей из вектора с двумя столбцами, где первый столбец целиком, а второй столбец ответы субъектов. Вектор y (itn) является зависимой переменной, состоящей из одного столбца с ответами субъектов. Пример (код R и набор данных) взят из http://www.polsci.ucsb.edu/faculty/glasgow/ps206/ps206.html, просто как пример для воспроизведения и выделения проблемы.

Я приложил оба кода (синтаксис Гаусса и 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

Р Гессиан:

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

Answer 1

Они просто масштабируются по-разному. Числа ГАУССА примерно в 1924 раза меньше, чем числа R.

Я думаю, что ГАУСС сохраняет числа в меньшем диапазоне для численной стабильности.