# 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. Jan 10, 2015 at 21:07
• @Daniel: I don't think so. Integrate behaves well (Version 8) (see my answer). Jan 10, 2015 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.) Jan 11, 2015 at 20:01
• I am mystified by the Hypergeometric1F1Regularized function. I would like to put this in an answer so you can see my text better, but I really have a question so think it would be wrong to post it as an answer. It is too long for one comment so I will break it into parts. Jul 26, 2015 at 17:15
• When I Integrate[p[q,C,n,s,[Nu],[Mu]],{q,0,1},Assumptions->{C,n,s,[Nu],[Mu]}[Element]Reals] it results in ConditionalExpression[C Gamma[4n[Mu]]Gamma[4n[Nu]]Hypergeometric1F1Regularized[4n[Nu],4n([Mu]+[Nu]),4ns],(n>0&&[Mu]>0&&[Nu]>0)||(n<0&&[Mu]<0&&[Nu]<0)]. Jul 26, 2015 at 17:15

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= 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= (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:= 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=
{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. Jan 10, 2015 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? Jan 11, 2015 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 Jan 10, 2015 at 23:22