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

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.

-
I suspect Integrate is coming up with a bad antiderivative. Will investigate. –  Daniel Lichtblau Jan 10 at 21:07
@Daniel: I don't think so. Integrate behaves well (Version 8) (see my answer). –  Dr. Wolfgang Hintze Jan 10 at 22:29
@Dr.WolfgangHintze Thanks for the observation and response, that might save me some trouble. (That, or double it, if both Integrater AND NIntegrate are derelict. We'll see.) –  Daniel Lichtblau Jan 11 at 20:01

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

-
Thank you. I used NIntegrate when I realized that the individual probabilities I could calculate would never sum up to 1 although Integrate said the opposite. I didn't know the concept of numerical stability. –  Remi.b Jan 10 at 23:16
Your answer already helps a lot but there's a point I still don't get. It seems to me that the concept of numerical stability only applies to the function NIntegrate and not to the function Integrate neither when I just calculate the probability $p(q \space |\space μ,ν,n,s)$. So I still don't fully understand why my plot shows a curve that seem to integrate to a value that is closer to the result of NIntegrate than to 1. Does my misunderstanding make sense? –  Remi.b Jan 11 at 16:27

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}}

-
Thanks a lot. That will help as well. I was kinda lost and I realized my post may contain two distinct questions. You answered the ML-method related question. +1 –  Remi.b Jan 10 at 23:22