Estimating parameters of a custom or constrained distribution

Question

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!

Answer 1

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

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

se = Sqrt[Diagonal[covMat]]
(* {0.000698432, 0.000786863, 0.930417} *)