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

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

[per isakson]

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

Answer 1

[Sean de Wolski]

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

[Hsinho]

As I understand, size(U) is len-by-R and size(Enew) is 1-by-R because the default setting of sum is to compute the sum of each column.

Could you explain more clearly why you think this line is the problem? If so, may I have your suggestion?

[Sean de Wolski]

The sum of uniformally random numbers fits a normal distribution. Think of it this way: if we have 10 random numbers between 0 and 1, what is likely to be their mean? What is less likely to be their mean?

[Sean de Wolski]

I don't know what that calculation is doing. I can say that summing uniform random numbers makes the result non-random.

[Hsinho]

Thanks for your explanation about what I observed. From your example "hist(sum(rand(1000)))", I imagine this histogram is the normal distribution of a population of 1000-by-1000 uniformly random numbers.

Use the same thinking, in my present code, the summing of the 130-by-5 matrix showing the normal distribution of a population of 130-by-5 uniformly random numbers should be expected.

Now, the solution I consider is to FORCE each column is a population. But, How?

I still wonder whether it is possible to use the concept of vectorization to write the code.

One applicable way is to use R=1 and use "parfor" to run several R=1 simultaneously. I guess "parfor" will take longer time than vectorization if vectorization works.

Is it possible to force each column is a population with the concept of vectorization?

[Sean de Wolski]

I have no idea what you are trying to say. Don't vectorize. Get it working and then vectorize if necessary.

[Hsinho]

Thank you! I will do so.

Answer 2

[Jan]

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.

[Hsinho]

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.

Answer 3

[Paul Metcalf]

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!

[Hsinho]

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]

Sorry but I can't see any line numbers in your post.

Could you please clarify which code is line 37-43?