Computing the pdf of two independent Gaussians conditioned on being inside the circle Suppose that $Z_1,Z_2$ are iid $N(0,1)$. For a fixed $r \ge 0$, I am interested in computing the conditional probability distribution: $P(Z_1\le x, Z_2 \le y \,| \, Z_1^2 + Z_2^2 \le r^2)$ for any $x,y \in \mathbb{R}^2$.
This entails computing the probability $P(Z_1\le x,\, Z_2 \le y, \,  Z_1^2 + Z_2^2 \le r^2)$. I am interested in computing the pdf of the conditional distribution above and so it is necessary that I consider different values of $x,y$. Upon doing the calculations, it turns out that the region of integration becomes messy as it greatly varies depending on $x,y$.
Are there any tricks that I can use to evaluate the above probability in a much simpler way? I am looking for a trick that easily extends to the case when I have $Z_1,Z_2,...,Z_n \sim N(0,1)$ iid with the same condition that $Z_1^2 + \dots + Z_n^2 \le r^2$. In addition, I am also interested in the case when the sphere has its center away from the origin.
Any suggestions?
 A: STANDARD NUMERICAL APPROACH
It is kinda bizarre in fact, but if I understand the situation you're trying to evaluate the integral
$$\iint_D f(u)f(v) \, dA$$
with $f(t)$ the standard normal density ($f(t)=\tfrac1{\sqrt{2\pi}}e^{-t^2/2}$), and while is not an issue (I assume) that integration has to be numerically done, it is an issue that $D$ is not nice to describe. Anyway, I believe that $\partial D$ is the union of segments and arcs of the circunference in such a way that the worst case is two of each type, and that the border is not that hard to parameterize since the arcs are the images of intervals —of the form $\left(\arcsin (y/r),\arccos (x/r)\right)$— and variations of this depending on each case by the function $\sigma(t)=(r \cos t,r \sin t)$.
So if parameterization of the border for each case is messy but not that much, you could consider using Green's Theorem and do
$$\iint_D f(u)f(v) \, dA=\frac12\int_{\partial D} F(x)f(y)\, dy - f(x)F(y)dx,$$
where $F(t)$ is the standard normal distribution function (a primitive of $f$), and since the last integral ends up being a Riemann integral, numerical integration is now quite straightforward.
Of course there are still a couple of different cases to study, but considering some symmetries and other things, is not that bad as before. For instance, all the cases in which $x_0^2+y_0^2\leq r^2$, no mater the quadrant on which $(x_0,y_0)$ lies correspond to integration over a region bounded by an arc of circumference closed by a horizontal segment above (at ordinate $y_0$) and a vertical segment to the right (at abscissa $x_0$). These are parameterized by:
$$\sigma_1(t)=(x_0,t),\quad y_1 \leq t\leq y_0,$$
$$\sigma_2(t)=(x_0-t,y_0),\quad 0\leq t\leq x_0-x_1,$$
$$\sigma_3(t)=(r\cos t,r\sin t),\quad t_y \leq t \leq t_x,$$
where
$$x_1=-\sqrt{r^2-y_0^2}$$
$$y_1=-\sqrt{r^2-x_0^2}$$
$$t_x=2\pi-\arccos(x_0/r)$$
$$t_y=\pi-\arcsin(y_0/r).$$
The only other cases worth analyzing occur when $x_0^2+y_0^2>r$ and are:


*

*$x_0>r$ and $y_0\leq r$

*$y_0>r$ and $x_0\leq r$ (by symmetry it can be approached with the previous case)

*$0<y_0\leq r$ and $0<x_0\leq r$ (here $\partial D$ has two arcs and two segments).


The remaining cases of course give $0$ or $1$ conditional probability.

MONTE CARLO APPROACH
If you prefer to use a Monte Carlo method, just generate two standard normal values $x_1$ and $y_1$ independently and check whether $x_1\leq x$, $y_1\leq y$ and $x_1^2+y_1^2\leq r^2$ are ALL true or not. Then generate $x_2$ and $y_2$ independently and standard normal and check again... and so on. The number proportion of 'successes' to repetitions will be an estimate of that probability. Remember to divide by $P(Z_1^2+Z_2^2\leq r^2)$ to calculate the conditional probability. That one can be analytically obtained or you can also calculate it's rate of success to repetitions in the previous simulation as an estimate, and then divide both estimations, or divide directly both numbers of 'successes'.
Do NOT try to generate random points of the disk at once, since most likely you will not respect the distribution of both coordinates and/or their independence. The safest (and perhaps essentially the only) way to simulate this model is to pick the two coordinates independently and eventually disregard those that fall not in the circle.

I haven't implemented anyone yet, but the Monte Carlo method is way easier to program. Anyway, it might be more time consuming than the other ---once programmed--- for a given precision (which would be an actual precision in the first case and just a good confidence precision in the second case).
If it were my choice to take, I would use the Monte Carlo method in the beginning of the investigation, but for a more profound and rigorous work I'd work on implementing the first option.
