Fix the seed in Matlab R2011b | Why RandStream failed?

I employed the following codes to fix the seed in Matlab R2011b

s = RandStream('mt19937ar','Seed',0);

RandStream.setGlobalStream(s);

Then I tried to draw some uniform/normal/wishart random numbers in thousands of iterations, but the random numbers are the same only for first iteration when I restart Matlab. The iterations are like:

for i = 1: n

    a = rand(1,1);
    b = randn(3,1);
    c = wishrnd(omega,t);

end

Though it doesn't work within my complicated codes of real interest, I experimented with a simpler version as the above and it works whenever I restart Matlab.

I am really confused what is the source of this failure. Any insights are appreciated. Many thanks.

Allen

Answer 1

[the cyclist]

The rng() function is available in R2011b, and is the recommended method for seeding the random number generators.

doc rng

for details.

[Allen Liu]

Actually, I've tried 'rng' as well, but the problem remains. It seems that somewhere of the code changed the stream in iterations, however, I'm pretty sure that I never define the stream again except that in the very beginning.

Do you have any idea that any operation may cause this kind of problem? Thanks.

Answer 2

[Peter Perkins]

Allen, I can't tell what you mean by, "the random numbers are the same only for first iteration when I restart Matlab." You may mean that you expect to be able to restart MATLAB, run the same code, and see exactly the same numbers in all iterations of the loop, and what you actually see is that they are the same only for the first iteration. If you are really executing those first two lines that set the global random number stream, then you should get repeatable results. In fact, since those two lines do exactly what MATLAB already does at startup, you should see repeatable results even if you don't execute them.

It would be good if you would provide more specific information.

If the above is what you mean, then there must be something going on in your code that you have not mentioned. Especially since you can't reproduce the prblem in a simple example. My suggestion would be to set a breakpoint somewhere in your code at the point when you think things have gone wrong. When execution stops at that breakpoint, check and see what the global stream is, because it's possible that code that you are calling is changing it.

Beyond that, you're going to have to be more specific. Hope this helps.

[Allen Liu]

Yes, you correctly understand my purpose. And actually, it works in my simple toy codes though failed in my code of interest.

I'll try to call the global stream to see whether that changed (of course, it has changed for some reason). But I am sure that after establishing the global setting at the very beginning, I did not use any command to alter the stream (but it happened). Do you have any idea which kind of operation can result in this problem?

Maybe I'll post a simplified codes later. Thanks.

Answer 3

[Allen Liu]

To follow up, I've tried my best to debug it, however, I'm still not able to sort out the source of varying random numbers. Actually I found the seed had never been altered! Weird (>.<)~~

I post the simplified codes below which can't be abridged further in order to capture the fundamental features of my codes. You may find it still a bit messy, sorry for that. Codes are based on Kalman filtering and Gibbs sampling, but it seems not essential to understand why random numbers are varying which are key to draw kdraw and Bdraw. The desired result must present the same kpmean and temp1_mat whenever you run the codes. Random numbers are generated ONLY in functions "ts_prior", "wish", "gck" and "carter_kohn1", and please just believe that codes are error free.

I can't thank you enough for your great help and invaluable time.

clear all;
clc;
%========= fix seed ===============================
%in Matlab R2011b
s = RandStream('mt19937ar','Seed',0);
RandStream.setGlobalStream(s);

% rng(1111);

% stream = RandStream.getGlobalStream;
% reset(stream,0);

%-----------------------------LOAD DATA------------------------------------
Y = xlsread('experiment.xlsx');  

% Number of observations and dimension of X and Y
t = size(Y,1);
M = size(Y,2);

% training sample | initialize Kalman filter
tau = 40;  % size of training sample, CAN BE ALTERED according to the size of data
p = M; 
plag = 2; 
numa = p*(p-1)/2;
ylag = mlag2(Y,plag); % Y is [T x m]. ylag is [T x (nk)]
ylag = ylag(plag+1:t,:); 
m = p + plag*(p^2); % dimension of parameters

% Create Z_t matrix
Z = zeros((t-tau-plag)*p,m);
for i = tau+1:t-plag
    ztemp = eye(p); 
    for j = 1:plag        
        xtemp = ylag(i,(j-1)*p+1:j*p);
        xtemp = kron(eye(p),xtemp);
        ztemp = [ztemp xtemp];  %#ok<AGROW>
    end
    Z((i-tau-1)*p+1:(i-tau)*p,:) = ztemp;
end

y = Y(tau+plag+1:t,:)'; % sample used for estimation
t=size(y,2);    % size of the sample used for estimation 

%----------------------------PRELIMINARIES---------------------------------
% Set some Gibbs - related preliminaries
nrep = 50;  % Number of replications
nburn = 0;   % Number of burn-in-draws
apart = 1;   % save every apart-th draw 
it_print = 10;  %Print in the screen every "it_print"-th iteration

%========= PRIORS:====================================
%========= PRIORS ON TRANSISION PROBS (Beta) 
a_prob = .5;
b_prob = .5;
ap_0 = a_prob*ones(3,1);    
bp_0 = b_prob*ones(3,1);    

ap=zeros(3,1);
bp=zeros(3,1); %parameters for Beta distribution: alpha, beta

%=========PRIORS ON TIME-VARYING PARAMETERS AND THEIR COVARIANCES
[B_OLS,VB_OLS,A_OLS,sigma_OLS,VA_OLS]= ts_prior(Y,ylag(1:tau,:),tau,p,plag);

% Set some hyperparameters here (see page 831, end of section 4.1)
k_Q = 0.01;

%-------- Now set prior means and variances (_prmean / _prvar)
% B_0 ~ N(B_OLS, 4Var(B_OLS))
B_0_prmean = B_OLS;
B_0_prvar = 4*VB_OLS;

% Note that for IW distribution: scale and shape parameters...
% Q ~ IW(k2_Q*size(subsample)*Var(B_OLS),size(subsample))
Q_prmean = ((k_Q).^2)*tau*VB_OLS;
Q_prvar = tau;
%========= INITIALIZE MATRICES=======================
% Specify covariance matrices for measurement and state equations
consQ = 0.0001;
consH = 0.01;
Qdraw = consQ*eye(m);
Qchol = sqrt(consQ)*eye(m);
Ht = kron(ones(t,1),consH*eye(p));
Htsd = kron(ones(t,1),sqrt(consH)*eye(p));
Bdraw = zeros(m,t);
Zs = kron(ones(t,1),eye(p));
kdraw = 1*ones(t,3);
pdraw = .5*ones(1,3);
kold = kdraw;
kmean = zeros(t,3);
kvals = ones(2,1);   % k is a binary number 0 or 1
kvals(1,1) = 0;
kprior = .5*ones(2,3);

% Storage matrices for posteriors 
B_post1 = zeros(5,t*nrep/apart); 
B_post2 = zeros(5,t*nrep/apart);
B_post3 = zeros(5,t*nrep/apart);
B_post4 = zeros(6,t*nrep/apart);
Q_post = zeros(m,m*nrep/apart);
k_post = zeros(t,3*nrep/apart);
kp_post = zeros(nrep/apart,3);
B_postmean = zeros(m,t);
Qmean = zeros(m,m);
kpmean = zeros(1,3);

%==========store generated random numbers ========================
temp1_mat = zeros(t-1,(nburn+nrep)/apart); 
%----------------------------- END OF PRELIMINARIES ------------------------

%=========== START SAMPLING ==========================

tic; % This is just a timer
disp('Number of iterations');
for irep = 1:nrep + nburn    % GIBBS iterations starts here
   % Print iterations
    if mod(irep,it_print) == 0
       disp(irep);toc;
    end

% STEP I: Sample B from p(B|y,A,Sigma,V)         
  % I.1: draw K1 index and related probabilities
pdrawa = zeros(3,1);
ap(1,1) = ap_0(1,1) + sum(kdraw(:,1));
bp(1,1) = bp_0(1,1) + t - sum(kdraw(:,1));
pdrawa(1,1) = betarnd(ap(1,1),bp(1,1));
pdraw(1,1) = pdrawa(1,1);
kprior(2,1) = pdrawa(1,1);
kprior(1,1) = 1 - kprior(2,1);

[kdrawa,lpyK,tempv_mat1] = gck(y,zeros(p,t),Z,Htsd,zeros(m,t),...
  kron(ones(t,1),eye(m)),Qchol,kold(:,1),t,B_0_prmean,B_0_prvar,2,kprior(:,1),kvals,p,m);
kdraw(:,1) = kdrawa;   % update kdraw
kold(:,1) = kdraw(:,1);
temp1_mat(:,irep) = tempv_mat1;  %store rand(1,1)'s for drawing k | target  1   needs to be checked

% draw parameters using Kalman Filter
[Bdrawc,log_lik1] = carter_kohn1(y,Z,Ht,Qdraw,m,p,t,B_0_prmean,B_0_prvar,kdraw(:,1));
Bdraw = Bdrawc; % update parameters: Bdraw

Btemp = Bdraw(:,2:t)' - Bdraw(:,1:t-1)';
sse_2 = zeros(m,m);
for i = 1:t-1
    sse_2 = sse_2 + Btemp(i,:)'*Btemp(i,:);
end

Qinv = inv(sse_2 + Q_prmean);
Qinvdraw = wish(Qinv,t+Q_prvar);
Qdraw = inv(Qinvdraw);
Qchol = chol(Qdraw)';

    %--------SAVE AFTER-BURN-IN DRAWS AND IMPULSE RESPONSES --------------
  if irep > nburn
   if mod(irep-nburn, apart) == 0
   B_post1(:,t*((irep-nburn)/apart-1)+1:t*(irep-nburn)/apart) = Bdraw(1:5,:);
   B_post2(:,t*((irep-nburn)/apart-1)+1:t*(irep-nburn)/apart) = Bdraw(6:10,:);
   B_post3(:,t*((irep-nburn)/apart-1)+1:t*(irep-nburn)/apart) = Bdraw(11:15,:);
   B_post4(:,t*((irep-nburn)/apart-1)+1:t*(irep-nburn)/apart) = Bdraw(16:21,:);
            Q_post(:,m*((irep-nburn)/apart-1)+1:m*(irep-nburn)/apart) = Qdraw;
            k_post(:,p*((irep-nburn)/apart - 1)+1:p*(irep-nburn)/apart) = kdraw;
            kp_post((irep-nburn)/apart,:) = pdraw;
            B_postmean = B_postmean + Bdraw;
            Qmean = Qmean + Qdraw;
            kmean = kmean + kdraw;
            kpmean = kpmean + pdraw;
   end
   end
  end

%===========GIBBS SAMPLER ENDS HERE============================

B_postmean = B_postmean./(nrep/apart);
Qmean = Qmean./(nrep/apart);
kmean = kmean./(nrep/apart);    
kpmean = kpmean./(nrep/apart);  % target 2 needs to be checked

%%-------- funtions are as follows -----------------------------------
function [Xlag] = mlag2(X,p)
[Traw,N]=size(X);
Xlag=zeros(Traw,N*p);
for ii=1:p
    Xlag(p+1:Traw,(N*(ii-1)+1):N*ii)=X(p+1-ii:Traw-ii,1:N);
end

%-----------------------------------------------------------------
function [bols,bvar,a0,ssig1,a02mo] = ts_prior(rawdat,x,tau,p,plag)

yt = rawdat(plag+1:tau+plag,:)';  
m = p + plag*(p^2);   
Zt = zeros(tau*p,m);

for i = 1:tau
    ztemp = eye(p); 
    for j = 1:plag        
        xtemp = x(i,(j-1)*p+1:j*p);
        xtemp = kron(eye(p),xtemp);
        ztemp = [ztemp xtemp];  %#ok<AGROW>
    end
    Zt((i-1)*p+1:i*p,:) = ztemp;
end

yy = reshape(rawdat(plag+1:tau+plag,:)',p*tau,1);
bols = (Zt'*Zt) \ Zt' * yy;   %B_OLS
sse2 = zeros(p,p);

y_hat = [];
for i = 1:tau
    zhat1 = Zt((i-1)*p+1:i*p,:);
    residual = yt(:,i) - zhat1*bols;  
    y_hat = [y_hat residual]; %#ok<AGROW>   
    sse2 = sse2 + residual*residual';
end
hbar = sse2./tau; 

Q = zeros(p*plag+1, p*plag+1);    
for i = 1:tau
    temp = [1 x(i,1:p) x(i,p+1:plag*p)]; 
    Q = Q + temp' * temp;
end
Q = (1/tau)*Q;      
bvar = kron(hbar, inv(Q)); 

achol = chol(hbar)';    
ssig = zeros(p,p); 
for i = 1:p
  ssig(i,i) = achol(i,i);    
end

ssig1 = zeros(p,1);
for i = 1:p
   ssig1(i,1) = log(ssig(i,i)^2);  % log (sigma_t^2)
end

achol = chol(hbar)';
ssig = diag(diag(achol));
achol = inv(achol/ssig);

numa = p*(p-1)/2;
a0 = zeros(numa,1);
ic = 1;
for ii = 2:p
  a0(ic:ic+ii-2,1) = achol(ii,1:ii-1);
  ic = ic + ii - 1;
end

% Simulate var-covar matrix of A
% Log sigma gives the state h(.) functions
ssig1 = 2*log(diag(ssig) );

hbar1 = inv(tau*hbar);  % inv(Scatter Matrix)
a02mo = zeros(numa,numa);
a0mean = zeros(numa,1);
for irep = 1:4000
  hdraw = wish(hbar1,tau);  % Wishart distribution
  hdraw = inv(hdraw); % Inverse wishart
  achol = chol(hdraw)';
  ssig = diag(diag(achol));
  achol = inv(achol/ssig);
  a0draw = zeros(numa,1);
  ic = 1;
  for ii = 2:p
    a0draw(ic:ic+ii-2,1) = achol(ii,1:ii-1);
    ic = ic + ii - 1;
  end
  a02mo = a02mo + a0draw*a0draw';
  a0mean = a0mean + a0draw; 
end

a02mo = a02mo./4000;
a0mean = a0mean./4000;
a02mo = a02mo - a0mean*a0mean';  % Var(X) = E(XX') - (EX)(EX)'

%--------------------------------------------------------------------
function A = wish(h,n)

A = chol(h)'*randn(size(h,1),n);
A = A*A';

%-----------------------------------------------------------------
function [kdraw,lpy2n,tempv_mat1] = gck (yg,gg,hh,capg,f,capf,sigv,kold,t,ex0,vx0,nvalk,kprior,kvals,p,kstate)

%==================================================================    
 tempv_mat1 = [];    % store the random number generated from rand(1,1)
%=================================================================---

% GCK's Step 1 on page 821
lpy2n=0;
mu = zeros(t*kstate,1);
omega = zeros(t*kstate,kstate);

for i = t-1:-1:1
  gatplus1 = kold(i+1,1)*sigv; 

    ftplus1 = capf(kstate*i+1:kstate*(i+1),:);
    cgtplus1 = capg(i*p+1:(i+1)*p,:);
    htplus1 = hh(i*p+1:(i+1)*p,:)';  
    rtplus1 = (htplus1'*gatplus1)*(htplus1'*gatplus1)' + cgtplus1*cgtplus1';
    btplus1 = gatplus1*gatplus1'*htplus1 / rtplus1; 
    atplus1 = (eye(kstate) - btplus1*htplus1')*ftplus1;

    cct1 = eye(kstate) - gatplus1'*(htplus1 / rtplus1)*htplus1'*gatplus1; 
    cct2= chol(cct1)';
    ctplus1 = gatplus1*cct2;

    otplus1 = omega(kstate*i+1:kstate*(i+1),:);
    dtplus1 = ctplus1'*otplus1*ctplus1 + eye(kstate);
    omega(kstate*(i-1)+1:kstate*i,:) = atplus1'*(otplus1 - otplus1*(ctplus1/dtplus1)*ctplus1'*otplus1)*atplus1 + ftplus1'*(htplus1/rtplus1)*htplus1'*ftplus1;
    satplus1 = (eye(kstate) - btplus1*(htplus1'))*f(:,i+1) - btplus1*gg(:,i+1);
    mutplus1 = mu(kstate*i+1:kstate*(i+1),:);
    mu(kstate*(i-1)+1:kstate*i,:) = atplus1'*(eye(kstate) - otplus1*(ctplus1/dtplus1)*ctplus1')*(mutplus1 - otplus1*(satplus1 + btplus1*yg(:,i+1))) + ftplus1'*(htplus1/rtplus1)*(yg(:,i+1) - gg(:,i+1) - htplus1'*f(:,i+1));  
  end

% GCKs Step 2 on pages 821-822
kdraw = kold;
ht = hh(1:p,:)';
ft = capf(1:kstate,:);
gat = zeros(kstate,kstate);
rt = ht'*ft*vx0*ft'*ht + ht'*(gat*gat')*ht+ capg(1:p,:)*capg(1:p,:)';
jt = (ft*vx0*ft'*ht + gat*gat'*ht) / rt;
mtm1 = (eye(kstate) - jt*ht')*(f(:,1) + ft*ex0) + jt*(yg(:,1) - gg(:,1));
vtm1 = ft*vx0*ft'+ gat*gat' - jt*rt*jt';

%store probabilities of k = 1
lprob = zeros(nvalk,1); 

for i = 2:t

    ht = hh((i-1)*p+1:i*p,:)';
    ft = capf(kstate*(i-1)+1:kstate*i,:);
    for j = 1:nvalk
        gat = kvals(j,1)*sigv;

        rt = ht'*ft*vtm1*ft'*ht + ht'*(gat*gat')*ht + capg((i-1)*p+1:i*p,:)*capg((i-1)*p+1:i*p,:)';

        jt = (ft*vtm1*ft'*ht + gat*gat'*ht)/ rt;
        mt = (eye(kstate) - jt*ht')*(f(:,i) + ft*mtm1) + jt*(yg(:,i) - gg(:,i));
        vt = ft*vtm1*ft' + gat*gat' - jt*rt*jt';        
        lpyt = -.5*log(det(rt)) - .5*(yg(:,i) - gg(:,i) - ht'*(f(:,i) + ft*mtm1))'*(rt\(yg(:,i) - gg(:,i) - ht'*(f(:,i) + ft*mtm1)));

        [V, lambda] = eig(vt);
        tt = V * sqrt(lambda); 

        ot = omega(kstate*(i-1)+1:kstate*i,:);
        mut = mu(kstate*(i-1)+1:kstate*i,:);
        tempv = eye(kstate) + tt'*ot*tt;
        lpyt1n = -.5*log(det(tempv)) -.5*(mt'*ot*mt - 2*mut'*mt - (mut - ot*mt)'*(tt/tempv)*tt'*(mut - ot*mt));

        lprob(j,1) = log(kprior(j,1)) + lpyt1n + lpyt;

        if i == 2;
            lpy2n = lpyt1n + lpyt;
        end
    end
    pprob = exp(lprob)./sum(exp(lprob));

    %========================================================--
    % generate random uniformly distributed number
    tempv = rand(1,1);
    tempv_mat1 = [tempv_mat1; tempv];
    %========================================================--

    tempu = 0;

    for j = 1:nvalk
        tempu = tempu + pprob(j,1);

        if tempu > tempv;
            kdraw(i,:) = kvals(j,:);
            break
        end

    end

    gat = kdraw(i,1)*sigv;  
    rt = ht'*ft*vtm1*ft'*ht + ht'*(gat*gat')*ht + capg((i-1)*p+1:i*p,:)*capg((i-1)*p+1:i*p,:)';
    jt = (ft*vtm1*ft'*ht + gat*gat'*ht) / rt;
    mtm1 = (eye(kstate) - jt*ht')*(f(:,i) + ft*mtm1) + jt*(yg(:,i) - gg(:,i));
    vtm1 = ft*vtm1*ft' + gat*gat' - jt*rt*jt';           
  end
  end

%---------------------------------------------------------------------------
function [bdraw,log_lik] = carter_kohn1(y,Z,Ht,Qt,m,p,t,B0,V0,kdraw)

% Kalman Filter
bp = B0;
Vp = V0;
bt = zeros(t,m);
Vt = zeros(m^2,t);
log_lik_ = 0;
for i=1:t
  R = Ht((i-1)*p+1:i*p,:);
  H = Z((i-1)*p+1:i*p,:);
  cfe = y(:,i) - H*bp;   
  f = H*Vp*H' + R;   
  inv_f = H'/f;
  log_lik_ = log_lik_ + log(det(f)) + (cfe'/ f)*cfe;
  btt = bp + Vp*inv_f*cfe;
  Vtt = Vp - Vp*inv_f*H*Vp;
  if i < t
      bp = btt;
      Vp = Vtt + kdraw(i,:)*Qt;
  end
  bt(i,:) = btt';
  Vt(:,i) = reshape(Vtt,m^2,1);
end

log_lik = - (t/2) * log(2*pi) - (1/2)*log_lik_; 

bdraw = zeros(t,m);

%=== generate multivariate normal random vector ========
bdraw(t,:) = btt' + randn(1,m)*chol(Vtt);
%=======================================================

for i=1:t-1
  bf = bdraw(t-i+1,:)';
  btt = bt(t-i,:)';
  Vtt = reshape(Vt(:,t-i),m,m);
  f = Vtt + kdraw(t-i,:)*Qt;
  inv_f = Vtt/f;
  cfe = bf - btt;
  bmean = btt + inv_f*cfe;
  bvar = Vtt - inv_f*Vtt;

    %=== generate multivariate normal random vector ========
    bdraw(t-i,:) = bmean' + randn(1,m)*chol(bvar)   
    %=======================================================
  end
  bdraw = bdraw';

%------------------------------------------------------------------------

[Honglei Chen]

Hi Allen, as I explained in the other forum where you posted the same question, it is unlikely that people here will have time to read your script especially it doesn't wrong without your Excel file. You need to provide simple reproduction steps. Peter has given excellent advice above regarding how to debug your script. I suspect that certain function in your script calls some other functions that alter the random number generator.

Answer 4

[Honglei Chen]

Do you have to set it to global stream? If not, why not generate those random numbers in form of

rand(s,1,n)

so you by pass the possibility that other things mess up the random number state behind the scene?

[Allen Liu]

I thought one can simply use any excel file instead, sorry for that.

Please forgive me for the next stupid question. What do you mean by

rand(s,1,n)

cause I never saw this before, does that mean one set the seed = s first, then generate an 1-by-n vector? or something else?

Actually, what I need is just a bunch of random numbers - uniformly / univariate or multivariate normally distributed, but it's better to fix these random numbers for convincing the reviewers, although essentially my results would not be overturned.

Thanks heaps.

[Honglei Chen]

Sorry for the confusion, s stands for stream. Something like this

rs = RandStream.create('mt19937ar');
rand(rs,1,10)

[Allen Liu]

Cool, I'll try it when I am in office tomorrow.

Thanks a bunch, dude.