1
$\begingroup$

In short, the question is how to construct the function you want to optimize?

link3 = {679, 531, 379, 272, 198, 143, 99, 67, 46, 32, 22, 14, 9, 5, 3, 1};
klink3 = 20;
B1997 = {14, 10, 6, 3, 3, 6, 2, 3, 2, 1, 1, 1, 0, 3, 2, 0, 3, 1, 1, 1, 0, 0, 1, 2, 0, 0, 1};
kB1997 = 50;

mylik[data_, kdata_] := Block[{K, mypdf, fj, pj, loglik},
    K = kdata;
    mypdf = ProbabilityDistribution[
        Binomial[K, j]*(w*p^j*(1 - p)^(K - j) + (1 - w)*
        Beta[alpha + j, beta + K - j]/Beta[alpha, beta]), {j, 0, K, 1}
    ];
    pj = Probability[pj == Range[0, K], pj \[Distributed] mypdf];
    fj = Prepend[PadRight[data, K], f0];
    loglik = LogGamma[Total[fj] + 1] - Total[LogGamma[fj + 1]] + fj.Log[pj]
];

mylik2[data_, kdata_, f0_?NumericQ, w_?NumericQ, p_?NumericQ, 
    alpha_?NumericQ, beta_?NumericQ] := 
    Block[{K, mypdf, fj, pj, loglik},
    K = kdata;
    mypdf = ProbabilityDistribution[
        Binomial[K, j]*(w*p^j*(1 - p)^(K - j) 
        + 
        (1 - w)* Beta[alpha + j, beta + K - j]/Beta[alpha, beta]), {j, 0, K, 1}];
    pj = Probability[pj == Range[0, K], pj \[Distributed] mypdf];
    fj = Prepend[PadRight[data, K], f0];
    loglik = LogGamma[Total[fj] + 1] - Total[LogGamma[fj + 1]] + fj.Log[pj]
];


cons = {f0 > 0, alpha > 0, beta > 0, 0 < p < 1, 0 < w < 1};
pars = {f0, p, w, alpha, beta};

Using different ways of constructing the likelihood function,

OptMethod = "RandomSearch"

NMaximize[{mylik[link3, klink3], cons}, pars, Method -> {OptMethod, "SearchPoints" -> Automatic}] // AbsoluteTiming
NMaximize[{mylik[link3, klink3], cons}, pars, Method -> {OptMethod, "SearchPoints" -> 100}] // AbsoluteTiming
NMaximize[{mylik[B1997, kB1997], cons}, pars, Method -> {OptMethod, "SearchPoints" -> Automatic}] // AbsoluteTiming
NMaximize[{mylik[B1997, kB1997], cons}, pars, Method -> {OptMethod, "SearchPoints" -> 100}] // AbsoluteTiming

Just under 10 seconds.

NMaximize[{mylik2[link3, klink3, f0, w, p, alpha, beta], cons}, pars, Method -> {OptMethod, "SearchPoints" -> Automatic}] // AbsoluteTiming
NMaximize[{mylik2[link3, klink3, f0, w, p, alpha, beta], cons}, pars, Method -> {OptMethod, "SearchPoints" -> 100}] // AbsoluteTiming
NMaximize[{mylik2[B1997, kB1997, f0, w, p, alpha, beta], cons}, pars, Method -> {OptMethod, "SearchPoints" -> Automatic}] // AbsoluteTiming
NMaximize[{mylik2[B1997, kB1997, f0, w, p, alpha, beta], cons}, pars, Method -> {OptMethod, "SearchPoints" -> 100}] // AbsoluteTiming

Almost a minute using the second method.

It seems there is clearly a difference. So what are the general rules to construct a function for optimization?

$\endgroup$
  • $\begingroup$ Some of the timing is overhead. For example when I removed the NumericQ tests in mylik2 it runs as fast as mylik $\endgroup$ – John Morganthau Dec 21 '14 at 2:17
  • $\begingroup$ But a lot of the time people are telling me I should leave ?NumericQ in the function as that is the "proper" way to write a function? It also ensures wired behavior when using optimization, such as "the function is not a real number at ........", sometimes, NMAximize does not replace the parameters using their numerical values. $\endgroup$ – Chen Stats Yu Dec 21 '14 at 14:58
  • $\begingroup$ I seriously doubt repeated calls to NumericQ is the problem. It is likely that without NumericQ NMaximize is able to come up with an optimized expression for your function symbolically which it caches. With NumericQ it can't optimize/evaluate the expression symbolically. $\endgroup$ – Andy Ross Dec 24 '14 at 18:37

Your Answer

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

Browse other questions tagged or ask your own question.