Getting standard errors (Hessian matrix) from optimization with constraints

Question

ClearAll["Global`*"]
link1 = {84, 54, 36, 21};
klink1 = 4;
skinks = {56, 19, 28, 18, 24, 14, 9};
kskinks = 7;


mylikBetaBin[data_, kdata_, f0_, alpha_, beta_] := Module[
    {K, pj, fj, j, N, loglik, above, below},

    K = kdata;
    pj = Table[
        Binomial[K, j]*
        Beta[alpha + j, K - j + beta]/
        Beta[alpha, beta], {j, 0, K}
    ];

    fj = Prepend[data, f0];

    N = Sum[fj[[j]], {j, 1, Length[fj]}];
    above = LogGamma[N + 1];
    below = Sum[LogGamma[fj[[j]] + 1], {j, 1, Length[fj]}];
    loglik = fj.Log[pj];

    loglik = above - below + loglik
]

(*mylikBetaBin[link1, klink1, 37, 0.8, 0.3]*)

testfun[f0_, alpha_, beta_] := 
mylikBetaBin[skinks, kskinks, f0, alpha, beta];

(*testfun[165.178,0.573,1.576]*)

NMaximize[
    {testfun[f0, alpha, beta], 
    f0 >= 0 && alpha >= 0 && beta >= 0},
    {f0, alpha, beta},
    Method -> {"RandomSearch", "SearchPoints" -> 250}
]

I am new to Mathematica. I need to do some (global) optimization with constraints. This is a simple example I just managed to get it work! (finally)

I want to know if I am doing the right thing? Is NMaximize the function to use?

But the main question is, is there an optimizer will give me the standard errors (Hessian evaluated at these points)?

If not, how would I get them?

Answer 1

ans=NMaximize[
    {testfun[f0, alpha, beta], 
    f0 >= 0 && alpha >= 0 && beta >= 0},
    {f0, alpha, beta},
    Method -> {"RandomSearch", "SearchPoints" -> 200}
]


he = D[testfun[f0, alpha, beta], {{f0, alpha, beta}, 2}]

he /. ans[[2]] // MatrixForm

Inverse[%]
Diagonal[-%]
Sqrt[%]

After some trial and error, I believe I found 'a' way to do this.

If anyone has a better solution, please let me know. Thanks!