Message Boards Message Boards

Is there a way to speed things up for this optimization with NIntegrate

Posted 10 years ago
mylik[data_, kdata_, model_] := 
      Block[{K0, pj, fj, N0, loglik, above, below},
       K0 = kdata;
       If[model == "LNB",
        pj = NIntegrate[
             N[PDF[BinomialDistribution[kdata, 1/(1 + Exp[-x])], #]
               *
               PDF[NormalDistribution[mu, sigma], x]],
             {x, -Infinity, Infinity}
             ] & /@ Range[0, kdata];
        ];
       fj = Prepend[PadRight[data, kdata], 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;

NMaximize[
 {mylik[link3, klink3, "LNB"], f0 > 0, sigma > 0},
 {f0, mu, sigma},
 MaxIterations -> 1000
 ]

At the moment, the code works fine. But it's just not fast enough.

Is there a way to Improve this?

Thanks!

POSTED BY: Casper YC

Not sure how much this will help but here are a few suggestions.

(1) Put all parameters into the objective function, including the ones being adjusted by the optimization.

(2) Insist that they be explicitly NumberQ to avoid any useless attempt at symbolic preprocessing (by NMaximize or Nintegrate).

(3) Tighten the bounds of integration, say to mu -+5. This seems not to detract from quality of the result and speeds the objective function evaluation.

(4) Maybe roll it down to 400 or so interations in the optimizer. Some experiments indicate that gives a result that does not differ by much from 200 iterations.

I also made the constraints weak since that's going to be done by the code anyway. With these changes, here is what I have.

Clear[mylik]
mylik[data_, kdata_, model_, f0_?NumberQ, mu_?NumberQ, 
   sigma_?NumberQ] := Block[{K0, pj, fj, N0, loglik, above, below},
   K0 = kdata;
   If[model == "LNB", 
    pj = NIntegrate[
         N[PDF[BinomialDistribution[kdata, 1/(1 + Exp[-x])], #]*
           PDF[NormalDistribution[mu, sigma], x]], {x, -Infinity, 
          Infinity}] & /@ Range[0, kdata];];
   fj = Prepend[PadRight[data, kdata], 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;

It takes a little over two minutes. I don't know if that would be viewed as a huge improvement but it seems to be faster than what I started out with.

Timing[
 NMaximize[{mylik[link3, klink3, "LNB", f0, mu, sigma], f0 >= 0, 
   sigma >= 0}, {f0, mu, sigma}, MaxIterations -> 400]]

(* Out[137]= {137.848000, {-47.234921188, {f0 -> 610.928712217, 
   mu -> -2.17379818051, sigma -> 0.987297937367}}} *)

Here is a useful speedup proposed by my colleague Michael Trott. We go back to integrating over the full range for this version.

Clear[mylik2, int]

int[kdata_, Y_, mu_, sigma_, x_] = 
  N[PDF[BinomialDistribution[kdata, 1/(1 + Exp[-x])], Y]*
    PDF[NormalDistribution[mu, sigma], x]];

mylik2[data_, kdata_, model_, f0_?NumberQ, mu_?NumberQ, 
   sigma_?NumberQ] := Block[{K0, pj, fj, N0, loglik, above, below},
   K0 = kdata;
   If[model == "LNB", 
    pj = NIntegrate[
         Evaluate[int[kdata, #, mu, sigma, x]], {x, -Infinity, 
          Infinity}, Method -> "DoubleExponential"] & /@ 
       Range[0, kdata];];
   fj = Prepend[PadRight[data, kdata], 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;

Timing[
 NMaximize[{mylik2[link3, klink3, "LNB", f0, mu, sigma], f0 >= 0, 
   sigma >= 0}, {f0, mu, sigma}, MaxIterations -> 400]]

(* Out[160]= {32.536000, {-47.2342781829, {f0 -> 610.950495085, 
   mu -> -2.17381590688, sigma -> 0.98731594371}}} *)
POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract