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?

Improve this question
$\endgroup$
8
  • 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
    Commented Feb 13, 2015 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
    Commented Feb 13, 2015 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
    Commented Feb 13, 2015 at 22:23
  • $\begingroup$ a PDF should be non-negative everywhere $\endgroup$
    – Manuel --Moe-- G
    Commented Feb 13, 2015 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
    Commented Feb 13, 2015 at 22:25

2 Answers 2

Reset to default
6
$\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}}

*)
Improve this answer
$\endgroup$
2
  • $\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
    Commented Feb 14, 2015 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
    Commented Feb 14, 2015 at 22:46
2
$\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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.