Could the simple iterative algorithm for generating random numbers be favouring digits of value 0

algorithmic-randomnessalgorithmscomputer sciencerandomrandom variables

I have been creating a random number generator, using Python.

My script starts by creating a one-dimensional $1 \times N$ matrix of integers ranging from 0 and 9. We can denote this matrix as $A$. I am using a novel method for this, which arguably involves pure random numbers. So far this seems to be successful and an analysis of 4000 of these integers shows that each integer has approximately ~10% representation in $A$, which indicates (although does not fully prove) that my random 0-9 integer generator is fair. So let us assume that the integers are randomly generated and distributed.

Next, I created an iterative algorithm that attempts to add an extra layer of chaos and generate random numbers between 0 to 1, using the elements of $A$. The algorithm is as follows:

Start with a pre-specified integer seed $n$.
Define $N_1 = A_{1,n}$ and use $N_1$ to create the first decimal place of a new number between 0 and 1 (so starting with "0.$N_1$").
Define $S_1 = A_{1,n+1}$ and use $S_1$ as a step size to the next desired element: $A_{1,n+1+S_1}$.
Define $N_2 = A_{1,n+1+S_1}$ and use $N_2$ to create the second decimal place.
Define $S_2 = A_{1,n+2+S_1}$ and use $S_2$ as the next step size.
Define $N_3 = A_{1,n+S_2}$ as the third decimal place and so on.
Finish building the number after 10 decimal places, then begin building a new number (resuming from the same element of $A$).
If the $j$ value of any desired $A_{1,j}$ element at any time exceeds $N$ (which would result in an error), then redefine that $j$ as: $j_{new} = j_{old} – N$ and continue.

In plain English, this can be described as: "Make the first decimal place the $n^{th}$ element (where $n$ is the seed), skip ahead to the next element, skip ahead by the value of that element, make the second decimal place the new element you reach, skip ahead to the next element from there and so on. Wrap around to the beginning of the matrix if the final element is surpassed. Start a new number after 10 decimal places have been generated."

The result of this is a set of supposedly random numbers, between 0 to 1, each having been constructed as: $0. N_1 N_2 N_3 … N_{10}$.

In order to gain some insight into how random the generated set is, I created a script that monitors the abundances of individual 0-9 digits in the decimal places. It turns out that consistently ~20% of the digits are now 0, which is a disproportionately high representation (and I made sure to leave out the "0." starts of the numbers!).

Given that there was a previously balanced representation of 0 digits in the integer matrix $A$, what about the algorithm could have favoured 0 digits ending up in the final random number set? How can such a question be mathematically approached? If there is no clear bias induced by the algorithm I have described, then does it indicate that the integers of $A$ were perhaps not randomly distributed, after all?

Best Answer

Yes, your algorithm is generating "extra" $0$s because it makes $S_i$ serve as the "skip distance" to the next digit. So, if $S_i=0$ there is no skip, forcing the next digit to be $0:$ $$S_1 = A_{1,n+1}=0\quad\implies\quad N_2 = A_{1,n+1+S_1}=A_{1,n+1}=S_1=0,$$ etc. Thus, $$\begin{align} P(N_{i+1}=0)&=\underbrace{P(N_{i+1}=0\mid S_i=0)}_{1}\underbrace{P(S_i=0)}_{0.1}+\underbrace{P(N_{i+1}=0\mid S_i\ne 0)}_{0.1}\underbrace{P(S_i\ne 0)}_{0.9}=0.19.\\ \end{align}$$

(I also verified this by generating $10^8$ uniformly distributed random digits and applying your algorithm.)

Note: I think there's a typo where you wrote $N_3 = A_{1,n+S_2}$, which presumably should be $N_3 = A_{1,(n+2+S_1)+S_2};$ i.e., the position of $N_3$ is meant to be $S_2$ places after the position of $S_2,$ etc., so that $S_i$ acts as a "skip distance".

The non-uniformity in the digit-distribution is easily removed by simply adding $1$ to the present "skip distance"; i.e., use $N_2 = A_{1,n+1+(S_1+1)},\ N_3 = A_{1,(n+2+S_1)+(S_2+1)},$ etc.

Related Solutions

[Math] Algorithm for tetration to work with floating point numbers

When I studied the various known matrices of combinatorical numbers I also looked at the following idea: what if we discuss functions with a set of results, not only one number? So for instance the concept of sine and cosine gets some special charme if we look not only at f(x) = sin(x), g(x)=cos(x) but at 2x2-matrices containing cos() and sin() and the input of matrix-multiplications has two parameters and the output as well. This helps much to understand algebra with complex numbers, for instance.
After that I thought: what if we look at matrices which allow not only to look at a function $f(x)=a_0 + a_1x + a_2x^2 + a_3x^3 + ...$ involving all the powers of $x$ as input, but also spit out all powers of $f(x)$ in the same manner. Then for instance

$$ V(x) \cdot P = V(x+1) $$ where $V(x)$ means the vector $V(x)=[1,x,x^2,x^3,x^4,...]$ and $P$ means the upper triangular Pascal matrix. Of course $$(V(x) \cdot P )\cdot P = V(x+2)$$ and of course it is easily to see that $$V(x) \cdot P^h = V(x+h)$$ first for integral powers of $P$ which any programmer can implement - and if you download Pari/GP you can just use the matrix-language of it.
It needs a certain ingenious mind to find out how to define fractional powers of $P$ to make the fractional operation of addition possible...

After that one can find, that a very well known combinatorical matrix $B$ performs $$ V(x) \cdot B = V(exp(x))$$
and logically $$ V(x) \cdot B^h = V(exp^{oh}(x))$$
for integer $h$ - but which means implementing tetration to integer iterates $h$.
To find now the fractional-h power of $B$ -which is your question- is not so simple, and using matrices practically means to truncate them to finite size.

Well, anyway meaningful approximations can be done and truncated versions of $B$ can be diagonalized and general powers of it can then be found by simply computing the general powers on each of the scalar numbers in the diagonal - and so one can introduce fractional powers of $B$ in the formula $$ V_{32}(x) \cdot B_{32x32}^h \approx V_{32}(exp^{oh}(x))$$
(for for instance $32$ and $32x32$ sized vectors and matrix) and this implements then the fractional iteration of the exponential, aka tetration.

The matrices $P$ and $B$ of which I talk here are known as "Carleman-matrices" and this whole concept can be applied to fractional iteration of many functions for which you have a power series. (Well, we'll have convergence issues and much more complication and non-uniqueness and what not - but that's not the focus of this answer)
So using "mateigen" in Pari/GP for the diagonalization allows the following Pari/GP code:

n=32
V(x,dim=n) = vector(dim,r,x^{r-1})      \\ define the function for generating V(x)
b = sqrt(2)
bl = log(b)
B = matrix(n,n,r,c,bl^(r-1)*(c-1)^(r-1)/(r-1)!)  \\ define the Carlemanmatrix 
    \\ for   V(x) * B = V(b^x) 
Y = V(1) * B  \\ in Y we'll have approximately V(b^x) 
    \\where the first two entries have the best approximation to the expected
    \\ value
    \\ now diagonalize
M = mateigen(B) \\ needs high precision of say 800 or 1600 digits for all operations
W = M^-1 
D = diag(W*B*M) 
   \\ after this we can do B^h = M * D^h * W 
   \\  and  D^h can be done by d_k,k^h on the elements of the diagonal
dpow(A,h=1)= for (r=1,#A, A[r,r]=A[r,r]^h);return(A) \\ define a function for powers of D      

BPow (h)= M*dpow(D,h)*W  \\ define a function for the fractional power of B

\\ after this we can do
Y = V(1) * BPow(0.5) 
\\ and have in Y[2] a nice approximation for the half-iterate or b^^0.5
\\where b=sqrt(2)

I you want to experiment with this I recommend to use dimension $n=16$ first (with much poorer approximation) but already you'll need default(realprecision,200) or so. To do mateigen on matrix $B$ of size $32x32$ one must compute $b$, $bl$, $B$ and so on already to even $1600$ or $2000$ internal digits precision and the mateigen can need several minutes to complete. I tried this also with size $n=64 $ and this was really a full afternoon computation... but is a nice exercise anyway.

Remark: I've discussed this a bit in the tetrationforum and have made an article which compares a couple of such naive approaches to fractional exponential-towers. The above method I've called there "polynomial approach" because using finite matrices without taking care for compatibility with the theoretically infinite sized matrices and their fractional powers means just involve polynomials for that approximative solutions. See the index and the link to the essay

Best Answer

Related Solutions

[Math] Algorithm for tetration to work with floating point numbers

Related Question