It seems crazy to me that there's no way to produce unbiased numbers in an arbitrary range without rejection sampling and a loop. Is there a proof of this?
Double-precision Floats have more values between 0.0 and 1.0, than between 1.0 and 2.0. In fact, roughly half of ALL double-precision floats exist between -1.0 and 1.0, a very small minority of them exist between 1.0 and 2.0.
To generate unbiased random numbers between 0.0 and 2.0, it therefore requires you to either reject a significant amount of numbers in the 0.0 to 1.0 range, or perform some kind of many-to-few mapping in the 1.0 to 2.0 range.
----------
With regards to arbitrary INTEGER ranges, the proof is even easier. A random bitstream has 2^number-of-bits possible random values. Which does NOT divide evenly into an arbitrary integer range.
For example, 5-random bits will represent 32-different values. There's no way to map 32-values and divide them evenly into 0-9 (10 numbers).
I'd expect it's possible by changing the generator at the lowest level, but it makes sense to me that you need a loop if you don't control the underlying generator.
Imagine you want to turn a random number in 1..4 into a random number in 1..3. The original is your only source of randomness, so the rest should be deterministic. Then each outcome in 1..4 has to map to exactly one number in 1..3, but there's no mapping that accepts all of 1..4 while still giving each of 1..3 an equal probability.
Could it cap the maximum number of generator calls though? Rejection sampling is technically O(infinity) because you could reject an unbounded number of times. This isn't a problem in practice but it sure is theoretically annoying. With a cap on the maximum number of calls, it would be O(1).
If you want a number in 1..3, and the generator provides numbers from 1..4, and you want a cap n, then you could model it as a generator that provides numbers in 1..4^n. There's never a way to split that space into three equal parts.
You always end up with unwanted extra possibilities that you still need to handle somehow.
You don't need rejection sampling, but you do need a loop. It is easier to see that a finite number of samples is not sufficient:
If you have only a probability distribution defined by a product space where all distinguishable events have probability p^i, for some finite i, then any subset of the distinguishable events accumulate to a probability r * p^i for some integral r. If your goal is a probability that is not an integral multiple of p^i, you are out of luck with a finite number of samples.
Proposition: There is no way, given some integer r and a sequence of iid random integer samples X1,X2,... uniformly distributed in [0,n) for some n>1, to always construct some finite k and function f : [0,n)^k -> [0,r) such that f(X1,...,Xk) is uniformly distributed over [0,r).
Proof: Suppose such a k and f exist. The distribution of U[0,n)^k is isomorphic to that of U[0,n^k) (just treat it like writing down a k-digit number in base n). And so f must be a function from a set of n^k things to a set of r things. By the pigeonhole principle there must be some integers x,y in [0,r) such that the preimage of x has size at least ceiling(n^k/ r) and the preimage of y has size at most floor(n^k/ r). By the fundamental theorem of arithmetic there exists (for any n,k) some r such that n^k/r is not an integer and so the probabilities of x,y (being proportional to the sizes of their fibres under f) are not equal.
————
The gist of this is that you always might need to loop (repeatedly sample) for some ranges and you might need to repeat arbitrarily many times.
One way is by rejection.
Another nice thing is if you want to generate a Bernoulli r.v. for some probability p a computable number, you lazily compute the bits of p simultaneously to a random sequence of bits (distributed as Bernoulli(1/2)) and compare the two. If your random sequence is definitely less than p then generate 0. If definitely greater then generate 1. If not sure then generate some more bits.
In this way any Bernoulli random variable may be generated from an infinite sequence of iid Bernoulli(1/2), and basically any probability space can be modelled in this way too. In this sense, all of probability can be built out of tosses of a fair coin)
