6
$\begingroup$

Context

I am interested in measuring the dark energy equation of state of the universe while fitting the PDF of density in cells. This involves fitting a one parameter PDF to some data.


Following this question I would therefore like to find the MaximumLikelihood solution of a customized PDF

I have a user defined one parameter PDF ρPDFns (given below, or described scientifically here) which I can evaluate for any value of its parameter σ

pl = Plot[
   Table[ρPDFns[ρ, σ, -3/2], {σ, 0.1, 0.4, 0.05}] // Release, {ρ,
     0.1, 4}, PlotRange -> {0., 5}]

Mathematica graphics

I would like to estimate σ using a MaximumLikelihood approach.

I follow the documentation and define

Clear[myDistribution]; 
myDistribution[σ_] := 
 ProbabilityDistribution[ρPDFns[ρ, σ, -3/2], {ρ,
    1/10, 4}, Assumptions -> (σ > 0 && σ < 5/10)];

For instance

Plot[PDF[myDistribution[2/10], x], {x, 1/10, 4}, PlotRange -> All]

Mathematica graphics

works fine. I can of course generate a sample from this customized PDF

dat = RandomVariate[myDistribution[2/10], 500];
dat // Histogram

Mathematica graphics

The problem

Maximumlikelihood estimation fails:

FindDistributionParameters[dat, 
 myDistribution[σ], {σ, 1/10},
 ParameterEstimator -> {"MaximumLikelihood", Method -> "NMaximize", 
   PrecisionGoal -> 1, MaxIterations -> 5}]

It produces this kind of error messages.

(* FindRoot::srect: Value [Sigma] in search specification {s$23324,[Sigma]} is not a number or array of numbers. >> *)

Mathematica graphics

Question

How can I make FindDistributionParameters work on one parameter PDF which cannot be evaluated before the parameter has a numerical value.


Here is the code for the PDF.

ζ[x_, ν_] = 1/(1 - x/ν)^ν;
τ[y_, ν_] =  Module[{x}, x /. Solve[y == ζ[x, ν], x][[1]] // Quiet];
σ[r_, ns_, α_] = Sqrt[2 *1/((r)^(3 + ns + α) + (r)^(3 +   ns - α))];
ψ[y_, ν_, ns_, α_] = 1/2 τ[y, ν]^2/σ[y^(1/3), ns, α]^2;
dψ[y_, ν_, ns_, α_] = D[ψ[y, ν, ns, α], {y, 1}];
ddψ[y_, ν_: 21/13, ns_, α_] =  D[ψ[y, ν, ns, α], {y, 2}];
Lowlogrho[x_, s_, ν_, ns_, α_] =  1/Sqrt[2 Pi]*1/s*
   Sqrt[ddψ[x, ν, ns, α] + 
dψ[x, ν, ns, α]/x] Exp[-1/s^2*ψ[x, ν,   ns, α]];
Clear[logs0, logs1, logs2];
logs0[s_, ν_, ns_, α_] :=  logs0[s, ν, ns, α] = 
   Module[{β},  NIntegrate[  Lowlogrho[β, s, ν, ns, α], {β, 0.1,   10}]];
logs1[s_, ν_, ns_, α_] :=  logs1[s, ν, ns, α] = 
   Module[{β},     NIntegrate[ β Lowlogrho[β, s, ν, 
       ns, α], {β, 0.1, 10}]];
logs2[s_, ν_, ns_, α_] :=   logs2[s, ν, ns, α] = 
  Module[{β},  NIntegrate[ β^2 Lowlogrho[β, s, ν, 
      ns, α], {β, 0.1, 10}]]; Clear[effsiglog]; 
effsiglog[sig_, ν_, ns_, α_] := effsiglog[sig, ν, ns, α] = 
  Module[{s},  s /. FindRoot[logs0[s, ν, ns, α] logs2[s, ν, ns, α]/
         logs1[s, ν, ns, α]^2 == 1 + sig^2, {s, sig}] //Quiet];
Options[ρPDFns] = {ν -> 21/13, α -> 0};
Clear[ρPDFns]; ρPDFns[ρ_, sigma_, ns_, OptionsPattern[]] :=
 logs1[effsiglog[sigma, OptionValue[ν], ns, 
      OptionValue[α]], OptionValue[ν], ns, 
     OptionValue[α]]/
    logs0[effsiglog[sigma, OptionValue[ν], ns, 
       OptionValue[α]], OptionValue[ν], ns, 
      OptionValue[α]]^2 Lowlogrho[ρ*
     logs1[effsiglog[sigma, OptionValue[ν], ns, 
        OptionValue[α]], OptionValue[ν], ns, 
       OptionValue[α]]/
      logs0[effsiglog[sigma, OptionValue[ν], ns, 
        OptionValue[α]], OptionValue[ν], ns, 
       OptionValue[α]], 
    effsiglog[sigma, OptionValue[ν], ns, OptionValue[α]], 
    OptionValue[ν], ns, OptionValue[α]] // Quiet;

PS: a simpler, in my opinion equivalent test problem, is given in this question.

$\endgroup$
  • $\begingroup$ I'm fairly sure that part of the problem is that part of the problem here is that during the processing of the problem, Mathematica is trying to symbolically evaluate functions that should only evaluate with numerical inputs. For example, your function effsiglog has a FindRoot in it, so it should never try to evaluate symbolically. To prevent this from happening, you can use the _?NumericQ filter pattern on your inputs: effsiglog[sig__?NumericQ, ν__?NumericQ, ns__?NumericQ, α__?NumericQ] :=... This alone may not solve the problem, but it might at least be a step in the right direction. $\endgroup$ – Sjoerd Smit Dec 22 '16 at 15:27
2
$\begingroup$

I can't directly answer your question

How can I make FindDistributionParameters work on a one parameter PDF which cannot be evaluated before the parameter has a numerical value?

but the maximum likelihood estimate can be obtained by first following @SjoerdSmit 's advice about using ?NumericQ in the definition of the distribution:

myDistribution[σ_?NumericQ] := ProbabilityDistribution[ρPDFns[ρ, σ, -3/2],
  {ρ, 1/10, 4}, Assumptions -> (σ > 0 && σ < 5/10)]

followed by the use of the LogLikelihood function:

dat = RandomVariate[myDistribution[2/10], 500];
FindMaximum[{LogLikelihood[myDistribution[σ], dat], 0 < σ < 5/10}, {σ, 0.1}]
(* {124.317, {σ -> 0.198984}} *)

I don't know why FindDistributionParameters fails (even with the use of ?NumericQ).

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