# FindDistributionParameters fails with custom distribution which cannot be evaluated symbolically

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}] 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] works fine. I can of course generate a sample from this customized PDF

dat = RandomVariate[myDistribution[2/10], 500];
dat // Histogram 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. >> *) 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][] // 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.

• 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. – Sjoerd Smit Dec 22 '16 at 15:27

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).