# Maximum likelihood estimation based on frequency table

I have a frequency table of grades, I'm trying to estimate the distribution as a truncated normal distribution (a grade below 65 is a failing grade and doesn't count).

Since I don't have the actual "data" samples needed for many Mathematica built-in functions, I have to define my own likelihood function.

Assuming I have $p_i$ grades in the $i$th bin - $[a_i,b_i]$, the likelihood is defined as:

$$L(\mu,\sigma) = \prod_{i} Prob(a_i<x<b_i | 65<x<100)^{p_i}$$ Where $X \sim N(\mu,\sigma)$.

To simplify the question, here's the likelihood of one grade to be in the $[70,80]$ bin, for certain $\mu,\sigma$:

In:= ClearAll[likelihoodf]
likelihoodf[m_?NumericQ, std_?NumericQ] :=
NProbability[70 < x < 80,
x \[Distributed]
TruncatedDistribution[{65, 100}, NormalDistribution[m, std]]]

In:= likelihoodf[75, 5]

Out= 0.698583


But when I'm trying to maximize it I'm receiving an error:

In:= NArgMax[likelihoodf[m, std], {m, std}]

During evaluation of In:= Power::infy: Infinite expression 1/0. encountered. >>

During evaluation of In:= NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>

During evaluation of In:= NormalDistribution::posprm: Parameter -0.13703 at position 2 in NormalDistribution[-0.535769,-0.13703] is expected to be positive. >>

During evaluation of In:= NormalDistribution::posprm: Parameter -0.13703 at position 2 in NormalDistribution[-0.535769,-0.13703] is expected to be positive. >>

During evaluation of In:= NormalDistribution::posprm: Parameter -0.13703 at position 2 in NormalDistribution[-0.535769,-0.13703] is expected to be positive. >>

During evaluation of In:= General::stop: Further output of NormalDistribution::posprm will be suppressed during this calculation. >>

During evaluation of In:= NArgMax::nnum: The function value -NProbability[70<x<80,x\[Distributed]TruncatedDistribution[{65,100},NormalDistribution[-0.535769,-0.13703]]] is not a number at {m,std} = {-0.535769,-0.13703}. >>

Out= NArgMax[likelihoodf[m, std], {m, std}]


And the following method also returns an error:

In:= FindMaximum[likelihoodf[m, std], {{m, 74}, {std, 3}}]

During evaluation of In:= NormalDistribution::posprm: Parameter -0.891788 at position 2 in NormalDistribution[75.3113,-0.891788] is expected to be positive. >>

During evaluation of In:= NormalDistribution::posprm: Parameter -0.891788 at position 2 in NormalDistribution[75.3113,-0.891788] is expected to be positive. >>

During evaluation of In:= FindMaximum::nrnum: The function value -NProbability[70<x<80,x\[Distributed]TruncatedDistribution[{65,100},NormalDistribution[75.3113,-0.891788]]] is not a real number at {m,std} = {75.3113,-0.891788}. >>

Out= {0.887236, {m -> 74., std -> 3.}}


How can I maximize such a likelihood function?

You can get the function explicitly :

Probability[70 < x < 80, x \[Distributed] TruncatedDistribution[{65, 100}, NormalDistribution[m, std]]]
(* (Erfc[(-80 + m)/(Sqrt std)] - Erfc[(-70 + m)/(Sqrt std)])/(2 (1/2 Erfc[(-100 + m)/(Sqrt std)] - 1/2 Erfc[(-65 + m)/(Sqrt std)])) *)


I think most of the numerical errors are due to the denominator becoming very small, so I am using an artificial cutoff :

cutoff at 10^-3:

NMaximize[{(Erfc[(-80 + m)/(Sqrt std)] - Erfc[(-70 + m)/(Sqrt std)])/(2 (1/2 Erfc[(-100 + m)/(Sqrt std)] -
1/2 Erfc[(-65 + m)/(Sqrt std)]) + 10^-3), 0 < std(*, 50<m<
100*)}, {m, std}, Reals]
(* {0.9995, {m -> 71.7539, std -> 0.000370381}} *)


cutoff at 10^-9:

NMaximize[{(
Erfc[(-80 + m)/(Sqrt std)] - Erfc[(-70 + m)/(Sqrt std)])/(
2 (1/2 Erfc[(-100 + m)/(Sqrt std)] -
1/2 Erfc[(-65 + m)/(Sqrt std)]) + 10^-9), 0 < std(*, 50<m<
100*)}, {m, std}, Reals]
(* {1., {m -> 71.7539, std -> 0.000370381}} *)


I also tried to play with WorkingPrecision and the like in order to avoid using the cutoff but I get messages about not enough memory.

• Thanks a lot for your tips! I've posted another answer. I think constrained optimization is better for the problem. But your response was really helpful! – User Jan 30 '15 at 11:39

Using b.gatessucks tips, I managed to get a better solution. I'm using constrained optimization with some reasonable limits on the parameters:

In:= ClearAll[likelihoodf]
likelihoodf[m_?NumericQ, std_?NumericQ] :=
Evaluate[Probability[70 < x < 80,
x \[Distributed]
TruncatedDistribution[{65, 100}, NormalDistribution[m, std]]]]

In:= likelihoodf[75, 5] // N

Out= 0.698583

In:= NArgMax[{likelihoodf[m, std],
65 < m < 100 \[And] 1 < std < 50}, {m, std}]

Out= {75., 1.}