# Sampling from conditional multivariate normal distribution ("conditioned on ellipses" / stratified sampling)

I have an issue related to stratified sampling. To explain, I need to start with

one-dimensional case:

Let $$N\sim \mathcal{N}(0,1)$$, i.e., a standard normal variable with cdf $$F$$. Partition $$\mathbb{R}$$ into $$m$$ disjoint intervals $$A^1=(-\infty,a_1], A^2=(a_1,a_2], \ldots, A^m=(a_{m-1},+\infty)$$ so that $$P(N\in(a_{j-1},a_j])=1/m$$. This can be done by setting $$a_k = F^{-1}(k/m), k=1,\ldots,m-1$$ (here I assume I have $$F^{-1}$$, havig this in computer packages is enough)

Task: How to simulate $$N | N\in A^k$$, i.e., $$N$$ conditioned that it is in the interval $$(a_{k-1},a_k$$)?

Solution: Let $$V=a_{k-1}+(a_k-a_{k-1})U,$$ where $$U\sim\mathcal{U}(0,1)$$. Then we have what we wanted: $$F^{-1}(V)=_d (N|N\in A^k)$$

multivariate-dimensional case: Let $$\bf{N}\sim\boldsymbol{\mathcal{N}}({\boldsymbol\mu},\Sigma)$$ be a $$d\geq 2$$ dimensional normal distribution. Let us stick to $$d=2$$. Partition $$\mathbb{R}^2$$ into $$m$$ disjoint sets $$A^1, \ldots, A^m$$ in a specific way: Namely, $$A^1$$ is the area (enclosed by an be an ellipse) such that $$P({\bf N}\in A^1)=1/m$$. Then $$A^2$$ is the area between two ellipses, such that $$P({\bf N}\in A^2)=1/m$$ and so on (last $$A^m$$ is "everything else"). I hope description is clear.

Task: How to simulate $${\bf N} | {\bf N}\in A^k$$ ??

It is clear how to proceed in case when $$\boldsymbol{\Sigma}=\boldsymbol{I}$$ (identity matrix). Then, still in case $$d=2$$, the density can be presented in polar coordinates $$f(r,\theta) = g(r)h(\theta)=r e^{-{1\over 2}r^2} {1\over 2\pi},\quad r\geq 0, \theta\in(0,2\pi)$$ One can calculate that (denote $$N=(N_1,N_2)$$) $$P( N_1^2+N_2^2< t)=1-e^{-{1\over 2} t^2}$$ and thus $$P\left( \sqrt{-2\ln\left(1-{k\over m}\right)} \leq N_1^2+N_2^2< \sqrt{-2\ln\left(1-{k+1\over m}\right)}\right)={1\over m},$$ thus one can sample $$\theta\sim\mathcal{U}(0,2\pi)$$ and $$r$$ correspondingly. Below m=5, thus 5 stratas are presented, 100 points simulated in stratas 1,2,3,4

However, I do not know how to proceed with general $$\boldsymbol{\mu}, \boldsymbol{\Sigma}$$. Any help appreciated.

• What about simulating for $\Sigma=\mathbf I_p$ and $\mu=0_p$, and then applying a location scale transform? Jan 9, 2020 at 22:18

Yes, actually, it was quite simple. Cholesky decomposition: $$\boldsymbol{\Sigma}=\boldsymbol{A}^T\boldsymbol{A}$$, then each point $$(x,y)$$ simulated from $$\boldsymbol{N}(\boldsymbol{0},\boldsymbol{I})$$ conditioned on strata $$A^i$$ (as described in question) is transformed to $$(x,y)\boldsymbol{A}+\boldsymbol{\mu}$$. Below example with $$\boldsymbol{\mu}=(0,0)$$ and $$\Sigma=\left(\begin{array}{ll}0.5 & 0.5 \\ 0.5 & 1\end{array}\right)$$.
Left image = points $$\boldsymbol{Z}\sim\boldsymbol{N}(\boldsymbol{0},\boldsymbol{I}), \qquad$$ right image = $$\boldsymbol{Z}\boldsymbol{A}$$
(100 points in each strata $$A^k, k=1,\ldots,5$$)