# 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. Jun 22, 2020 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 Jun 23, 2020 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 Jun 23, 2020 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 Jun 23, 2020 at 1:11
• I'll expand the question, spelling out the nested integrals, to make my point clear Jun 23, 2020 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 Jun 23, 2020 at 2:17