# Why the difference in these results of estimating the parameters of Erlang distribution?

Investigating (in 13.2 on Windows 10) the FindDistributionParameters command for learning purpose, I consider a small sample {1.1, 1.9, 2.1, 3.1, 3.9} from a possible statistical population distributed as ErlangDistribution. I estimate the parameters, making use of different methods.

FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9},
ErlangDistribution[k, \[Lambda]], ParameterEstimator -> "MaximumLikelihood"]


{k -> 6, \[Lambda] -> 2.47934}

The above is identical to the result of

NMaximize[{LogLikelihood[ErlangDistribution[k, \[Lambda]], {1.1, 1.9, 2.1, 3.1, 3.9}],
k \[Element] PositiveIntegers && \[Lambda] >= 0}, {k, \[Lambda]},
Method -> {"DifferentialEvolution", "ScalingFactor" -> 1}]


and

FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9},  ErlangDistribution[k, \[Lambda]],
ParameterEstimator -> {"MaximumLikelihood", Method -> "FindMaximum"}]


clearly says that k must be a positive integer number.

Then

FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9},  ErlangDistribution[k, \[Lambda]],
ParameterEstimator -> "MethodOfMoments"]


{k -> 4.91309, \[Lambda] -> 2.0302}

and

FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9},  ErlangDistribution[k, \[Lambda]],
ParameterEstimator -> "MethodOfCentralMoments"]


produces the same. The value of k is not integer, but the documentation says "Method-of-moment-based estimators may not satisfy all restrictions on parameters".

Let us write down

eq1 = Moment[{1.1, 1.9, 2.1, 3.1, 3.9}, 1] == Moment[ErlangDistribution[k, \[Lambda]], 1]
eq2 = Moment[{1.1, 1.9, 2.1, 3.1, 3.9}, 2] == Moment[ErlangDistribution[k, \[Lambda]], 2]


,equating sample moments and distribution moments, and solve the system

Solve[{eq1, eq2}, {k, \[Lambda]}, Reals]


{{k -> 6.14136, \[Lambda] -> 2.53775}}

We see the difference in the values: {k -> 4.91309, \[Lambda] -> 2.0302} vs {k -> 6.14136, \[Lambda] -> 2.53775}. How to explain it?

PS. [CASE:5013654]

• FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9}, ErlangDistribution[k, \[Lambda]], ParameterEstimator -> {"MethodOfMoments", "MomentOrders" -> {1, 3}}] performs {k -> 6.56024, \[Lambda] -> 2.71084}. Commented Mar 13, 2023 at 20:07
• Solve[CentralMoment[{1.1, 1.9, 2.1, 3.1, 3.9}, 1] == CentralMoment[NormalDistribution[\[Mu], \[Sigma]], 1] && CentralMoment[{1.1, 1.9, 2.1, 3.1, 3.9}, 2] == CentralMoment[NormalDistribution[\[Mu], \[Sigma]], 2], {\[Mu], \[Sigma]}, Reals] produces {{\[Sigma]->-0.976524},{\[Sigma]->0.976524}} and Solve[CentralMoment[{1.1, 1.9, 2.1, 3.1, 3.9}, 3] == CentralMoment[NormalDistribution[\[Mu], \[Sigma]], 3] && CentralMoment[{1.1, 1.9, 2.1, 3.1, 3.9}, 2] == CentralMoment[NormalDistribution[\[Mu], \[Sigma]], 2], {\[Mu], \[Sigma]}, Reals] produces {}. Commented Mar 14, 2023 at 5:23
• In the above we see that defect occurs for other distributions too. Commented Mar 14, 2023 at 5:24
• The problem is that you don't understand that the first central moment for all sample data and distributions with a mean is zero and therefore does not contain any information for estimating the parameters. Now it certainly doesn't help that for two-parameter distributions Mathematica uses the first raw moment (the mean) and the second central moment when one chooses "MethodOfCentralMoments" if one doesn't specify the moment orders to use.
– JimB
Commented Mar 14, 2023 at 6:10
• Even much worse the result of Solve[Moment[{1.1, 1.9, 2.1, 3.1, 3.9}, 1] == Moment[TriangularDistribution[{min, max}, c], 1] && Moment[{1.1, 1.9, 2.1, 3.1, 3.9}, 2] == Moment[TriangularDistribution[{min, max}, c], 2], {min, max, c}, Reals] , though FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9}, TriangularDistribution[{min, max}, c], ParameterEstimator -> "MethodOfMoments"] produces {min -> 0.246228, max -> 4.98257, c -> 2.0312}. Commented Mar 14, 2023 at 6:14

I can explain what causes the difference but not why.

Your use of the raw moments can be rewritten in the following manner:

data = {1.1, 1.9, 2.1, 3.1, 3.9};
Solve[{Mean[data] == Mean[ErlangDistribution[k, λ]],
data.data/Length[data] == Moment[ErlangDistribution[k, λ], 2]}]
(* {{k -> 6.14136, \[Lambda] -> 2.53775}} *)


But what FindDistributionParameters does is essentially one or both of the following:

Solve[{Mean[data] == Mean[ErlangDistribution[k, λ]],
Variance[data] == Variance[ErlangDistribution[k, λ]]}, {k, λ}]
(* {{k -> 4.91309, λ -> 2.0302}} *)

Solve[{Mean[data] == Mean[ErlangDistribution[k, λ]],
Variance[data] + Mean[data]^2 == Moment[ErlangDistribution[k, λ], 2]}]
(* {{k -> 4.91309, λ -> 2.0302}} *)


The differences in the approaches will get smaller with larger sample sizes.

• Thank for your explanation. It seems that FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9}, ErlangDistribution[k, \[Lambda]], ParameterEstimator -> "MethodOfCentralMoments"] returning the same {k -> 4.91309, λ -> 2.0302} also does not use any method of central moments. Commented Mar 14, 2023 at 4:37
• Also FindDistributionParameters[{1.1, 1.9, 2.1, 3.1, 3.9}, ErlangDistribution[k, \[Lambda]], ParameterEstimator -> "MethodOfCumulants"] returns the same {k -> 4.91309, \[Lambda] -> 2.0302}. Commented Mar 14, 2023 at 4:40
• Therefore, the documentation to FindDistributionParameters does not correspond to reality. I find it dangerous. Commented Mar 14, 2023 at 4:48
• The documentation gives the standard warning that you posted "Method-of-moment-based estimators may not satisfy all restrictions on parameters". The use of moments also doesn't always use all of the available information to estimate parameters. The warnings are there.
– JimB
Commented Mar 14, 2023 at 6:00
• If one just uses the first two cumulants, those correspond exactly to the mean and variance, respectively. It would be only if more or other cumulants are selected, would "MethodOfCumulants" provide a different set of estimates.
– JimB
Commented Mar 14, 2023 at 6:03

If k is constrained to be a positive integer, "MethodOfMoments" is equivalent to "MaximumLikelihood"

Clear["Global*"]

data = {1.1, 1.9, 2.1, 3.1, 3.9};

dist = ErlangDistribution[k, λ];

param1 = FindDistributionParameters[data, dist]

(* {k -> 6, λ -> 2.47934} *)

DistributionParameterAssumptions[dist]

(* k ∈ Integers && k > 0 && λ > 0 *)


Constrain k to be a positive integer

Off[Solve::ratnz]

kVal = Round[SolveValues[
Moment[dist, #] == Moment[data, #] & /@ {1, 2},
k, {λ}][[1]]]

(* 6 *)

param2 = Prepend[
FindDistributionParameters[data,
ErlangDistribution[kVal, λ],
ParameterEstimator -> "MethodOfMoments"],
k -> kVal]

(* {k -> 6, λ -> 2.47934} *)

param3 = Prepend[Solve[
{Moment[ErlangDistribution[kVal, λ], 1] ==
Moment[data, 1], λ > 0}, λ][[1]], k -> kVal]

(* {k -> 6, λ -> 2.47934} *)

param1 == param2 == param3

(* True *)

• Thank you for your work. Unfortunately, that you use is not any method of moments. In particular, it is known that rounding may not produce an optimal solution. Therefore, your claim "If k is constrained to be a positive integer, "MethodOfMoments" is equivalent to "MaximumLikelihood"" does not correspond to reality. Commented Mar 14, 2023 at 5:14
• If k is not an integer then the distribution is not a valid ErlangDistribution and does not "correspond to reality". To enforce an integer parameter, an integer function is required (e.g., Round, Ceiling, or Floor) with the method of moments used to determine k. And both param2 and param3 continued to use a method of moments to find the remaining parameter. Commented Mar 14, 2023 at 5:42
• As @JimB demonstrates in his answer, in fact, ParameterEstimator -> "MethodOfMoments"` likely does not use any method of moments, but you apply this option. This is a point. Commented Mar 14, 2023 at 5:55
• Just a note that the Erlang distribution is a special case of the gamma distribution with the restriction that $k$ is an integer. When $k$ is not an integer, plugging that into the Erlang pdf formula results in a legitimate gamma pdf.
– JimB
Commented Mar 14, 2023 at 6:30