2
$\begingroup$

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]] - 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?

$\endgroup$
  • 1
    $\begingroup$ have you tried calculating a few dozen points of pdf[s] and generating substituting a Interpolation function to replace use of pdf[s]? $\endgroup$ – Manuel --Moe-- G Feb 13 '15 at 21:45
  • $\begingroup$ 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 $\endgroup$ – Remi.b Feb 13 '15 at 22:02
  • $\begingroup$ 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}} $\endgroup$ – Manuel --Moe-- G Feb 13 '15 at 22:23
  • $\begingroup$ a PDF should be non-negative everywhere $\endgroup$ – Manuel --Moe-- G Feb 13 '15 at 22:23
  • $\begingroup$ 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}] ] $\endgroup$ – Manuel --Moe-- G Feb 13 '15 at 22:25
4
$\begingroup$

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}}

*)
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Remi.b Feb 14 '15 at 22:25
  • $\begingroup$ @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 ;-) ). $\endgroup$ – ciao Feb 14 '15 at 22:46
1
$\begingroup$

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.

linear interpolation of pdf

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