How to create a sample of a desired size with rejection sampling from a custom distribution?

Question

I have implemented rejection sampling on a customized distribution in R. The code appears to work, however due to the nature of rejection sampling, there is a certain proportion of values that get rejected so I was wondering how I could create a sample of a desired size?

By this I mean the following. Suppose I desire a sample of size $n$. When I perform the rejection sampling procedure on $n$ values, I end up with a sample of size $m$ where $m < n$ from my custom distribution. I could try running the procedure $N$ times where $N > n$, but then I would end up with a sample of size $M$ where $M > n$ (for a choice of sufficiently large $N$). If I end up with $m$ points, I have fewer points than I need. If I end up with $M$ points, then I have more points than I need and I don't know which points to eliminate.

Can someone advise me on how to solve this please?

Thank you for looking at my post.

Answer 1

When the rejection rate is $q$ and you make $N$ attempts, the number of values you generate has a Binomial distribution with parameters $N$ and $p=1-q.$

Usually $N$ is so large that the Normal approximation to this distribution is a good description. Its mean is $\mu = pN$ -- that's how many points you expect to get -- and its standard deviation is $\sigma = \sqrt{Np(1-p)}$ -- that's the amount by which you expect the actual number to differ from $Np.$

Because deviations of more than a few standard deviations are extremely rare in a Normal distribution, increase $N$ so that $\mu-Z\sigma$ exceeds the desired count $n$. Here, $Z$ is a "few." That is, choose $N$ to satisfy

$$pN - Z\sqrt{Np(1-p)} \ge n.$$

For sufficiently large $Z$ you are extremely likely to generate too many values: just throw the last extra ones away.

Because this value of $N$ will be relatively only a little greater than $n/p,$ approximate the solution by

$$N \approx \frac{n}{p} + Z\sqrt{\frac{n}{p}\,p(1-p)}=\frac{n}{p} + Z\sqrt{n(1-p)}.$$

Since the additional amount $Z\sqrt{n(1-p)}$ is on the order of $n^{1/2},$ which is small compared to $n$ (for $n$ large), you can afford to make $Z$ fairly big. If the Binomial distribution were truly Normal, $Z=6$ or so would suffice. But the Normal approximation doesn't work too well in the tails, so use a larger value of $Z,$ perhaps $Z=16$ or even $Z=20.$

---

For example, suppose two-thirds of the attempts are going to be rejected on average, so that $q=2/3$ and $p=1-2/3 = 1/3.$ If you would like $n=1000$ values, run

$$N = \frac{1000}{1/3} + 16\sqrt{1000(2/3)} = 3414$$

trials. The relative cost of this insurance is $414/3000$ or an increase of about $14\%$ over the expected cost of $3000.$ (The cost is negligible for very large samples, such as $n=10^6.$

Before committing to a particular value of $Z$, you can check the true Binomial probability that this algorithm will fail to generate at least $n$ values. In the example, the chance is less than $0.2$ per million. In other words, if you generate $3414$ trials, you expect to end up with fewer than $1000$ non-rejected values only once in every five million times you do this. (For added protection, you can verify the number of non-rejected values is adequate and, if not, loop back to generate the very few more that are needed.)

Answer 2

I have more points than I need and I don't know which points to eliminate.

It doesn't matter! They are all i.i.d., so any unbiased rule (not based on looking at the values themselves) is fine -- e.g., throw away the ones generated last, or first. You could even generate random indexes to throw away but there's no need to bother. All of these give the same statistical properties for the resulting size-$n$ sample.

whuber's answer notes this -- and also does a lot of other analysis to estimate how large an initial sample is needed to generate $n$ points -- however, I don't think that's the core problem you had. You can always generate more points if you don't get enough in the first sample. It seems the only thing you were really "missing" is that you can freely throw away points when you have too many.