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}

{μ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[2])] *)

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

(* 0.419500831140737615758538073412 *)

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

(* 0.415478873239436619718309859154930 *(