Moments of Rejection Sampling

112 Views Asked by At

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<b_2\dots<b_m$ 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$

1

There are 1 best solutions below

0
On BEST ANSWER

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