1
$\begingroup$

I am trying to estimate the parameters of a custom distribution (an approximation of the skew normal), but I can't get Mathematica to use DistributionFitTest, EstimatedDistribution, or FindDistributionParameters on any custom distribution or a distribution where I have applied a constraint like

SkewNormalDistribution[μ, σ, (Sqrt[π] μ)/Sqrt[-π μ^2 + 2 σ^2]]

What I have currently is

data = RandomVariate[SkewNormalDistribution[-.003, .01, 3], 150];

SNAPDF[x_, α_] := Piecewise[{
{0, x < -3/α},
{1/(8 Sqrt[2*Pi])Exp[-x^2/2] (9*α*x + 3*α^2 x^2 + 
    1/3 α^3 x^3 + 9), -3/α <= x < -1/α}, 
{1/(4 Sqrt[2*Pi])Exp[-x^2/2] (3*α*x - 1/3 α^3 x^3 + 
    4), -1/α <= x < 1/α},
{1/(8 Sqrt[2*Pi])Exp[-x^2/2] (9*α*x - 3*α^2 x^2 + 
    1/3 α^3 x^3 + 7), 1/α <= x < 3/α}, 
{Sqrt[2/Pi] Exp[-x^2/2],3/α <= x}
}];

SNA[μ_, σ_, α_] := 
  ProbabilityDistribution[
   SNAPDF[(x - μ)/σ, α]/σ, {x, -Infinity, Infinity}, 
   Assumptions -> μ ∈ Reals && σ > 
      0 && α ∈ Reals];

EstimatedDistribution[
 data, 
 SNA[μ, σ, α], {{μ, -.01}, {σ, .01}, {α, 3}}]

I've looked through quite a few similar posts about topics like this but haven't found anything to solve my specific problems. Thanks!

$\endgroup$
6
  • $\begingroup$ Add code fill data variable. $\endgroup$
    – Edmund
    Aug 12, 2021 at 19:54
  • $\begingroup$ Do you not need $\frac{\text{SNAPDF}\left(\frac{x-\mu }{\sigma },\alpha \right)}{\sigma }$ rather than just $\text{SNAPDF}\left(\frac{x-\mu }{\sigma },\alpha \right)$ in the definition of SNA ? $\endgroup$
    – JimB
    Aug 13, 2021 at 5:45
  • $\begingroup$ Good point! I have now fixed that in the post, but both of the problems still persist $\endgroup$
    – John Smith
    Aug 13, 2021 at 15:03
  • $\begingroup$ I wonder if you should remove the reference to the constraints as it doesn't appear to have anything to do with estimating the parameters. Also, with $\mu=-0.01$ and $\sigma=0.01$, the constraint has $\alpha=\frac{\sqrt{\pi } \mu }{\sqrt{2 \sigma ^2-\pi \mu ^2}}$ being $\alpha=i \sqrt{\frac{\pi }{\pi -2}}$ (an imaginary number). $\endgroup$
    – JimB
    Aug 13, 2021 at 17:07
  • $\begingroup$ I will edit the parameters inside the RandomVariate such that α remains real in the context of applying the constraint to the sample data! The end goal is to try to fit a constrained version of the SNA distribution to the data, but Mathematica won't fit a constrained nor user defined distribution. $\endgroup$
    – John Smith
    Aug 13, 2021 at 19:34

1 Answer 1

1
$\begingroup$

I understand that you need to be able to use constraints but I don't see how that is related to the estimation problem you present. But what you present is a real problem.

To obtain the maximum likelihood estimates I've slightly modified your code. I think that your complicated Piecewise function where the conditions are highly dependent on the parameters (as opposed to conditions not involving the parameters) is bogging down the calculations in EstimatedDistribution. (At least that's my guess.)

(* Generate data *)
SeedRandom[12345];
data = RandomVariate[SkewNormalDistribution[-.01, .01, 3], 150];

(* Initial values for parameters based on skew normal distribution *)
{μ0, σ0, α0} = {μ, σ, α} /. FindDistributionParameters[data, SkewNormalDistribution[μ, σ, α]];

(* Define log of the likelihood *)
logL[data_, μ_?NumericQ, σ_?NumericQ, α_?NumericQ] := Module[{f, x},
  f = Piecewise[{{0, x < -3/α},
    {1/(8 Sqrt[2*Pi]) Exp[-x^2/2] (9*α*x + 3*α^2 x^2 + 1/3 α^3 x^3 + 9), -3/α <= x < -1/α},
    {1/(4 Sqrt[2*Pi]) Exp[-x^2/2] (3*α*x - 1/3 α^3 x^3 + 4), -1/α <= x < 1/α},
    {1/(8 Sqrt[2*Pi]) Exp[-x^2/2] (9*α*x - 3*α^2 x^2 + 1/3 α^3 x^3 + 7), 1/α <= x < 3/α},
    {Sqrt[2/Pi] Exp[-x^2/2], 3/α <= x}}]/σ /. x -> (x - μ)/σ // Simplify;
  Log[f /. x -> #] & /@ data // Total
  ]

(* Maximize the likelihood *)
mle = FindMaximum[{logL[data, μ, σ, α], (Min[data] - μ)/σ >= -3/α, σ > 0}, 
  {{μ, μ0}, {σ, σ0}, {α, α0}}]
(* {554.03, {μ -> -0.0101752, σ -> 0.00993595, α -> 3.62204}} *)

(* Display data histogram and fit *)
Show[Histogram[data, "FreedmanDiaconis", "PDF"],
 Plot[SNAPDF[(x - μ)/σ, α]/σ /. mle[[2]], {x, Min[data], Max[data]}]]

Histogram and fit

I'm also not understanding the purpose of performing a goodness-of-fit test. You know for a fact that your Piecewise function is not a skewed normal. No need to test for that. What you would seem to want is a goodness-of-fit summary statistic or two or three summary statistics (certainly something other than a P-value). In other words, what is it that would say that your estimation of the approximate density is close enough to what you can get by using a skewed normal?

Obtaining estimates of standard errors for the parameter estimates.

f = Piecewise[{{0, x < -3/α}, 
  {1/(8 Sqrt[2*Pi]) Exp[-x^2/2] (9*α*x + 3*α^2 x^2 + 1/3 α^3 x^3 + 9), -3/α <= x < -1/α}, 
  {1/(4 Sqrt[2*Pi]) Exp[-x^2/2] (3*α*x - 1/3 α^3 x^3 + 4), -1/α <= x < 1/α}, 
  {1/(8 Sqrt[2*Pi]) Exp[-x^2/2] (9*α*x - 3*α^2 x^2 + 1/3 α^3 x^3 + 7), 1/α <= x < 3/α}, 
  {Sqrt[2/Pi] Exp[-x^2/2], 3/α <= x}}]/σ /. x -> (x - μ)/σ // Simplify;

(covMat = -Inverse[((D[Log[f], {{μ, σ, α}, 2}] /. 
  mle[[2]]) /. x -> #) & /@ data // Total]) // MatrixForm

Covariance matrix

se = Sqrt[Diagonal[covMat]]
(* {0.000698432, 0.000786863, 0.930417} *)
$\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.