1
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – User Jan 30 '15 at 11:39
1
$\begingroup$

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.}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.