# Performance in calculating maximum-likelihood- based estimates

My probability density function is a complicated one for which numerical estimation is necessary. Here my pdf:

pdf[s_?NumericQ] :=
combn  NIntegrate[
N[q^k (1 - q)^(n - k), 100] (
E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]))/
NIntegrate[
E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]), {q, 1/(
2 n + 1), 1 - 1/(2 n + 1)}, MaxRecursion -> 12] , {q, 1/(
2 n + 1), 1 - 1/(2 n + 1)}, MaxRecursion -> 15]


Given some (realistic) values for the other parameters:

n = 25000;
k = 24991;
\[Mu] = 10^-4;
\[Nu] = 10^-4;
q = k/n;
combn = Binomial[n, k];


I can find the maximum likelihood estimate (MLE) over a given range of values with

MLE = FindMaximum[pdf[s], {s, 0.1, 0.6}] # Be carefull those calculations are slow!


Then, I can find the 95% confidence interval by doing

t = Table[pdf[s], {s, 0.1, 0.6, 0.01}]
Table[s, {s, 0.1, 0.6, 0.01}][[Flatten[
Position[t, _?(# > MLE[] - 1.92 &)]]]]


But this is really not optimal because I have to recalculate pdf[s] for all s. If I were to calculate MLE from my t only, then my estimate wouldn't be as accurate as with FindMaximum because FindMaximum calculate pdf[s] over smaller increments as it get closer to MLE (at least that is what I assume, let me know if this is wrong).

What better (faster) solution can I use to get a good estimate of MLE and not-to-bad estimates of the boundaries of the confidence interval?

• have you tried calculating a few dozen points of pdf[s] and generating substituting a Interpolation function to replace use of pdf[s]? Feb 13, 2015 at 21:45
• Sounds like an interesting idea. I don't quite know what precision I would get with this method compare to using FindMaximum. I also don't know how to implement that in Mathematica Feb 13, 2015 at 22:02
• something is strange: sample points of your pdf[s]: {{0.005, -28.5331}, {0.0670833, -6.87258}, {0.129167, -3.72334}, \ {0.19125, -2.68104}, {0.253333, -2.37156}, {0.315417, -2.38988}, \ {0.3775, -2.57275}, {0.439583, -2.84284}, {0.501667, -3.15976}, \ {0.56375, -3.50092}, {0.625833, -3.85318}, {0.687917, -4.20862}, \ {0.75, -4.56239}} Feb 13, 2015 at 22:23
• a PDF should be non-negative everywhere Feb 13, 2015 at 22:23
• used this bit of code: pdf0Table = Block[ {sLow = 0.005, sHigh = 0.75, sStepNum = 12, sStep}, sStep = (sHigh - sLow)/sStepNum; Table[{s, pdf[s]}, {s, sLow, sHigh, sStep}] ] Feb 13, 2015 at 22:25

Change the method used by NIntegrate:

pdf[s_?NumericQ] :=
combn NIntegrate[
N[q^k (1 - q)^(n - k),
100] (E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]))/
NIntegrate[
E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]), {q,
1/(2 n + 1), 1 - 1/(2 n + 1)}, MaxRecursion -> 12,
Method -> "DoubleExponential"], {q, 1/(2 n + 1),
1 - 1/(2 n + 1)}, MaxRecursion -> 15,
Method -> "DoubleExponential"]


Note - I did nothing else to your PDF function. It is a seriously bad idea to have the internals of a function depend on external symbols, in general. Consider rewriting this...

That said, timings (on a netbook, using your example parameters):

pdf[#] & /@ {.05, .1, .3, .6, .9} // Timing

OLDpdf[#] & /@ {.05, .1, .3, .6, .9} // Timing

(*

{0.920406,{0.000157606,0.00856213,0.0938841,0.0245793,0.00454049}}

{26.644971,{0.000157606,0.00856213,0.0938841,0.0245793,0.00454049}}

*)

• Thank you! When you say It is a seriously bad idea to have the internals of a function depend on external symbols does it mean that I should just pass the internal symbols as arguments to the pdf, that is push pdf to play with local variable rather than global variable? To the exception of combn, I don't have any constant (independent of s) in my pdf that I could take out of the function. Feb 14, 2015 at 22:25
• @Remi.b: Yes, pass them all - makes things more readable (e.g., my first glance, seeing combin had me asking "what's this? Some other function?..."). Also makes things more maintainable. Fine to ignore for quick&dirty (I do), but if you're going to be using it for more than throw-away, better to do so (not to mention, makes reading the question easier ;-) ).
– ciao
Feb 14, 2015 at 22:46

Far from what the OP asked... Crude table based linear interpolation approach to learning more about the PDF:

n = 25000;
k = 24991;
\[Mu] = 10^-4;
\[Nu] = 10^-4;
q = k/n;
combn = Binomial[n, k];
pdf[s_?NumericQ] :=
combn
NIntegrate[N[q^k (1 - q)^(n - k), 100]
(E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]))/
NIntegrate[E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]),
{q, 1/(2 n + 1), 1 - 1/(2 n + 1)}, MaxRecursion -> 12],
{q, 1/(2 n + 1), 1 - 1/(2 n + 1)}, MaxRecursion -> 15]

pdf0Table = Block[{sLow = 0.001, sHigh = 3., sStepNum = 12, sStep},
sStep = (sHigh - sLow)/sStepNum;
Table[{s, Chop[pdf[s]]}, {s, sLow, sHigh, sStep}]]

pdf0 = Interpolation[pdf0Table, InterpolationOrder -> 1]

Plot[pdf0[s], {s, 0.001, 3.}, PlotRange -> Full]

NIntegrate[pdf0[s], {s, 0.001, 3.}]

FindMaximum[pdf0[s], {s, 0.5}]


My weakness (and the slowness of computing points in pdf[s]) prevents me taking this further, sorry. 