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

Question

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

Answer 1

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