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.
