# Using NIntegrate to reproduce NProbability over joint Gaussian distribution

Consider a random vector {s,c} with a bivariate normal distribution. For a vector of positive scalars {a, ß, σz}, I'm interested in calculating (numerically) the probability

NProbability[c < (1 - CDF[NormalDistribution[ ß*(s - c), σz], a]),{s,c} \[Distributed] BinormalDistribution[{μs, μc}, {σs, σc}, ρ]]


Is there a way to write this same calculation using only NIntegrate?

What I've done so far

I've tried re-writing the probability, solving for s on one side of the inequality, and nesting the integrals:

f1[c_?NumericQ,μs_, μc_, σs_, σc_, σz_, ρ_, a_, ß_]:=NIntegrate[PDF[BinormalDistribution[{μs, μc}, {σs, σc}, ρ],{s,c}],{s, c + (a-σz*InverseCDF[NormalDistribution[],1-c])(ß)^-1,\[Infinity]},]

f2[μs_, μc_, σs_, σc_, σz_, ρ_, a_, ß_]:=NIntegrate[f1[c,μs, μc, σs, σc, σz, ρ, a, ß],{c,-\[Infinity],\[Infinity]}]


This approach unfortunately doesn't work because the computation gets stuck with InverseCDF[NormalDistribution[],1-c] for c below zero or above one.

Parameter values

The scalars and distribution parameters are not important. Here is a starting set of values that can be used for reference:

{μs, μc, σs, σc, σz, ρ, a, ß} = {.35, .5, 1.1, 1.2, 1.3, .25, 1, .5}

• Do you have values for σs, σc, ρ, ß, a, σz, μs, μc... etc? NIntegrate expects these to be numerical - as does NProbability. – flinty Jun 22 '20 at 23:48
• @flinty Yes, I do have values for these. I get NProbability to give me the required calculation - but I’d like to reproduce that calculation with NIntegrate – OO_SE Jun 23 '20 at 0:11
• Can you give the values then? The code in your NIntegrate doesn't work or has been formatted wrongly. I get ... not a valid limit of integration – flinty Jun 23 '20 at 0:42
• Thanks, I've edited the question and added an example set of parameter values, in case that's helpful. The code for NIntegrate, as written above, doesn't work - as you say, it rejects the limit of integration. But if you write this up nesting two integrals (one over s and a second one over c) then you see that if fails due to c<0 and c>1 – OO_SE Jun 23 '20 at 1:11
• I'll expand the question, spelling out the nested integrals, to make my point clear – OO_SE Jun 23 '20 at 1:14

{μs, μc, σs, σc, σz, ρ, a, ß} =
{.35, .5, 1.1, 1.2, 1.3, .25, 1, .5} // Rationalize;

ineq = c < (1 - CDF[NormalDistribution[ß*(s - c), σz], a]) //
FullSimplify

(* 2 c < Erfc[(5 (2 + c - s))/(13 Sqrt)] *)


NIntegrate is an available Method for NProbability

NProbability[ineq,
{s, c} \[Distributed]
BinormalDistribution[{μs, μc}, {σs, σc}, ρ],
WorkingPrecision -> 30,
Method -> {"NIntegrate",
{MinRecursion -> 15, MaxRecursion -> 25}}] (* 0.419500831140737615758538073412 *)


However, the "MonteCarlo" method does not result in warning messages.

NProbability[ineq,
{s, c} \[Distributed]
BinormalDistribution[{μs, μc}, {σs, σc}, ρ],
WorkingPrecision -> 30,
Method -> {"MonteCarlo", "SamplingIncrement" -> 10^4}]

(* 0.415478873239436619718309859154930 *(

• Thanks for this. I do get NProbability to give me the result, as you point out her, but what I'd like is to be able to reproduce this writing down the calculation explicitly in terms of NIntegrate, rather than the NProbability function – OO_SE Jun 23 '20 at 2:17