# Numerical integration with NIntegrate

pj[mu_, sigma_] :=
Table[
NIntegrate[
(1/(1 + Exp[-x]))^j
*
(1 - 1/(1 + Exp[-x]))^(10 - j)
/sigma/Sqrt[2*Pi]
*
Exp[-(x-mu)^2/2/sigma^2],
{x, -Infinity, Infinity}
],
{j, 0, 10}
]

pj[-1.9575, 0.3432]
pj[-4, 9]


I need this expression for a global optimization. But for some values, NIntegrate fails. I have tried a couple of options, such as MaxPoints and MaxRecursion, with no luck.

During the optimization, I still get all kinds of warnings and, sometimes, it just failed.

Below is a copy of the full code, adjusted to helpful comments and suggestion. But I am finding it hard to modify the pj term to be able to evaluate quickly (hopefully), without any warnings.

ClearAll["Global*"]

link1 = {84, 54, 36, 21};

skinks = {56, 19, 28, 18, 24, 14, 9};
kskinks = 7;

taxicabsA:={142, 81, 49, 7, 3, 1};
ktaxicabsA:= 10;

B1999:={11,12,10,4,4,1,4,2,3,3,0,2,4,1,1,0,1,1,1,2,0,1,0,0,0,0,0,0,0,1,3};
kB1999:= 50;

(* Modified Version *)
mylikLNB2[data_, kdata_, f0_, mu_, sigma_] := Module[
{K, pj, fj, j, N0, loglik, above, below},

K = kdata;

pj[muin_, sigmain_] :=

Table[
With[
{integrand=Simplify[
(1/(1 + Exp[-x]))^j
*
(1 - 1/(1 + Exp[-x]))^(K - j)/sigmain/Sqrt[2*Pi]
*
Exp[-(x - muin)^2/2/sigmain^2]
]},
NIntegrate[integrand, {x, -Infinity, Infinity}]
]
,{j, 0, K}
];

pj = pj[mu,sigma]*Table[Binomial[K, j],{j, 0, K}];

N0 = Sum[fj[[j]], {j, 1, Length[fj]}] ;
above = LogGamma[N0 + 1] ;
below = Sum[LogGamma[fj[[j]] + 1], {j, 1, Length[fj]}] ;
loglik = fj.Log[pj] ;

loglik = above - below + loglik

]

(* Could try other datasets *)
testfun[f0_, mu_, sigma_] :=
mylikLNB2[taxicabsA, ktaxicabsA, f0, mu, sigma];

(* Desired Result *)
testfun[106.1337,-1.9575,0.3432]

ans=NMaximize[
{testfun[f0, mu, sigma],
f0 >= 0 && sigma>0 },
{f0, mu, sigma},
Method -> {"RandomSearch", "SearchPoints" -> 50}
]

(* Gets the Desired Result, but not quite 'fast' enough *)


Anyone has a solution?

• What sort of values does it fail for? What warnings do you get? Sep 30 '14 at 19:03
• Say with pj[-4,9]. I get "SystemException["MemoryAllocationFailure" Sep 30 '14 at 19:11
• Could set Method->"NewtonCotesRule". The default behavior is not pretty. I'll file a bug report on it. (I'm not certain it's a bug but seems like it's worth a report at least.) Sep 30 '14 at 19:32
• @DanielLichtblau No success for Method->"NewtonCotesRule" on MMA10 . But thanks! Sep 30 '14 at 19:46
• It worked for me. But I had some formatting issues when I cut-pasted your input and I may have guessed incorrectly as to which sigmas went where. I did get the memory exception with default handling, and I did get numerical results with that explicit method setting. Best I can suggest is maybe try other methods and see if things improve. Sep 30 '14 at 20:11

In V10, the kernel seems unhappy with the exact quantities being given to NIntegrate. A work-around is to define pj with inexact coefficients.

pj[j_Integer, mu_?NumericQ, sigma_?NumericQ] :=
NIntegrate[
(1./(1. + Exp[-x]))^j*
(1. - 1./(1. + Exp[-x]))^(10 - j)/sigma/Sqrt[2.*Pi]*
Exp[-(x - mu)^2/2./sigma^2],
{x, -Infinity, Infinity}]

pj[#, -4., 9.] & /@ Range[0, 10]

{0.551423, 0.00482382, 0.000591354, 0.000165584, 0.0000811702,0.0000636039,
0.0000776634, 0.000150861, 0.00050613, 0.00371322, 0.226463}


This approach has the advantage of being very fast.

• I dont exactly understand the term "inexact". But it seems to be working. I will try it with my complete code. How do you find it "being very fast"? Oct 1 '14 at 21:39
• @m_goldberg: I have modified a full copy of the problem I am trying to achieve, if you dont mind taking a further look. It works too, but if I want to use more starting points for the optimization, it is not quite 'fast'. Thanks. Oct 1 '14 at 22:21
• @ChenStatsYu. In Mathematica, integers and rationals like 1 and 1/10 are exact, and decimals such as 1. and 0.1 are by default machine floating point numbers. The latter will (usually) produce faster but less accurate calculations. It is a oversimplification, but it can be useful to think of NIntegrate as being designed to work with inexact numbers and Integrate as being designed to work with exact ones. Oct 2 '14 at 8:19
• @ChenStatsYu. I don't see any modified full copy to look at. Oct 2 '14 at 8:21
• @m_goldberg Sorry that I didnt point out where the full code is. I posted it as an answer, but then modified it, with the whole problem I am trying to archive. I added a comment and @ you, but didnt think it works. Thanks for the explanation on 'exact'. Oct 2 '14 at 11:09
ClearAll[pj];
pj[mu_?NumericQ, sigma_?NumericQ] :=
Table[NIntegrate[(1/(1 + Exp[-x]))^ j*(1 - 1/(1 + Exp[-x]))^(10 - j)/sigma/Sqrt[2*Pi]*
Exp[-(x - mu)^2/2/sigma^2], {x, -Infinity, Infinity}], {j, 0., 10.}]
pj[-1.9575, 0.3432]
pj[-4, 9]

(* {0.273692, 0.0356199, 0.0051466, 0.000824808, 0.000146484, \
0.0000288024, 6.26425*10^-6, 1.50561*10^-6, 3.99548*10^-7,
1.16965*10^-7, 3.77401*10^-8}

{0.551423, 0.00482382, 0.000591354, 0.000165584, 0.0000811702, \
0.0000636039, 0.0000776634, 0.000150861, 0.00050613, 0.00371322, \
0.226463}*)

• I have exactly the same result with you. But pj[-4, 9] still does not work on my MMA 10. Sep 30 '14 at 19:24
• @ChenStatsYu The above code works here. But not yours. Sep 30 '14 at 19:25
• Which MMA do you have? I am using MMA 10 on Win7 X64. Sep 30 '14 at 19:26
• @ChenStatsYu v9. But please check that you're using both the ?NumericQ tests and the modified {j, 0., 10.} range Sep 30 '14 at 19:29
• I have added a picture to my OP at the top, I hope you can see it. It shows the result from my MMA10 on win7 X64. Sep 30 '14 at 19:29

If you examine what is going on with the option EvaluationMonitor :> Print[x], you discover that the last value of x is -6.337380909737406*^37, which causes an overflow. It's not the first x for which an overflow occurs, so overflow itself cannot be the reason. However the expression for which there is overflow is

(1 - 1/(1 + E^(-x)))^10


(or more specifically E^(-x)) which can be simplified to

1/(1 + E^x)^10


For this expression we get underflow instead of overflow. But in integration at machine precision, this may be treated as zero (or, more precisely, a negligible error). This formula happens to be achieved with Simplify.

Thus, with the following definition, I get pj[-4, 9] to evaluate:

pj[mu_, sigma_] := Table[
With[{integrand = (1/(1 + Exp[-x]))^j * (1 - 1/(1 + Exp[-x]))^(10 - j)/sigma/Sqrt[2*Pi]*
Exp[-(x - mu)^2/2/sigma^2] // Simplify},
NIntegrate[integrand, {x, -Infinity, Infinity}]
],
{j, 0, 10}]

pj[-4, 9]
(*
{0.551423, 0.00482382, 0.000591354, 0.000165584, 0.0000811702,
0.0000636039, 0.0000776634, 0.000150861, 0.00050613, 0.00371322, 0.226463}
*)

ClearAll["Global*"]

link1 = {84, 54, 36, 21};

skinks = {56, 19, 28, 18, 24, 14, 9};
kskinks = 7;

taxicabsA:={142, 81, 49, 7, 3, 1};
ktaxicabsA:= 10;

B1999:={11,12,10,4,4,1,4,2,3,3,0,2,4,1,1,0,1,1,1,2,0,1,0,0,0,0,0,0,0,1,3};
kB1999:= 50;

(* Modified Version *)
mylikLNB2[data_, kdata_, f0_, mu_, sigma_] := Module[
{K, pj, fj, j, N0, loglik, above, below},

K = kdata;

pj[j_Integer, mu, sigma] :=
NIntegrate[
(1./(1. + Exp[-x]))^j
*
(1. - 1./(1. + Exp[-x]))^(K - j)/sigma/Sqrt[2.*Pi]
*
Exp[-(x - mu)^2/2./sigma^2],
{x, -Infinity, Infinity}
];

pj = pj[#, mu, sigma] & /@ Range[0, K];

pj = pj*Table[Binomial[K, j],{j,0,K}];

N0 = Sum[fj[[j]], {j, 1, Length[fj]}] ;
above = LogGamma[N0 + 1] ;
below = Sum[LogGamma[fj[[j]] + 1], {j, 1, Length[fj]}] ;
loglik = fj.Log[pj] ;

loglik = above - below + loglik

]

(* Could try other datasets *)
testfun[f0_, mu_, sigma_] :=
mylikLNB2[taxicabsA, ktaxicabsA, f0, mu, sigma];

(* Desired Result *)
testfun[106.1337,-1.9575,0.3432]

ans=NMaximize[
{testfun[f0, mu, sigma],
f0 >= 0 && sigma>0 },
{f0, mu, sigma},
Method -> {"RandomSearch", "SearchPoints" -> 50}
]

(* Gets the Desired Result, but not quite 'fast' enough *)


This is a modified version based on @m_goldberg suggestion. I have presented a full copy of the code. I am trying to fit the dataset(s) with a distribution.

The NMaximize optimizer gets the answer with just about 20 multistart, sometimes have to have more starts, like 50? 100? 200? (or more). But it is still not quite 'fast' enough.With 5+ "SearchPoints", it takes significantly more time. And it still presents warnings during the optimization.

Below is a version that does not gives any warning, but really really slow to run.

pjNEW[j_Integer, mu_?NumericQ, sigma_?NumericQ, K_] :=
NIntegrate[(1./(1. + Exp[-x]))^
j*(1. - 1./(1. + Exp[-x]))^(K - j)/sigma/Sqrt[2.*Pi]*
Exp[-(x - mu)^2/2./sigma^2], {x, -Infinity, Infinity},
MinRecursion -> 9];

(*Modified Version*)
mylikLNB2[data_, kdata_Integer, f0_?NumericQ, mu_?NumericQ,
sigma_?NumericQ] :=
Module[{K, pj, fj, j, N0, loglik, above, below}, K = kdata;
pj = pjNEW[#, mu, sigma, K] & /@ Range[0, K];
pj = pj*Table[Binomial[K, j], {j, 0, K}];
N0 = Sum[fj[[j]], {j, 1, Length[fj]}];
above = LogGamma[N0 + 1];
below = Sum[LogGamma[fj[[j]] + 1], {j, 1, Length[fj]}];
loglik = fj.Log[pj];
loglik = above - below + loglik]

(*Could try other datasets*)
testfun[f0_, mu_, sigma_] :=
mylikLNB2[taxicabsA, ktaxicabsA, f0, mu, sigma];

(*Desired Result*)
testfun[106.1337, -1.9575, 0.3432]

ans = NMaximize[{testfun[f0, mu, sigma], f0 >= 0 && sigma > 0}, {f0,
mu, sigma}]

(*Gets the Desired Result,but not quite'fast' enough*)
`
• @m_goldberg Sorry if I didnt point out where the full code is. Thanks for the explanation on the 'exact'. Oct 2 '14 at 11:07