Generate normal distribution with LGC random issue

hi as you can see im trying to compere values btween my function randi and built in function rand but i can't find a way to make myrandi plot normal distribution plot

clc
%N=input('sequence_length')
N=5000

N = 5000

 la=65539

la = 65539

MinSec = fix(clock); 
seed = 100*(MinSec(5) + MinSec(6));
Urand = myrandi(seed,N);  
  subplot(3,2,1)
  hist(Urand) 
  %% Task 1: : the Linear Congruential Generator (LCG) 
v = rand(1,N);
%u=u./max(u);
subplot(3,2,2);
hist(v,100);
title(subplot(3,2,2),'Uniform random: rand()');
%% Task 2: Generate exponential distribution with LGC random 
u=myrandi(seed,N);
e=-log(1-u)./la;
subplot(3,2,3);
hist(e,100);
title(subplot(3,2,3),'Exponential random: myrand()');

%% Task 2: Generate exponential distribution with rand() random 
u=rand(1,N);
f=-log(1-u)./la;
subplot(3,2,4);
hist(f,100);
title(subplot(3,2,4),'Exponential random: rand()');
%% Task 2: Generate normal distribution with LGC random 
u1 = myrandi(seed,N);
u2 = myrandi(seed,N);
a  = (-2*log(u1));
b  = 2*pi*u2;
n  = sqrt(a.*sin(b));
m  = real(sqrt(a).*cos(b));

M = mean(m,'all'); %mean
s = std(m); %standard deaviation
T = normrnd(M, s);
subplot(3,2,5);
hist(T);
title(subplot(3,2,5),'Normal random: myrand()');
%% Task 2: Generate normal distribution with rand() function
u1 = rand(1,N);
u2 = rand(1,N);
a  = (-2*log(u1)).^.5;
b  = 2*pi*u2;
p  = real(a.*sin(b));%p and q are both
q  = real(a.*cos(b));%normal random variables
subplot(3,2,6);
hist(q,100);
title(subplot(3,2,6),'Normal random: rand()');

%Function that produce an uniform random numbers
%Using the technic of LCG r(k) = [lambda*r(k-1)]modulo(P)
function Urand = myrandi(seed,sequence_length)

   la=65539;
    N=sequence_length;
    P=2^31;
%change N accordingly 
%part a 
Urand=zeros(1,N);
seed(1)=seed;
for i=2:N 
    Urand(i)= mod((la*seed*(i-1)),P);
    
end
Urand = Urand./(P-1);
end

Answer 1

[John D'Errico]

I might check that you REALLY are using the same expressions for both of the transformations to Normal. Have you done so? (STRONG HINT!) Just take a scan through that part of your codes. It popped out at me. I think you should see what you did wrong. (Sorry that I'm only giving you a hint here, but to be honest, it looks like your coding skills should be strong enough that you just missed something obvious. Once you look carefully, comparing the two code segments...)

[Faisal Al-Wazir]

you mean this u1 = myrandi(seed,N);

u2 = myrandi(seed,N)

yes i just reilezed it's not right

[Faisal Al-Wazir]

the issue is i don't know what to change in any of them because it's math

[John D'Errico]

Nope. I was wondering why might you take a sqrt in one instance, but not the other?

a = (-2*log(u1));

a = (-2*log(u1)).^.5;

If u1 is in both instances, a uniform pseudo-random variable, then won't that change the character of your result?

By the way, sqrt is a better idea then .^.5. It is far more readable at the very least, and readable code is important.

I've not looked at the specific algorithm to decide if it is correct to produce a Normal. I could do that if necessary, but this seems a significant difference.

[John D'Errico]

But now that I look more carefully at your code, I see this:

T = normrnd(M, s);
subplot(3,2,5);
hist(T);

Why, in the use of myrandi to generate your random numbers, do you then use normrnd, and plot a SCALAR variable as a histogram? :)

Again, look carefully at your code.

[Faisal Al-Wazir]

so what can change here to make it right ?

%% Task 2: Generate normal distribution with LGC random 
u1 = myrandi(seed,N);
u2 = myrandi(seed,N);
a  = (-2*log(u1));
b  = 2*pi*u2;
m  = real(sqrt(a).*cos(b));
M = mean(m,'all'); %mean
s = std(m); %standard deaviation
T = normrnd(M, s);
subplot(3,2,5);
hist(T);

[Faisal Al-Wazir]

well john i've been stuck in this problem for some time , so i try to do random things and see what will get me closer to the answer

[John D'Errico]

Looking more closely yet, your myrandi is a god-awfully bad random number generator. :) But that is a characteristic of those simple LCRNGs.

For example,

MinSec = fix(clock); 
seed = 100*(MinSec(5) + MinSec(6))

seed = 7700

N = 100000;

x = myrandi(seed,N);
histogram(x,100,'norm','pdf')

That is actually far too uniform looking. Do you see how consistently flat it is? But worse, there is a nice little pattern there. Instead, you would want to see more like random noise., like that below.

histogram(rand(1,N),100,'norm','pdf')

I thought I saw some strange behavior in the other plots too. Again, this is probably a characteristic of a bad random number generator. So not really your problem, as you are I assume looking at this RNG for a specific reason.

%Function that produce an uniform random numbers
%Using the technic of LCG r(k) = [lambda*r(k-1)]modulo(P)
function Urand = myrandi(seed,sequence_length)
la=65539;
N=sequence_length;
P=2^31;
%change N accordingly 
%part a 
Urand=zeros(1,N);
seed(1)=seed;
for i=2:N 
    Urand(i)= mod((la*seed*(i-1)),P);
    
end
Urand = Urand./(P-1);
end

[Faisal Al-Wazir]

thanks for the cretique, but i'm still not clear on what to do for the 5th plot

[John D'Errico]

That was merely an observation, as to why people stopped using simple LCRNGs like that many, many years ago. They are terrible for multiple reasons. I suppose, possibly, this may be why you are having a problem. But let me now look carefully at the algorithm, making sure you are doing the correct thing in that code.

A quick check shows you are using the Box-Muller method to generate normal deviates.

https://en.wikipedia.org/wiki/Normal_distribution

https://en.wikipedia.org/wiki/Box–Muller_transform

So if u1 and u2 are uniform deviartes on (0,1), then we expect that

sqrt(-2*log(u1))*cos(2*pi*u2)

is a normally distributed random variable, as is also

sqrt(-2*log(u1))*sin(2*pi*u2)

We can test that...

N = 1000000;
u1 = rand(1,N);
u2 = rand(1,N);

r = sqrt(-2*log(u1)).*cos(2*pi*u2);
hist(r,1000,'norm','pdf')

It looks reasonable to me.

Now we can try your RNG.

MinSec = fix(clock); 
seed = 100*(MinSec(5) + MinSec(6))

seed = 5600

U1 = myrandi(seed,N);
U2 = myrandi(seed,N);

r = sqrt(-2*log(U1)).*cos(2*pi*U2);
hist(r,1000,'norm','pdf')

And it looks like warmed up crapola on a shingle. ;-)

So what did you do? I see you used the same starting seed for BOTH CASES, so both U1 and U2!!!! This means U1 and U2 will be the same sequence of numbers. For example,

U1 = myrandi(seed,10)

U1 = 1×10

0 0.1709 0.3418 0.5127 0.6836 0.8545 0.0254 0.1963 0.3673 0.5382

U2 = myrandi(seed,10)

U2 = 1×10

0 0.1709 0.3418 0.5127 0.6836 0.8545 0.0254 0.1963 0.3673 0.5382

So what you wrote has U1 and U2 as completely correlated sequences. Even if they really were even poorly generated pseudo-random numbers, they don't qualify for the Box-Muller transformation to Normal.

At the very least, you need to use DIFFERENT starting seeds for each sequence. So now let me try it with at least ifferent seeds. I'm too lazy to use anything intelligent to come up with two distinct seeds, so I'll just try one sequence, then split it into two sequences.

U1 = myrandi(seed,2*N);
U2 = U1(1:2:end);
U1 = U1(2:2:end);

r = sqrt(-2*log(U1)).*cos(2*pi*U2);
hist(r,1000,'norm','pdf')

And yet worse crapola, but do you see the pattern developing? This may be getting into the issue of the poorness of that RNG, as not truly very random. So, next, let me look at the RNG code you wrote. I'll do this in another comment.

function Urand = myrandi(seed,sequence_length)
la=65539;
N=sequence_length;
P=2^31;
%change N accordingly 
%part a 
Urand=zeros(1,N);
seed(1)=seed;
for i=2:N 
    Urand(i)= mod((la*seed*(i-1)),P);
    
end
Urand = Urand./(P-1);
end

Answer 2

[John D'Errico]

After looking at your code in my other answer, the problem clearly reduces to your RNG. Is that even a valid linear congruential RNG? GOD NO! For example, I'll use a simple starting seed to show what is happening.

seed = 1000;

myrandi(seed,10)

ans = 1×10

0 0.0305 0.0610 0.0916 0.1221 0.1526 0.1831 0.2136 0.2442 0.2747

Do you see that each such random number is just a multiple of the first one, and that it starts out at 0? We could have gotten the same result from linspace. You had this:

Urand(i)= mod((la*seed*(i-1)),P);

Effectively, you have not much more than a call to linspace, at least for the first few thousand random numbers you generated. And even after that, it is not much different.

Was that really your plan? Were you directed to write the RNG like that?

https://en.wikipedia.org/wiki/Linear_congruential_generator

Try out this generator.

myslightlybetterrand(seed,10)

ans = 1×10

0.0305 0.1831 0.8240 0.2959 0.3597 0.4950 0.7322 0.9386 0.0416 0.8024

Do you see these already look now more random?

histogram(myslightlybetterrand(seed,100000),1000,'norm','pdf')

And again, this is starting to look more random. Still a completely crappy RNG, as I could surely show with some effort.

N = 1000000;
seed1 = 56554;
seed2 = 343453;
u1 = myslightlybetterrand(seed1,N);
u2 = myslightlybetterrand(seed2,N);
r = sqrt(-2*log(u1)).*cos(2*pi*u2);
hist(r,1000,'norm','pdf')

And that works. Your problem was in the basic RNG.

function Seq = myslightlybetterrand(seed,N)
  Seq = zeros(1,N);
  % don't use the seed as the first element of your sequence.
  % make it effectively the 0'th element instead.
  Seq0 = seed;
  la = 65539;
  P = 2^31;
  
  for n = 1:N
    Seq0 = mod(Seq0*la,P);
    Seq(n) = Seq0;
  end
  
  Seq = Seq/(P-1);
end

function Urand = myrandi(seed,sequence_length)
la=65539;
N=sequence_length;
P=2^31;
%change N accordingly 
%part a 
Urand=zeros(1,N);
seed(1)=seed;
for i=2:N 
    Urand(i)= mod((la*seed*(i-1)),P);
    
end
Urand = Urand./(P-1);
end

[Faisal Al-Wazir]

Were you directed to write the RNG like that?

the answer is yes there are a lot of things that the dr instructs us to not change in this code.

thank you so much john I learned a lot from you today I hope you have a wonderful day

[John D'Errico]

Ok then. I was wondering just a wee bit about that. What you cannot control... :)