When analytical and numerical methods do not agree - Case study with Maximum Likelihoods methods

Question

Here is the probability distribution I am interested in:

$$P(q)=C e^{4 n s q} q^{4 n \nu - 1} (1 - q)^{4 n \mu - 1}$$

, where $e$ is the constant of Euler and $C$ is constant so that the whole thing integrates to 1. $n$, $\mu$, $\nu$ and $s$ are four parameters of the distribution. My goal is to use a maximum likelihood method in order to estimate the most likely value of the parameter $s$ given the other parameters and given the variable $q$. But I am encountering some issues related with different results obtained when using analytical versus numerical methods. Can you help me to solve this issue?

---

Below is what I tried...

I first tried to find the value for $C$ so I did:

1/(Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
     4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q) q^(
   4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1) /. s -> \[Mu]/10 /. 
 q -> 0.75

and got

1/(Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
         4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q) q^(
       4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1) /. s -> \[Mu]/10 /. 
     q -> 0.75

Given that conditions are always respected for the range of parameters I want to consider, I am satisfy with this answer. As a quick check, I made sure that it indeed integrates to one for a realistc choice of parameters by running.

\[Mu] = 10^-6
\[Nu] = \[Mu]
n = 100/(40 \[Nu])
Integrate[
 Re[1/(Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
      4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q) q^(
    4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1) /. s -> \[Mu]/10], {q, 0, 
  1}] 

And indeed it integrates to 1. However, if I use NIntegrate...

NIntegrate[
     Re[1/(Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
          4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q) q^(
        4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1) /. s -> \[Mu]/10], {q, 0, 
      1}]

it does not at all integrate to one (it integrates to 0.0000538285). If I calculate for the probability for a given $q$

1/(Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
     4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q) q^(
   4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1) /. s -> \[Mu]/10 /. 
 q -> 0.75

I get some tiny probability (0.0000181795 in this case). This probability corresponds to what can be seen when we graph the probability distribution:

Plot[1/(Gamma[4 n \[Mu]] Gamma[
     4 n \[Nu]] Hypergeometric1F1Regularized[4 n \[Nu], 
     4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q ) q^(
   4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1) /. s -> \[Mu]/10, {q, 0, 1}]

Looking at this graph, it is obvious that it actually does not integrate to 1.

Below is my attempt to visualize the log-likelihood curve (for each $s$ considered)

searchS = Table[s, {s, 0, 1, 0.001}];
likelihoods = List[];
Do[s = searchS[[i]];
 likelihoods = 
  Append[likelihoods, 
   1/(Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
      4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s]) E^(4 n s q ) q^(
    4 n \[Nu] - 1) (1 - q)^(4 n \[Mu] - 1)]
 , {i, 1, Length[searchS]}]
loglikelihoods = Log[likelihoods];
ListPlot[Transpose[{searchS, loglikelihoods}], 
 AxesLabel -> {"s", "log(likelihood)"}]

which is definitely not what I expected.

Answer 1

The problem is numerically unstable for some parameter ranges.

We shall show a simple example.

Your normalized distribution is given by

p[q_, n_, \[Mu]_, \[Nu]_, s_] := 
 Exp[4 n s q] q^(4 n \[Nu] - 
     1) (1 - q)^(4 n \[Mu] - 1)/(Gamma[4 n \[Mu]] Gamma[
      4 n \[Nu]] Hypergeometric1F1Regularized[4 n \[Nu], 4 n (\[Mu] + \[Nu]), 
      4 n s])

Check normalization: ok for all parameters if {n>0,[Mu]>0,[Nu]>0}

Integrate[p[q, n, \[Mu], \[Nu], s], {q, 0, 1}, 
 Assumptions -> {n > 0, \[Mu] > 0, \[Nu] > 0}]

(* Out[3]= 1 *)

Calculate the k-th moment (k>=0)

Integrate[q^k p[q, n, \[Mu], \[Nu], s], {q, 0, 1}, 
 Assumptions -> {n > 0, \[Mu] > 0, \[Nu] > 0, k >= 0}]

(* Out[6]= (Gamma[k + 4 n \[Nu]] Hypergeometric1F1Regularized[k + 4 n \[Nu], 
  k + 4 n (\[Mu] + \[Nu]), 4 n s])/(
Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[4 n \[Nu], 4 n (\[Mu] + \[Nu]), 
  4 n s]) *)

Now look at the normalization with NIntegrate for very simple values of the parameters

Table[NIntegrate[p[q, n, 1, 1, 1], {q, 0, 1}], {n, 1, 20}]

During evaluation of In[63]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in q near {q} = {0.673813}. NIntegrate obtained 0.8149444193581905and 0.0004882716173466382 for the integral and error estimates. >>

(* 
Out[63]= 
{1., 1., 1., 1., 1., 1., 1., 1., 1.00009, 1.00112, 1.02324, 1.25255, 2.32575, 
-0.0207066, 0.00601814, 0.814944, 1., 1., 1., 1.} 
*)

Hence I recommend to look very carefully at the parameter values of your problem, and check first if NIntegrate gives 1.

PS: why do you attempt to use NIntegrate at all?

Regards, Wolfgang

Answer 2

Not sure what is going on with the different results of Integrate and NIntegrate. This does not mean that the analytic form of $C$ is erroneous.

Note that plotting the likelihood function (using the expression of $C$ provided by Integrate and the parameter values you used) over a resticted range of $s$ (instead of $[0,1]$) clearly shows that the log-likelihood is maximized for a non-zero value.

Plot[Log[(E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]))/(
   Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
     4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s])] //. {q -> 0.75, 
   n -> 100/(40 \[Nu]), \[Nu] -> \[Mu], \[Mu] -> 10^(-6)}, {s, 
  10^(-7), 10^(-5)}, AxesLabel -> {"s", "log(likelihood)"}]

The exact value can be found by FindMaximum (a good initial estimate is required)

FindMaximum[
 Log[(E^(4 n q s) (1 - q)^(-1 + 4 n \[Mu]) q^(-1 + 4 n \[Nu]))/(
   Gamma[4 n \[Mu]] Gamma[4 n \[Nu]] Hypergeometric1F1Regularized[
     4 n \[Nu], 4 n (\[Mu] + \[Nu]), 4 n s])] //. {q -> 0.75, 
   n -> 100/(40 \[Nu]), \[Nu] -> \[Mu], \[Mu] -> 10^(-6)}, {s, 
  10^(-6)}]

{1.66, {s -> 2.75*10^-6}}