%% Random Numbers
%% MATLAB's Built-in Functions
% x = rand: single uniformly distributed random number in (0,1).
% x = rand(n): n by n matrix ...
% x = rand(n,m): n by m matrix

clear
R=rand(10000,1);
histogram(R)

%
% Generate random integers in a range.

randi([10,20],2,3)

%
% Normally distributed random numbers (\mu=0, \sigma=1).

R=randn(1000,1);
histogram(R)

%
% Change the desired mean and standard deviation.

mu  = 5;
sig = 2;
R = mu + sig*randn(1000,1);
histogram(R)
[mean(R), std(R)]

%
% Random numbers are not that random. Need to change the seed.

rng(120)
rand()
rand()

rng(120)
rand()
rand()

rng(50)
rand()
rand()

rng('shuffle')
rand()
rand()

%
% One can create a random permutation of numbers.

randperm(10)

%
%% Middle Square Method
% How are random numbers generated? They can only be generated using a deterministic process in a computer unless physically determined from random physical processes. Thus, one inputs a seed, does a computation, and outputs a new seed to generate the next number.
% One such process is the middle square method. One specifies a number of a given length, n. Square the seed and extract a middle string of n digits of the given length. If the number has fewer digits, then pad it on the left with 0's. This is the new seed. The process is repeated. The list of outputs are the pseudorandom numbers.

clear
seed = 1234;
seed^2
n = 4;                    % Length
lmt = 1000;               % Number of prns 

for i=1:lmt
  temp=num2str(seed^2);
  s=size(temp);
  while strlength(temp) <2*n
      temp = "0"+temp;
  end
  s1=floor(n/2)+1;
  s4=s1+n-1;
  x(i)=extractBetween(temp,s1,s4);
  %fprintf('%d %s\n',i,x(i));
  seed=str2num(x(i));
  y(i)=seed;                          % Could normalize by 9999
end
idx = y > 0;
histogram(y(idx))

% The histogram indicates how random this is.
%
% One can take the process and create a function.

seed = 1234;
seedmax = 10^numel(num2str(seed))-1;
n = 4;
Num = 10;
for i=1:Num
    seed = RandMidSquare(seed,n)
    if seed==0
        i
        break
    end
end

%
% Then, one can create a grayscale picture and look for patterns or periodicities in the pattern.

seed = 1234567890;
n = 5;
seedmax = 10^n-1;
nrows=100;
ncols=100;
M=zeros(nrows,ncols);
for i=1:nrows
    for j=1:ncols
        seed=RandMidSquare(seed,n);
        r = seed/seedmax;
        M(i,j)=r*255;
    end
end
imagesc(M)
colormap(gray)
set(gca,'XTick',[], 'YTick', [])

%
%% Linear Congruential Generator
% Another method for generating pseudorandom numbers is the Linear Congruential Generator. This generates seeds using a simple formula: Let s be the seed. The next seed is given by s = a s + c \mbox{ mod } m, where a is the multiplier, c the increment, and m the modulus.
% There have been studies producing good values for the constants. Here is an example with a histogram. [See Wikipedia entry]

a=1103515245;
c=12345;
m=2^31;
seed=1234;
N = 10000;
for k=1:N
    seed=nextRand(seed,a,c,m);
    r(k) = seed/m;
end
histogram(r)

%
% We can also produce a grayscale image showing there is a hidden pattern.

nrows=500;
ncols=500;                         % Try changing to 550
a=1103515245;
c=12345;
m=2^31;
seed=1234;
M=zeros(nrows,ncols);
for i=1:nrows
    for j=1:ncols
        seed=nextRand(seed,a,c,m);
        r(k) = seed/m;
        M(i,j)=r(k)*255;
    end
end
imagesc(M)
colormap(gray)
set(gca,'XTick',[], 'YTick', [])

%
%% Mersenne Twister
% This begs the question - How well does MATLAB's random number generator fair? What type of generator does MATLAB use?
% MATLAB's default is the Mersenne Twister generator with seed 0. The Mersenne twister was developed in 1997.
% This PRNG has a period length is chosen to be a Mersenne prime. Many software programs use {\displaystyle 2^{19937}-1}. [See more at the links.] This was the 24th known Mersenne prime consisting of 6002 digits and was found in 1971.

nrows=1000;
ncols=1000;
M=zeros(nrows,ncols);
for i=1:nrows
    for j=1:ncols
        M(i,j)=rand()*255;
    end
end
imagesc(M)
colormap(gray)
set(gca,'XTick',[], 'YTick', [])

% How good are your generators? One could us the diehard tests, statistical tests for the quality of a PRNG.
%
%% Mean and Variance
%  Using the uniform random numbers from the last example, we can create a set of random numbers with a given mean and variance. This is based on the fact that rand generates numbers uniformly distributed on the interval (0,1) with mean 0.5 and variance 1/12.

figure
M2=reshape(M/255,1,[]);
m = 2;
v = .1;
M3=sqrt(12*v) * (M2-0.5) + m;
histogram(M3)
axis tight
disp(['mean: ', num2str(mean(M3))])
disp(['variance: ', num2str(var(M3))])
disp(['st. deviation: ', num2str(std(M3))])

%
%% The Box-Muller Method
% One often wants a different distribution, such as a normal or Gaussian distribution. This is the typical bell-shaped distribution with a given mean and standard deviation.
% One can generate a normal distribution from a uniform distribution. One algorithm is given by the Box-Mueller Method.

% BOX-MUELLER METHOD
u = rand(2,100000);
r = sqrt(-2*log(u(1,:)));
theta = 2*pi*u(2,:);
x = r.*cos(theta);
y = r.*sin(theta);
 
figure
% X SAMPLES
subplot(121);
hist(x,100);
axis square
title(sprintf('Box-Muller Samples x\n Mean = %1.2f\n Variance = %1.2f',mean(x),var(x)))
xlim([-6 6])
 
% Y SAMPLES
subplot(122);
hist(y,100);
axis square
title(sprintf('Box-Muller Samples y\n Mean = %1.2f\n Variance = %1.2f',mean(y),var(y)))
xlim([-6 6])

%
% The Box-Muller Method gives a distribution of mean 0 and variance 1. For distributions of other shapes, one needs to rescale and shift the output. If we have a normal distribution X\sim N(0,1),  then we can create a new distribution,Y=\sigma X+\mu \sim N(\mu ,\sigma ^2), with mean \mu and variance \sigma^2.

m = 2.0;
s = 0.3;

x2 = s*x+m;
y2 = s*y+m;

figure
% X SAMPLES
subplot(121);
hist(x2,100);
axis square
title(sprintf('Box-Muller Samples x\n Mean = %1.2f\n St. Dev. = %1.2f',mean(x2),std(x2)))
xlim([-6 6])
 
% Y SAMPLES
subplot(122);
hist(y2,100);
axis square
title(sprintf('Box-Muller Samples y\n Mean = %1.2f\n St. Dev. = %1.2f',mean(y2),std(y2)))
xlim([-6 6])

%
%% Blum Blum Shub Algorithm
% This is yet one other PRNG called the Blum Blum Shub Algorithm. This still should be looked at visually similar to other PRNGs.

clear
p=11; q=19;
m=p*q;

% The seed should be an integer co-prime
% to m and not 1 or 0.
seedg=5;
N=10;
for k=1:N
    rn(k)=mod(seedg^2,m);
    seedg=rn(k);
end;
%The Random Numbers are
rn/m

%
%% Monte Carlo Method
% One can use random numbers to compute \pi.We randomly toss darts at a 2\times 2square. The area of the inscribed unit circle is \pi, Since the area of the unit circle is the fraction of darts landing in the circle times the area of the circumscribed square, we get an estimae for \pi as 4 times the fraction of darts.

clear;
N=100; % Square 100 pixels by 100 pixels.
r=1;   % Uses N^2 random numbers
n=0;  
M=255*ones(N);
for i=1:N
    for j=1:N
        x=-r+2*r*rand();
        y=-r+2*r*rand();
        if ((x^2+y^2)<=r^2)
            n = n+1;
            M(floor(x*N/2)+N/2+1,floor(y*N/2)+N/2+1)=0;
        end
    end
end

imagesc(M)
colormap(gray)
set(gca,'XTick',[], 'YTick', [])
rectangle('Position',[0,1,N,N-1])
axis equal
axis off

pi_sim=4*n/N^2

%
%% RANDU PRNG
% IBM's RANDU is considered to be one of the most ill-conceived random number generators. It was used in the 1960s and 1970s. It has a multiplier of 65539 and modulus of 2^{31}.For more, see https://en.wikipedia.org/wiki/RANDU.

clear
a=65539;
c=0;
m=2^31;
seed=123;
nrows=1000;
ncols=1000;
M=zeros(nrows,ncols);
for i=1:nrows
    for j=1:ncols
        seed=nextRand(seed,a,c,m);
        M(i,j)=seed/m*255;
    end
end
imagesc(M)
colormap(gray)
set(gca,'XTick',[], 'YTick', [])

%
% Lets try a 3D plot. You can rotate the plot and look for patterns. Find any?

clear
a=65539;
c=0;
m=2^31;
seed=123;
nrows=20;
ncols=20;
ndeep=21;
M3=zeros(nrows,ncols,ndeep);
for i=1:nrows
    for j=1:ncols
        for k=1:ndeep
            seed=nextRand(seed,a,c,m);
            M3(i,j,k)=seed/m*255;
        end
    end
end

X=reshape(M3,[],3);
plot3(X(:,1),X(:,2),X(:,3),'.');

%
% Now try a plot of (x,y,z) where each coordinate is randomly generated. How does this plot look?

clear
a=65539;
c=0;
m=2^31;
seed=1233;
N=1000;
X=zeros(1,N);
Y=X;
Z=X;
for i=1:N
    seed=nextRand(seed,a,c,m);
    X(i)=seed/m;
    seed=nextRand(seed,a,c,m);
    Y(i)=seed/m;
    seed=nextRand(seed,a,c,m);
    Z(i)=seed/m;
end
plot3(X,Y,Z,'.')
view(-10,20)

%
% Finally, we can change the viewing angle. We see a definite pattern. It is known that LCGs with modulus m can generate points in n- dimensional space that lie in parallel hyperplanes.

view(-22.7200,-65.2084)

% Try other seeds. Maybe we can apply this 3D perspective with other PRNGs to see why some are better than others.
%
% This is typical of Linear Conguential Generators as pointed out in 1968 by Marsaglia, "Random Numbers Fall Mainly in Planes."
% For RANDU we have x_{k+1}=(2^{16}+3)x_k \mbox{ mod } 2^{31}. Then, iterating, x_{k+2}=(2^{16}+3)^2x_k \mbox{ mod } 2^{31}, or x_{k+2}-6x_{k+1}+9x_k = 0 \mbox{ mod } 2^{31}. This means x_{k+2}-6x_{k+1}+9x_k =  2^{31} \ell.
% Dividing by 2^{31}, we obtain numbers between 0 and 1. Calling the new numbers x,y,z, we obtain the planes z-6y+9x=  \ell.

% Plot plane 9x - 6y+z = 0 mod m
% Ax + By + Cz + D = 0
ell=1;
A = 9; B = -6; C = 1; D=0;
[x y] = meshgrid(0.2:0.1:0.7); 
%[x y] = meshgrid((-ell-1)/9:0.01:ell/9,-ell/6:.01:(1-ell)/6); 
z = -1/C*(A*x + B*y + D); 
surf(x,y,z) 

%
% Plotting the plane with the "random" numbers, we can see the plane lines up with a subset of points after a 3D rotation.

hold on
plot3(X,Y,Z,'.')
hold off
%view(-22.7200,-65.2084)

%
% This notebook's functions are below.

function y = nextRand(seed,a,c,m)
y = mod(a*seed+c , m);
end

function y = RandMidSquare(seed,n)
temp=num2str(seed^2);
s=size(temp);
while strlength(temp) <2*n
    temp = "0"+temp;
end
s1=floor(n/2)+1;
s4=s1+n-1;
x=string(extractBetween(temp,s1,s4));
y=str2num(x);   
end