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

Posted 9 years ago
6741 Views
|
|
0 Total Likes
|
 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!
 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}}} *)