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

3 views (last 30 days)

Show older comments

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)];

##### 1 Comment

per isakson
on 15 Jul 2012

It is not possible to just copy&paste and run your code. A few values are missing.

### Accepted Answer

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)))

### More Answers (2)

Jan
on 15 Jul 2012

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
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?

### See Also

### Categories

### Products

### Community Treasure Hunt

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

Start Hunting!