## Strange random number generator

A random number generator outputs a sequence of n random numbers {r1, r2, … , rn} from a uniform distribution on (0,1). This sequence can be used to define an integer p with n-1 binary digits according to the following algorithm:

     for (i in 1:(n-1))
{
if (r[i+1] > r[i])
then set bit i of p equal to 1
else set bit i of p equal to 0
}


Let’s call p the pattern of the sequence. For example, the sequence {0.5703, 0.1617, 0.2629, 0.5404, 0.3860} has pattern 0110; the sequence {0.6441, 0.6715, 0.0802, 0.8980, 0.5447, 0.4748, 0.5214, 0.0110} has pattern 1010010. Since a set of size n corresponds to a pattern of n-1 binary digits, the number of possible patterns is 2n-1.

What does the distribution of p look like? For a fixed n, the probability that a set of n random numbers will correspond to a particular pattern p is given by:

$(-1)^w \int_0^1 \int_{p_1}^{\tau_1} \int_{p_2}^{\tau_2} \cdots \int_{p_{n-1}}^{\tau_{n-1}} d \tau_n d \tau_{n-1} \cdots d \tau_2 d \tau_1$

where w is the Hamming weight of p (i.e., the number of 1s it contains).

Here is some Mathematica code to generate the distribution for a given value of n:

n = -1;
While[(n < 1) || (n > 12),
n = Input["Enter number of binary digits: "]];
func = {1};
ProgressIndicator[Dynamic[pr], {1, n}]
For[i = 1, i ≤ n, i++,
pr = i;
fl = Integrate[func, {x, 0, x}];
fr = Integrate[func, {x, x, 1}];
func = Join[fl, fr];];
freq = Integrate[func, {x, 0, 1}];
BarChart[freq]


Here is a barchart of the frequency distribution with n=6. The x-axis runs from 00000 to 11111.

The distribution is self-similar, which I guess makes it a fractal.