How to obtain this optimization result?

Question

mylik[data_, kdata_, model_: "BBB"] := 
Block[{K, pj, fj, N0, loglik, above, below},
    K = kdata;
    If[model == "BBB", 
        pj = Table[
        w*PDF[BinomialDistribution[K, p], j] 
        +
        (1 - w)*PDF[BetaBinomialDistribution[alpha, beta, K], j], 
        {j, 0, K}];
    ];
    fj = Prepend[PadRight[data, K], f0];
    N0 = Total[fj];
    above = LogGamma[N0 + 1]; 
    below = Total[LogGamma[fj + 1]]; 
    loglik = fj.Log[pj]; 
    loglik = above - below + loglik; 
    loglik
];

link3 = {679, 531, 379, 272, 198, 143, 99, 67, 46, 32, 22, 14, 9, 5, 3, 1};
klink3 = 20;
liklink3 = mylik[link3, klink3] // N;

NMaximize[{liklink3, f0 > 0, alpha > 0, beta > 0, 0 < w < 1, 0 < p < 1}, {f0, alpha, beta, w, p}, MaxIterations -> 1000]

NMaximize[{liklink3, f0 > 0, alpha > 0, beta > 0, 0 < w < 1, 0 < p < 1}, {f0, alpha, beta, w, p}, Method -> {"SimulatedAnnealing", "SearchPoints" -> 200}, MaxIterations -> 1000]

NMaximize[{liklink3, f0 > 0, alpha > 0, beta > 0, 0 < w < 1, 0 < p < 1}, {f0, alpha, beta, w, p}, Method -> {"RandomSearch", "SearchPoints" -> 200}, MaxIterations -> 1000]

NSolve[{0.1206 == alpha/(alpha + beta), 
0.116 == 1/(alpha + beta)}, {alpha, beta}];
mletmp = Flatten[{%, p -> 0.5396, w -> (1 - 0.9892), f0 -> 871.3}]
liklink3 /. mletmp

A lot of the time, I am only getting -46.0252 for the target function. I want to get -45.7408. enter image description here

This difference of the target function value may not be very big, but the estimate of f0 is significantly different.

I have tried all possible methods, with many starting points with no success. Is there anything else could be done?

Try differential evolution, as per this. Your code produces a torrent of errors on my computer, so I can't tell you experimentally what parameter values will be most suitable for differential evolution, but those given in that answer should give acceptable results in many cases. — Oleksandr R., Dec 16 '14 at 00:10
@OleksandrR. I did try differential evolution with 500 SearchPoints. It did not give the desired estimates. I tried to change some other setting, such as working precision. I also tried all methods with 800 searchpoints. No luck. I just dont really see that's necessary to use more than 1000 searchpoints. — Chen Stats Yu, Dec 16 '14 at 00:19
Generally it is not necessary to use extremely large numbers of search points (at least not with differential evolution; with RandomSearch this may be beneficial). However, the values of the tuning parameters $F$ ("ScalingFactor") and $C$ ("CrossProbability") are absolutely critical to the success or otherwise of differential evolution. Did you adjust those? — Oleksandr R., Dec 16 '14 at 01:32
I believe I had played around with ScalingFactor. I will read a bit on CrossProbability and have a go. It's just a bit frustrating that I can't get it to the desired values (when i know exactly what they are). — Chen Stats Yu, Dec 16 '14 at 01:43
In the first instance you might like to try "ScalingFactor" between 0.5 and 1, and "CrossProbability" either 0 to 0.1 or 0.9 to 1. (Try both possibilities. Small values are better for separable functions; large ones, for unseparable functions. Intermediate values do not work well in general, but despite this, the default in Mathematica is 0.5.) — Oleksandr R., Dec 16 '14 at 01:51

How to obtain this optimization result?

0 Answers0