0
$\begingroup$
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?

$\endgroup$
0
$\begingroup$
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!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.