FindDistributionParameters gives an error for a mixture of user defined sinh-arcsinh distributions

Question

I have some data (1593 values). I first binned it using the Knuth rule and fitted a function that is a linear combination of two sinh-arcsinh functions. I then rounded the resulting parameters to perform a FindDistributionParameters search of parameters using "MaximumLikelihood" as a ParameterEstimator. The problem is that after a lengthy computation I get an error message:

FindMaximum::eit: "The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 100 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {8.33678,100.,4.16864}, is returned.

Unfortunately, there's no such thing like MaxIterations->Infinity in FindDistributionParameters.

The suitable code below; file.txt available here

import = Import["file.txt", "Table"];
data = Flatten[import];
data = Log[10, data];
data = Delete[data, Position[data, RankedMin[data, 1]][[1, 1]]];
data = Delete[data, Position[data, RankedMin[data, 1]][[1, 1]]];
data = Delete[data, Position[data, RankedMax[data, 1]][[1, 1]]];

sas[μ_, σ_, skew_, kurt_] := 
ProbabilityDistribution[((1 + ((z - μ)/σ)^2)^(-(1/2))
kurt Cosh[
kurt ArcSinh[(z - μ)/σ] - 
skew] Exp[-(1/2) Sinh[
kurt ArcSinh[(z - μ)/σ] - skew]^2])/(
Sqrt[2 π] σ), {z, -Infinity, Infinity}]

dist = MixtureDistribution[{p, 
1 - p}, {sas[μ1, σ1, skew1, kurt1], 
sas[μ2, σ2, skew2, kurt2]}]

param = FindDistributionParameters[data, 
dist, {{μ1, -0.2}, {σ1, 0.7}, {skew1, 0.3}, {kurt1, 
1.}, {μ2, 1.5}, {σ2, 0.5}, {skew2, 0.}, {kurt2, 
1.}, {p, 0.25}}, 
ParameterEstimator -> {"MaximumLikelihood", 
Method -> "NMaximize"}]

I'm not sure whether the error may be ignored, or it implies that there is room for a significant improvement of the obtained result. If the error may be dealt with, what's the solution?

Answer 1

I rewrote quite a few things... manually specifying the mixture, and setting very low iteration and precision goals because it hangs for what seems like an eternity. Importantly I also specified assumptions on the distribution:

sas[μ_, σ_, skew_, kurt_, 
  z_] := ((1 + ((z - μ)/σ)^2)^(-(1/2)) kurt Cosh[
     kurt ArcSinh[(z - μ)/σ] - 
      skew] Exp[-(1/2) Sinh[
        kurt ArcSinh[(z - μ)/σ] - skew]^2])/(Sqrt[
     2 π] σ)

(*dist=MixtureDistribution[{p,1-p},{sas[μ1,σ1,skew1,kurt1],\
sas[μ2,σ2,skew2,kurt2]}]*)
dist[p_, μ1_, μ2_, σ1_, σ2_, skew1_, skew2_, 
  kurt1_, kurt2_] := 
 ProbabilityDistribution[
  p*sas[μ1, σ1, skew1, kurt1, z] + (1 - p)*
    sas[μ2, σ2, skew2, kurt2, z], {z, -Infinity, Infinity},
   Assumptions -> {σ1 > 0, σ2 > 0, skew1 >= 0, 
    skew2 >= 0, kurt1 >= 0, kurt2 >= 0, 0 <= p <= 1}]


param = FindDistributionParameters[data, 
  dist[p, μ1, μ2, σ1, σ2, skew1, skew2, kurt1, 
   kurt2], {{μ1, -0.2}, {σ1, 0.7}, {skew1, 0.3}, {kurt1, 
    1.}, {μ2, 1.5}, {σ2, 0.5}, {skew2, 0.}, {kurt2, 
    1.}, {p, 0.25}}, 
  ParameterEstimator -> {"MaximumLikelihood", 
    Method -> {"NMaximize", PrecisionGoal -> 1, MaxIterations -> 5}}]
Show[Histogram[data, Automatic, "PDF"], 
 Plot[PDF[dist[p, μ1, μ2, σ1, σ2, skew1, skew2, 
     kurt1, kurt2] /. param, 
   z], {z, -2, 3}], PlotRange -> All]

Answer 2

If you only need an estimate of the density function as opposed to estimates and standard errors for the coefficients of a specified model structure, then obtaining a kernel density estimate might be all you need.

(* Define distribution *)
sas[μ_, σ_, skew_, kurt_] := ProbabilityDistribution[
  ((1 + ((z - μ)/σ)^2)^(-(1/2)) kurt Cosh[
      kurt ArcSinh[(z - μ)/σ] - skew]
      Exp[-(1/2) Sinh[
         kurt ArcSinh[(z - μ)/σ] - skew]^2])/(Sqrt[
      2 π] σ),
  {z, -Infinity, Infinity}]

(* Mixture distibution *)
dist = MixtureDistribution[{p, 1 - p}, {sas[μ1, σ1, skew1, kurt1], 
    sas[μ2, σ2, skew2, kurt2]}];

(* Generate some data *)
data = RandomVariate[
   dist /. {p -> 0.25, μ1 -> -0.2, σ1 -> 0.7, skew1 -> 0.3,
      kurt1 -> 1, μ2 -> 1.5, σ2 -> 0.5, skew2 -> 0, kurt2 -> 1}, 1593];

(* Histogram *)
p1 = Histogram[data, Automatic, "PDF"];

(* Nonparametric density estimate *)
d = SmoothKernelDistribution[data];
p2 = Plot[PDF[d, x], {x, -2, 3.2}];

(* Show both figures *)
Show[{p1, p2}]

This is not an answer to the specific question you have but does provide an alternative which avoids the error you see (not to mention avoiding the lengthy computation time).