Sampling from a distribution with arithmetic coding

Question

I'm struggling to implement a sampling scheme where arithmetic decoding is used to generate the samples.

This is based on excercise 6.3 from David MacKay's book, where you compare rejection sampling with arithmetic decoding sampling. For each symbol you want to generate with rejection sampling you have to generate an (uniformly distributed) integer from $1$ to $2^{32}$, rescale it to be in $\left(0,1\right)$ and then if your random number is less than 0.99 emit a 0 (or 1 accordingly).

To generate symbols by arithmetic coding, you feed standard random bits into an arithmetic decoder equipped with the appropriate model.

(Note: I think this works because when we compress data we go from symbols generated under P to symbols under Q - P,Q are probability distributions - and decompression works the opposite way. If we choose Q to be an easy distribution to sample from then by decoding symbols from Q we can generate symbols from P.)

I've implemented the methods in the code below:

clear all;
pdist = [0.99 0.01] ;
counts = [99 1];
N=100;
U = randi(2^32, [1 N]);
U = U/2^32;
rejection_samples(1:N) = 0;

for i=1:N
    if U(i) < pdist(1) ;
       rejection_samples(i) = 0;
    else
       rejection_samples(i) = 1;
    end
end

seq = randsrc(1, N, [0 1;0.5 0.5]);
arith_samples = arithdeco(seq, counts, N);

I'm confused becuase the arith_samples variable sometimes contains a 2 (but I guess this just maps onto (0,1) as the output never seems to contain a 0 - just 1 and 2).

The answers say that the rejection sampling method takes 32 bits per generated symbol, but that arithmetic coding only takes 0.081 = $H_2\left(0.01\right)$. I can't really see that from my code.

Is there some way I could measure 'bits used'?

Also this seems a bit too good to be true - surely there are better ways of generating random numbers than just rejection sampling?

The relevant chapter is found here:

http://www.cs.toronto.edu/~mackay/itprnn/ps/105.124.pdf

Answer 1

While this answer may have come a bit too late for OP, here it is.

Let's start by replacing the phrase 'bits used' with 'queries to an unbiased bit generator'. Then assume that the only random number generator you have at your disposal is an unbiased bit generator, that generates $0$ or $1$ with probabilities $p_0=p_1=0.5$. e.g. something like the following (well, this is probably not the best way to implement an unbiased bit generator in octave/matlab, but suffices for our purpose):

function b = unbiased_bit_generator() b=randi([0, 1], 1,1); end;

Now your goal is to use this function as less as possible (1 use = 1 query) to generate a given number of iid biased bits with a specified bias.

In case of rejection sampling, you would use $\texttt{unbiased_bit_generator()}$ 32 times to get a 32-bit integer $i$, then check if $\frac{i}{2^{32}}<p_0$ to generate one biased bit. Thus you require $32$ queries/biased-bit. Here is a biased bit generator based on 32-bit rejection sampling:

function b=biased_bit_generator1(p0)
    for i=1:32
        bin_array(i)=unbiased_bit_generator();
    end;
    str=num2str(bin_array);
    str=str(str~=' ');
    if bin2dec(str)/2^32<p0
        b=0;
    else
        b=1;
    end;
end;

In case of arithmetic decoder based sampling, you feed the decoder with a sequence $s$ of $k$ unbiased bits (generated using $\texttt{unbiased_bit_generator()}$ $k$ times) and get a sequence $d$ of of length $n$ biased bits. Thus you needed $\frac{k}{n}$ queries/biased-bit. Because arithmetic coding is a variable-length coding technique, two encoded sequences of same length may correspond to two message sequences of different length. Thus, if you were to end up with sequence $s_k'\neq s_k$ after using the unbiased bit generator $k$ times and feed it to the arithmetic decoder, then you may get as output a sequence $c_{n'}$ s.t. $n\neq n'$. In this case, you needed $\frac{k}{n'}$ queries/biased-bit.

The number 0.081 specified in the book is the expected number of queries/biased-bit where the expectation is over all encoded input sequences to the decoder. As you keep feeding your decoder with more and more unbiased bits and keep generating more and more biased bits the queries/biased-bit will approach this number.

Also this seems a bit too good to be true

I assume you are talking about 0.081 queries/biased-bit. I think this statement is equivalent to saying that Shannon's source coding theorem is too good to be true ;)