Do I meet the limitation of Matlab? A unwanted pattern keeps showing from a random simulation.

3 views (last 30 days)
Hsinho
Hsinho on 15 Jul 2012
I always get a unwanted pattern from a random simulation no matter how I change the way of random number generation. In the beginning, I thought random number generation causes the problem. After I tried all possible ways to generate the random numbers, the unwanted pattern remains.
Do I meet the limitation of Matlab? How can I overcome it?
The simulation result of E variables in the concatenated E matrix
4.112 9.494 11188.4671 9.2810 1.8651
5.058 9.584 11188.7245 9.2860 1.9200
The values of the 3rd column of E matrix are always higher than others. The values of the first and last columns are always lower than others. It looks like a kind of symmetric distribution from column 1 to 5. Only the first column of E is the expected result. The values in the last column are too low and not expected. This kind of pattern also can be observed when there are more columns (e.g. R=10 or 20).
What I want to do is to repeat what I got as in the first column of E matrix for any R variable I assign. In my code, each column means a independent simulation. In the 2.1e7 for-loop, I want to simultaneously do independent simulations as many as the R assigned.
However, when R=1, it is fine, but when R>1, I got the problem as described above. The calculation of U and E don't cause the problem since they are done by matrix operation.
I hope enthusiastic Matlab masters/experts/elites/lovers can help me to perceive the zen of Matlab to accomplish what I want to do.
The following lines are my simplified code. There are two versions to show how I generate the random numbers for every single either rand or randi call for the random stream.
k=a constant;
T=a constant;
len=130;
R=5;
y=zeros(len,R);
N=zeros(1,R);
% version 1: Use the default global stream
for n=1:2.1e7;
yold=y;
N=(0:len:(R-1)*len)+randi(len,1,R);
y(N)=random('unif',-1,8,1,R);
ynew=y;
U=myfun(y) % size(U) is a len-by-R matrix
Enew=sum(U); % size(Enew) is a 1-by-R matrix
try
E=cat(1,Eold,Enew);
catch
Eold=Enew;
end
try
negative=find(E(2,:)-E(1,:)<0);
positive=find(E(2,:)-E(1,:)>0);
if isempty(negative)~=1
Eold(:,negative)=Enew(:,negative);
yold(:,negative)=ynew(:,negative);
end
if isempty(positive)~=1
ind=find((exp(-(E(2,positive)-E(1,positive))/(k*T)))> rand(1,length(positive));
if isempty(ind)~=1
Eold(:,ind)=Enew(:,ind);
yold(:,ind)=ynew(:,ind);
end
end
y=yold;
catch
end
end
% version 2: The random streams were specified either by [s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15]=RandStream.create('mrg32k3a','NumStreams',15);
% or by the seeds selected by randi(2^32,1,15)
s1=RandStream('mt19937ar','Seed',3499200000);
s2=RandStream('mt19937ar','Seed',3890300000);
s3=RandStream('mt19937ar','Seed',545400000);
s4=RandStream('mt19937ar','Seed',3922900000);
s5=RandStream('mt19937ar','Seed',2716000000);
s6=RandStream('mt19937ar','Seed',418900000);
s7=RandStream('mt19937ar','Seed',1196100000);
s8=RandStream('mt19937ar','Seed',2348800000);
s9=RandStream('mt19937ar','Seed',4112500000);
s10=RandStream('mt19937ar','Seed',4144200000);
s11=RandStream('mt19937ar','Seed',676900000);
s12=RandStream('mt19937ar','Seed',4168700000);
s13=RandStream('mt19937ar','Seed',4111000000);
s14=RandStream('mt19937ar','Seed',2084700000);
s15=RandStream('mt19937ar','Seed',3437200000);
% either rand or randi calls for a specific stream assigned by a either way in version 2 above.
for n=1:2.1e7;
yold=y;
N1=(0:len:(R-1)*len);
N2=[randi(s1,len),randi(s2,len),randi(s3,len),randi(s4,len),randi(s5,len)];
N=N1+N2;y(N)=[-1+9*rand(s6),-1+9*rand(s7),-1+9*rand(s8),-1+9*rand(s9),-1+9*rand(s10)];
ynew=y;
U=myfun(y) % size(U) is a len-by-R matrix
Enew=sum(U); % size(Enew) is a 1-by-R matrix
try
E=cat(1,Eold,Enew);
catch
Eold=Enew;
end
try
negative=find(E(2,:)-E(1,:)<0);
positive=find(E(2,:)-E(1,:)>0);
if isempty(negative)~=1
Eold(:,negative)=Enew(:,negative);
yold(:,negative)=ynew(:,negative);
end
if isempty(positive)~=1
randE=[rand(s11),rand(s12),rand(s13),rand(s14),rand(s15)];
ind=find((exp(-(E(2,positive)-E(1,positive))/(k*T)))> randE(1,positive));
if isempty(ind)~=1
Eold(:,ind)=Enew(:,ind);
yold(:,ind)=ynew(:,ind);
end
end
y=yold;
catch
end
end
I even further change the assignment of s streams for rand or randi.
All the efforts still fail.
1. each column of y matrix and E matrix use the same stream
N2=[randi(s1,len),randi(s2,len),randi(s3,len),randi(s4,len),randi(s5,len)];
N=N1+N2; y(N)=[-1+9*rand(s1),-1+9*rand(s2),-1+9*rand(s3),-1+9*rand(s4),-1+9*rand(s5)];
randE=[rand(s1),rand(s2),rand(s3),rand(s4),rand(s5)];
2. each major random number generation step use the same stream
N2=randi(s1,len,1,R);
N=N1+N2;
y(N)=-1+9*rand(s2,1,R);
randE=rand(s3,1,R);
3. change s1 and s5 from both sides to the center (aiming for changing the pattern but still fail)
N2=[randi(s3,len),randi(s5,len),randi(s1,len),randi(s4,len),randi(s2,len)];
N=N1+N2; y(N)=[-1+9*rand(s3),-1+9*rand(s5),-1+9*rand(s1),-1+9*rand(s4),-1+9*rand(s2)];
randE=[rand(s3),rand(s5),rand(s1),rand(s4),rand(s2)];

Accepted Answer

Sean de Wolski
Sean de Wolski on 16 Jul 2012
Edited: Sean de Wolski on 16 Jul 2012
Just taking a quick glance:
Enew=sum(U);
This and other operations like it could certainly be the issue.
More
For example:
The following normal distribution is completely expected:
hist(sum(rand(1000)))
  6 Comments

Sign in to comment.

More Answers (2)

Jan
Jan on 15 Jul 2012
I'm convinced (this means I do not have a firm argument), that your problems are not caused by the random number generator, but by your algorithm. You have changed the RNG part efficiently, but still get the same unexpected behavior. This is a strong hint, that your program does not reflect your expectations.
  1 Comment
Hsinho
Hsinho on 16 Jul 2012
Thanks for your response. One mystery for me is that if the program goes wrong, why can I get the expected result when R=1? The code between Line 37 to 43 as I commented on the question above should be correct no matter what R is because all the columns of y_n=y(1:len-1,:); and y_n_plus_1=y(2:len,:) were selected. Hope we can solve this mystery.

Sign in to comment.


Paul Metcalf
Paul Metcalf on 16 Jul 2012
Edited: Walter Roberson on 16 Jul 2012
What's in myfun()?
You mention this: 'I want to simultaneously do independent simulations as many as the R assigned'.
Are you doing any parallel processing? If so, a poorly conceived parallel algorithm could be your problem... I am prepared to say I am 50% certain this is your problem.
PS - kudos for the formatting of your post!
  2 Comments
Paul Metcalf
Paul Metcalf on 18 Jul 2012
Sorry but I can't see any line numbers in your post.
Could you please clarify which code is line 37-43?

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!