# Moments of Rejection Sampling

I had tried out the following method for Gaussian sample generation using Uniform samples. I don't understand how the variance scales in my simulation results. The setup is described below.

Assume that $$X_1,X_2,\dots,X_n \sim Unif(-1,1)$$ iid samples are given with mean $$\mu_X$$ . For any $$m>1$$ fixed, define an increasing sequence $$b_1 such that $$b_i \in (-1,1)$$ for every $$i$$. Now consider the partition of $$\mathbb{R}$$ formed out of these $$b_i$$'s $$\{(-\infty,b_1),(b_1,b_2),\dots,(b_m, \infty)) \}$$. Denote the $$m+1$$ intervals as $$J_1,J_2,\dots J_{m+1}$$. Define for each $$i=1,2,\dots,m+1$$, $$\pi_i=\displaystyle\int_{J_i} \mathcal{N}(0,1) dx$$ where $$\mathcal{N}(0,1)$$ is the standard normal distribution.

The idea now is for each $$i$$, if $$X_i \in J_j$$, define $$Y_i=\begin{cases} X_i, \text{ w.p } \text{ } \pi_j \\ \mu_X, \text{ w.p } \text{ } 1-\pi_j \end{cases}$$ where w.p means with probability.

In other words the samples are accepted with probability $$\pi_j$$ if $$X_i$$ is in the $$j^{\text{th}}$$ partition. I am interested in $$\mathbb{E}(Y_i)$$ and the $$\mathbb{E}(Y_i-\mathbb{E}(Y_i))^2$$

• I am confused by what you mean by "accepted". Say $X_1 \in J_5, X_2 \in J_{17}, X_3 \in J_0, X_4 \in J_6$ then what are $Y_1, Y_2, Y_3, Y_4$ respectively? In particular I was guessing "w.p. $\pi_j$" means "with probability $\pi_j$" but then you didn't specify what happens otherwise (i.e. with prob $1 - \pi_j$). – antkam Oct 13 '19 at 13:15
• You are right. w.p is to be read as "with probability". You are right, the problem seems to be not mathematically precise. You can imagine discarding as setting $Y_i= \mu_X$ w.p $1-\pi_j$. – lebesgue Oct 13 '19 at 14:21
• ah, "discarded"... so there are actually 2 different models: (1) $Y_i = E[X]$ w.p. $1-\pi_j$ is the easier thing to describe. or (2) take the sequence $X_i$ and skip those where a separate Bernoulli trial comes up with fallure ($1 - \pi_j$). In this case, it might be $Y_1 = X_1$ then skip $X_2, X_3$ and $Y_2 = X_4$ etc? – antkam Oct 13 '19 at 15:28
• Yes you are right. Model 2 is what was intended to my original problem stated. Skipping those samples which give a failure in a n independent Bernoulli draw. – lebesgue Oct 13 '19 at 16:24

Partial answer only (too tired to calculate the Variance)

Model $$2$$ is actually easier to solve, so I will start with that.

Let $$A$$ be the event that a particular $$X_i$$ is accepted, i.e. its associated $$Bernoulli(\pi_j)$$ results in success.

What is $$Y_k$$? It is the next accepted $$X_i$$ trial. In other words it is distributed like $$X_i$$ but conditioned on $$X_i$$ being accepted. If I may abuse notation a bit, $$Y \sim X \mid A$$. This means we use Bayes Theorem for the case when one variable is continuous and obtain the pdf of $$Y$$:

$$\forall z \in J_j: f_Y(z) = f_{X\mid A}(z) = {f_X(z) P(A \mid X=z) \over P(A)} = {\frac12 \pi_j \over P(A)}$$

IMHO it is easier to understand this geometrically:

• $$X$$ is uniform in $$(-1, 1)$$ so its pdf $$f_X(x) = \frac12$$, a constant function (within the support or range of $$X$$, i.e. $$(-1,1)$$).

• Now you obtain the pdf of $$Y$$ in two steps:

• (Step 1) chop up $$(-1,1)$$ into the $$J_1, \dots J_{m+1}$$ regions, and in each region $$J_j$$, rescale the $$f_X$$ by $$\pi_j$$. This gives you a piecewise-constant function.

• (Step 2) the piecewise-constant function does not have area $$=1$$, so we need to rescale the whole thing by whatever constant necessary to make its area $$=1$$ s.t. it can be the pdf of $$Y$$. That rescaling constant is $${1 \over P(A)}$$.

The main equation above gives the full pdf of $$Y$$, so you can calculate $$E[Y]$$ and $$Var[Y]$$ from it. The rest is just algebraic manipulations. First:

$$P(A) = \int_{-1}^1 P(A \mid X=z) f_X(z) dz = \frac12 \sum_{j=1}^{m+1} length(J_j) \pi_j$$

where $$length(J_j) = b_j - b_{j-1}$$ is the length of that interval which intersects with the range of $$X$$. (We use the convention $$b_0 = -1, b_{m+1} = 1$$.) The above $$P(A)$$ simply calculates the area of the piecewise-constant function after Step 1, s.t. when you rescale it by $${1 \over P(A)}$$ you get the same "shape" but now with area $$=1$$ and it can be a pdf.

Next we go for $$E[Y]$$. We can write the integral, but in this specific problem it is easier to condition on the region, because conditioned on $$Y$$ (i.e. the original $$X$$) being $$\in J_j$$, it is uniform within that region and therefore has mean equal to the center of the region. So:

$$E[Y] =\sum_{j=1}^{m+1} E[Y \mid Y\in J_j] P(Y \in J_j) =\sum_{j=1}^{m+1} {b_j + b_{j-1} \over 2} {\frac12 length(J_j) \pi_j \over P(A)}$$

Sorry I am too tired to explicitly calculate $$Var[Y]$$ right now...

A note on Model 1, where $$Y = \mu_X$$ if the Bernoulli trial fails. This actually makes $$Y$$ a mixture of a continuous variable (when $$X$$ is accepted) and a discrete variable (when $$X$$ is discarded). In a typical continuous variable (like $$X$$), the probability density is well defined and $$>0$$ in the range of $$X$$, but the probability that the variable equals any particular value is zero, e.g. $$\forall x: P(X=x) = 0$$. In the mixed $$Y$$, we have $$P(Y = \mu_X) = 1- P(A) \neq 0$$.

It's not a problem, really, but just makes everything tedious. E.g. let $$\bar{A}$$ denote the complement of $$A$$, i.e. $$X$$ is discarded, then:

$$\begin{array}{} E[Y] &=E[Y \mid \bar{A}] P(\bar{A}) + \sum_{j=1}^{m+1} E[Y \mid (X\in J_j) \cap A ] P((X \in J_j) \cap A)\\ &= \mu_X(1-P(A)) + \sum_{j=1}^{m+1} {b_j + b_{j-1} \over 2} \frac12 length(J_j) \pi_j \end{array}$$