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 enter image description here

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

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

1 Answer 1


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$)

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.