Derivative of bivariate normal CDF with common mean parameters I am trying to calculate a derivative of the form $\frac{d}{dz}\Phi_2(\mu_1(z),\mu_2(z),\rho)$, where $\Phi_2$ is the standard bivariate normal CDF. 
I am thinking it might be an application of a double Leibniz rule (example here) but the lower limit of integration in the CDF is not finite.
Edit: To be clear, I define the standard bivariate normal CDF as
$$\Phi_2(\mu_1(z),\mu_2(z),\rho) = \int_{-\infty}^{\mu_1(z)}\int_{-\infty}^{\mu_2(z)}\frac{1}{2 \pi \sqrt{1-\rho^2}}\exp[-\frac{1}{2(1-\rho^2)}(x^2+y^2-2\rho x y)]\,dy\,dx$$ where $\mu_1(z)$ and $\mu_2(z)$ are functions differentiable with respect to z.
 A: We can obtain a nice closed-form answer simply by applying definitions and the most basic result of linear regression theory: no calculation is needed.
First consider more generally what happens to $\Phi_2(x,y,\rho)$ when $x$ is changed to $x+\mathrm{d}x$ for a positive infinitesimal $\mathrm{d}x.$  By definition, $\Phi_2(x,y,\rho)$ is the total probability in the square where $X\le x$ and $Y\le y.$  The difference therefore is the total probability in the infinitesimal half-infinite vertical strip $S(x,\mathrm{d}x,y)$ delimited at the left by $X=x,$ the right by $X=x+\mathrm{d}x,$ and the top by $y.$ 

Colors denote values of the bivariate density with $\rho=2/5.$  Its regression line $y=\rho x$ is shown in white. The colored vertical strip on the right is $S(x,\mathrm{d}x,y).$  The colored horizontal strip at the top is the analog of this strip after exchanging $X$ and $Y.$
The theory of linear regression teaches us that in the linear regression line for $\Phi_2,$ the conditional distribution of $Y\mid X=x$ has mean $\rho x$ and variance $1-\rho^2.$  Since the probabilities in this vertical strip do not appreciably change across it from left to right, and must be proportional to the standard Normal density $\phi(x),$ then in terms of the standard normal CDF $\Phi$ they must therefore equal
$$\eqalign{
\Pr(S(x,\mathrm{d}x,y)) &= \Pr(Y\le y\mid X\in[x,x+\mathrm{d}x))\Pr(X\in[x,x+\mathrm{d}x)) \\
&= \Phi\left(\frac{y-\rho x}{\sqrt{1-\rho^2}}\right)\,\phi(x)\mathrm{d}x .
} \tag{*}$$
It's clear that the same result holds for negative infinitesimal $\mathrm{d}x.$
Reversing the roles of $X$ and $Y$ changes nothing in the reasoning and only swaps $x$ and $y$ in the result: from the symmetric expression for the bivariate Normal density, the situation is identical. 
Specializing to the question, let $\mathrm{d}z$ be an infinitesimal change in $z.$  By definition of the derivative, this induces simultaneous infinitesimal changes in $x$ and $y$ given by
$$\eqalign{
\mathrm{d}x &= \mathrm{d}\mu_1(z) = \mu_1^\prime(z)\mathrm{d}z; \\
\mathrm{d}y &= \mathrm{d}\mu_2(z) = \mu_2^\prime(z)\mathrm{d}z.
}$$
Together this creates two infinitesimal strips between the rectangles $X\le \mu_1(z), Y\le \mu_2(z)$ and $X\le \mu_1(z+\mathrm{d}z), Y \le \mu_2(z + \mathrm{d}z),$ as shown in the figure. Their total area is given by two applications of $(*)$ upon substituting $(x,y)=(\mu_1(z),\mu_2(z))$ and $(\mathrm{d}x,\mathrm{d}y) = (\mu_1^\prime, \mu_2^\prime)\mathrm{d}z:$
$$\eqalign{
& \mathrm{d}\Phi_2(\mu_1(z),\mu_2(z),\rho) \\ &=\Phi\left(\frac{\mu_2(z)-\rho \mu_1(z)}{\sqrt{1-\rho^2}}\right)\phi(\mu_1(z))\mu_1^\prime(z)\mathrm{d}z + \Phi\left(\frac{\mu_1(z)-\rho \mu_2(z)}{\sqrt{1-\rho^2}}\right)\phi(\mu_2(z))\mu_2^\prime(z)\mathrm{d}z.
}$$
Dividing both sides by $\mathrm{d}z$ gives the desired derivative.
A: This post here seems relevant. For the case $\mu_1(z)=\mu_2(z)=z$, what you are asking for is essentially the density of $\max(X,Y)$ as is shown here. 
Let $(X,Y)\sim\mathcal{B\,N}(\mu_x=0,\mu_y=0,\sigma^2_x=1,\sigma^2_y=1,\rho)$ having joint pdf $f_{X,Y}$. 
Using $\Phi$ and $\phi$ to denote CDF and PDF respectively of a standard normal variate as usual.
Then, 
\begin{align}
\Phi_2(\mu_1(z),\mu_2(z))&=\int_{-\infty}^{\mu_2(z)}\int_{-\infty}^{\mu_1(z)}f_{X,Y}(x,y)\,\mathrm{d}x\,\mathrm{d}y
\\&=\Pr(X\leqslant\mu_1(z),Y\leqslant\mu_2(z))
\\&=\int_{-\infty}^{\mu_1(z)}\Pr(Y\leqslant\mu_2(z)\mid X=x)\,\phi(x)\,\mathrm{d}x
\\&=\int_{-\infty}^{\mu_1(z)}\Phi\left(\frac{\mu_2(z)-\rho x}{\sqrt{1-\rho^2}}\right)\,\phi(x)\,\mathrm{d}x
\end{align}
The last equality follows from the fact that $[Y\mid X=x]\sim\mathcal N(\rho x,1-\rho^2)$.
Using Leibniz rule, $\frac{\mathrm{d}}{\mathrm{d}z}\Phi_2(\mu_1(z),\mu_2(z))$ equals $$\int_{-\infty}^{\mu_1(z)}\frac{\partial}{\partial z}\Phi\left(\frac{\mu_2(z)-\rho x}{\sqrt{1-\rho^2}}\right)\phi(x)\,dx+\mu_1'(z)\Phi\left(\frac{\mu_2(z)-\rho \mu_1(z)}{\sqrt{1-\rho^2}}\right)\phi(\mu_1(z))$$
The 'contribution' due to the lower limit of the integral is zero as the lower limit is not a function of $z$. I think the final answer simplifies further in terms of $\Phi(\cdot)$ and $\phi(\cdot)$. 
