Maximum likelihood estimation based on frequency table

Question

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[1096]:= ClearAll[likelihoodf]
likelihoodf[m_?NumericQ, std_?NumericQ] := 
 NProbability[70 < x < 80, 
  x \[Distributed] 
   TruncatedDistribution[{65, 100}, NormalDistribution[m, std]]]

In[1098]:= likelihoodf[75, 5]

Out[1098]= 0.698583

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

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

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

During evaluation of In[1102]:= 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[1102]:= NormalDistribution::posprm: Parameter -0.13703 at position 2 in NormalDistribution[-0.535769,-0.13703] is expected to be positive. >>

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

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

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

During evaluation of In[1102]:= 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[1102]= NArgMax[likelihoodf[m, std], {m, std}]

And the following method also returns an error:

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

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

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

During evaluation of In[1101]:= 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[1101]= {0.887236, {m -> 74., std -> 3.}}

How can I maximize such a likelihood function?

Answer 1

You can get the function explicitly :

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

Answer 2

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

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

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

Out[1118]= 0.698583

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

Out[1119]= {75., 1.}