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

1 次查看（过去 30 天）
Hsinho
Hsinho 2012-7-15
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)];

请先登录，再回答此问题。

采纳的回答

Sean de Wolski
Sean de Wolski 2012-7-16
编辑：Sean de Wolski 2012-7-16
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 个评论
显示 4更早的评论

请先登录，再进行评论。

更多回答（2 个）

Jan
Jan 2012-7-15
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 个评论
Hsinho
Hsinho 2012-7-16
编辑：Hsinho 2012-7-16
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.

请先登录，再进行评论。

Paul Metcalf
Paul Metcalf 2012-7-16
编辑：Walter Roberson 2012-7-16
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 个评论
Hsinho
Hsinho 2012-7-16
Thank you for answering. The myfun(y) was posted as above in the comment of the question between line 37 to 43.
I think I am doing a kind of parallel processing by vectorization but not by "parfor." The way of vectorization I did is as the line 37 to 43.
Please help me to identify the problem. Thank!
Paul Metcalf
Paul Metcalf 2012-7-18
编辑：Paul Metcalf 2012-7-18
Sorry but I can't see any line numbers in your post.
Could you please clarify which code is line 37-43?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Just for fun 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by