as the distribution
<span class="math display">\[
Z \sim \mathsf{multivariate\ normal}\left(
\begin{bmatrix}
0 \\
0
\end{bmatrix}, \
\begin{bmatrix}
1 &amp; \rho
\\
\rho &amp; 1
\end{bmatrix}
]\right).
\]</span></p>
<p>The following Stan program implements this function,</p>
<pre><code>real binormal_cdf(real z1, real z2, real rho) {
  if (z1 != 0 || z2 != 0) {
    real denom = fabs(rho) &lt; 1.0 ? sqrt((1 + rho) * (1 - rho)) : not_a_number();
    real a1 = (z2 / z1 - rho) / denom;
    real a2 = (z1 / z2 - r