1
$\begingroup$

In here, I asked a related (similar) question. Using ?NumericQ seems have slowed down the entire calculation.

But I now have a more serious problem:

B1998 = {14, 9, 7, 4, 2, 4, 4, 0, 4, 2, 0, 1, 0, 3, 0, 3, 1, 0, 0, 1, 1, 0, 1, 1, 2, 0, 1, 0, 1, 1};
kB1998 = 50;

myBBBlik[data_, kdata_, f0_?NumericQ, p_?NumericQ, w_?NumericQ, 
    mu_?NumericQ, theta_?NumericQ] :=
    Block[
    {K, fj, pj, loglik},
    K = kdata;
    fj = Prepend[PadRight[data, K], f0];
    pj = Table[
    (w*p^j*(1 - p)^(K - j)
    + (1 - w)
    *NProduct[mu + r*theta, {r, 0, j - 1}]
    *NProduct[1 - mu + r*theta, {r, 0, K - j - 1}]
    /NProduct[1 + r*theta, {r, 0, K - 1}]
    )*Binomial[K, j], {j, 0, K}];

    loglik = 
    LogGamma[Tr[fj] + 1] - Tr[LogGamma[fj + 1]] + Tr[fj*Log[pj]];
    loglik
];

NMaximize[{myBBBlik[B1998, kB1998, f0, p, w, mu, theta], f0 > 0, 0 < p < 1, 0 < w < 1, 0 < mu < 1, theta > 0}, {f0, p, w, mu, theta}] // AbsoluteTiming


{170.991900, {-49.4846, {f0 -> 0.0602345, p -> 0.0478278, w -> 0.522311, mu -> 0.259876, theta -> 0.11778}}}

More than three minutes.

myBBBlik2[data_, kdata_] :=
    Block[
    {K, fj, pj, loglik},
    K = kdata;
    fj = Prepend[PadRight[data, K], f0];
    pj = Table[
    (w*p^j*(1 - p)^(K - j)
    + (1 - w)
    *Product[mu + r*theta, {r, 0, j - 1}]
    *Product[1 - mu + r*theta, {r, 0, K - j - 1}]
    /Product[1 + r*theta, {r, 0, K - 1}]
    )*Binomial[K, j], {j, 0, K}];

    loglik = 
    LogGamma[Tr[fj] + 1] - Tr[LogGamma[fj + 1]] + Tr[fj*Log[pj]];
    loglik
];

NMaximize[{myBBBlik2[B1998, kB1998], f0 > 0, 0 < p < 1, 0 < w < 1, 0 < mu < 1, theta > 0}, {f0, p, w, mu, theta}] // AbsoluteTiming

{12.480022, {-49.4846, {f0 -> 0.0602345, p -> 0.0478278, w -> 0.522311, mu -> 0.259876, theta -> 0.11778}}}

Note that in the second case, I have change from NProduct to Product.

So maybe I should keep leaving the parameters out of the function? But then in here (a different model), it seems that sometimes, we just have to put the parameters in the function, as their arguments AND use ?NumericQ.

Now the extreme case:

NMaximize[{myBBBlik[B1998, kB1998, f0, p, w, mu, theta], f0 > 0, 0 < p < 1, 0 < w < 1, 0 < mu < 1, theta > 0}, {f0, p, w, mu, theta}, Method -> "RandomSearch"] // AbsoluteTiming   

I have added Method -> "RandomSearch", I waited and waited, and waited,

In this case, it seems that ?NumericQ has made it "impossible" to optimize????

enter image description here

Note the time stamp in the screenshot. I aborted it after nearly half an hour.

Why does ?NumericQ cause so much trouble??

What's the best way to construct function for NMaximize?? Do we leave variables out of the function, or in the function as arguments?

$\endgroup$
11
  • $\begingroup$ Could you perhaps try to come up with a MWE to demonstrate the problem? At the moment the question seems quite localized. $\endgroup$
    – Yves Klett
    Dec 24 '14 at 8:45
  • $\begingroup$ @YvesKlett I am not sure how "minimal" you would like it to be? I have given a function that I want to optimize, and compared with a variation of the function. And showed the timing. $\endgroup$ Dec 24 '14 at 14:03
  • 1
    $\begingroup$ You might get better answers if your question and code are concise and address a general problem. At the moment, parsing your code is not entirely straightforward. $\endgroup$
    – Yves Klett
    Dec 24 '14 at 14:19
  • 2
    $\begingroup$ It often happens that by searching for a simple example, you can discover the essence of what is causing the problem. Is it because of the form of the function you are maximizing? Is it the number of variables being maximized over? Is it the ranges of the parameters? If you can narrow this down, you may find that you can state the problem in a compelling manner. $\endgroup$
    – bill s
    Dec 24 '14 at 16:26
  • 1
    $\begingroup$ You are welcome - switching over to holiday mode now though over here ;-) $\endgroup$
    – Yves Klett
    Dec 24 '14 at 16:32
1
$\begingroup$

The reason is that you have an inefficient inner loop when you use myBBBlik. Every time you call myBBBlik, it takes a certain amount of time to calculate the result (on my computer, about 1/2 second of mucking around because your likelhood function is a bit complex).

But myBBBlik2 creates an algebraic expression once (?or twice?) and can substitute the numerical values in relatively quickly. So the "Prepend, PadRight, Table" stuff is skipped.

To demonstrate this, add the following to your code in both myBBBlik and myBBBlik2:

Dynamic@callCounter
callCounter = 1;

Inside the Block, add

...
Block[
{K, fj, pj, loglik},
callCounter = callCounter + 1;
K = kdata;
...

You will see that in myBBBlik, callCounter is incremented roughly every 1/2 second, as myBBBlik is called about that often. But when you use myBBBlik2, it is called only twice, because myBBBlik2 has symbolic output.

$\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.