# Accept-reject and subsets of iid samples

I have some confusion about subsets of iid samples being distributed as the original sample. As an illustration, consider the accept-reject algorithm to produce iid samples from a pdf $$f(x)$$. We draw, say $$M=100$$, samples from a proposal $$g(x)$$, perform a test for each of these individual samples and say we end up with $$N=30$$ samples that passed the test, and hence they represent an iid sample from $$f(x)$$.

Now, my question is: these $$N=30$$ samples are iid from $$f(x)$$ since they resulted from the accept-reject algorithm, but at the same time, are they iid from $$g(x)$$ since they were all drawn independently from it?

More generally, is any subset of an iid sample (from $$g(x)$$), also an iid sample (from $$g(x)$$)?

Your final 30 points are not an i.i.d. sample from $$g$$ (since they are an i.i.d. sample from $$f$$), the accept-reject step does modify the distribution. So the answer to your question is: no, a subsample from an i.i.d. sample with distribution $$g$$ does not have the same distribution $$g$$.

You can see that very easily on a simpler example. Suppose $$B_1, ..., B_n$$ are iid draw from a $$\mathcal{B}ernouilli(p)$$ (with $$p$$ not $$0$$ nor $$1$$), and select only the points which are equal to $$1$$. Your subsample is obviously not following a $$\mathcal{B}ernouilli(p)$$ since it is all equal to 1.

I guess what can be puzzling here is that if the subsample was draw uniformly and independently from the sample (in a bootstrap spirit for instance), then you would have that the distribution of your subsample is following the same distribution than the sample (without conditionning on the realization of the sample).

The accepted and rejected samples are no longer distributed from $$g$$ because an event depending on the realisations $$Y_i$$ occurred (acceptance with probability $$f(y_i)/Mg(y_i)$$) or did not occur (rejection with probability $$1-f(y_i)/Mg(y_i)$$) which, marginally, changed their distribution. While the accepted $$X_i$$'s are distributed from $$f$$ rather than $$g$$ and iid, since $$g(x_i)\times\frac{f(x_i)}{Mg(x_i)} \propto f(x_i)$$ the rejected $$Z_i$$'s are distributed from $$g(z_i)\times\left\{1-\frac{f(z_i)}{Mg(z_i)}\right\} \propto \frac{g(z_i)-Mf(z_i)}{1-M}$$ and independent conditional on $$N$$.

Here is an excerpt from Monte Carlo Statistical Methods taking advantage of the distinction between accepted and rejected subsamples. Itself borrowed from our 1996 Biometrika Rao-Blackwellisation paper.

Consider an Accept-Reject method based on the instrumental distribution with density $$g$$. If the original sample produced by the algorithm is $$(X_1,\ldots,X_m)$$, it can be associated with two iid samples, $$(U_1,\ldots,U_N)\quad\text{ and }\quad(Y_1,\ldots,Y_N)$$ with corresponding distributions $${\cal U}_{[0,1]}$$ and $$g$$; $$N$$ is then the stopping time associated with the acceptance of $$m$$ variables $$Y_j$$. An estimator of $$\mathbb E_f[h]$$ based on $$(X_1,\ldots,X_m)$$ can therefore be written $$\delta_1 = {1\over m} \; \sum_{i=1}^m \; h(X_i) = {1\over m}\; \sum_{j=1}^N\; h(Y_j)\; \mathbb I_{U_j\leq w_j}\,,$$ with $$w_j = f(Y_j)/Mg(Y_j).$$

A reduction of the variance of $$\delta_1$$ can be obtained by integrating out the $$U_i$$'s, which leads to the estimator $$\delta_2 = {1\over m} \; \sum_{j=1}^N \; \mathbb E[\mathbb I_{U_j \leq w_j} | N,Y_1,\ldots,Y_N] \; h(Y_j) = {1\over m} \sum_{i=1}^N \rho_i h(Y_i),$$ where, for $$i =1, \ldots, n-1$$, $$\rho_i$$ satisfies \begin{align*} \rho_i &= \mathbb{P}(U_i\le w_i|N=n,Y_1,\ldots,Y_n) \\ &= w_i \frac{\sum_{(i_1,\ldots,i_{m-2})} \prod_{j=1}^{m-2} w_{i_j} \prod_{j=m-1}^{n-2} (1-w_{i_j})}{\sum_{(i_1,\ldots,i_{m-1})} \prod_{j=1}^{m-1} w_{i_j} \prod_{j=m}^{n-1} (1-w_{i_j})}, \tag{1} \end{align*} while $$\rho_n = 1$$. The numerator sum is over all subsets of $$\{1,\ldots,i-1, i+1, \ldots, n-1\}$$ of size $$m-2$$, and the denominator sum is over all subsets of size $$m-1$$. The resulting estimator $$\delta_2$$ is an average over all the possible permutations of the realized sample, the permutations being weighted by their probabilities. The Rao-Blackwellized estimator is then a function only of $$(N,Y_{(1)},\ldots,Y_{(N-1)}, Y_N)$$, where $$Y_{(1)},\ldots,Y_{(N-1)}$$ are the order statistics.

Although the computation of the $$\rho_i$$'s may appear formidable, a recurrence relation of order $$n^2$$ can be used to calculate the estimator. Define, for $$k\le m < n$$, $$S_k(m) = \sum_{(i_1,\ldots,i_k)} \prod_{j=1}^{k} w_{i_j} \prod_{j=k+1}^{m} (1-w_{i_j}),$$ with $$\{i_1,\ldots,i_m\} = \{1,\ldots,m \}$$, $$S_k(m)=0$$ for $$k>m$$, and $$S^i_k(i)=S_k(i-1)$$. Then we can recursively calculate \begin{align*} S_k(m) &= w_mS_{k-1}(m-1)+(1-w_m)S_k(m-1), \\ S^i_k(m) &= w_m S^i_{k-1}(m-1)+(1-w_m)S^i_k(m-1) \end{align*} and note that weight $$\rho_i$$ of (1) is given by $$\rho_i =w_i\; S^i_{t-2}(n-1)\big/S_{t-1}(n-1) \qquad (i