Home > wafo > wsim > mcsim.m

mcsim

PURPOSE ^

Simulates a Markov chain.

SYNOPSIS ^

x=mcsim(P,T,x0)

DESCRIPTION ^

 MCSIM   Simulates a Markov chain.
 
  CALL: x = mcsim(P,T)
        x = mcsim(P,T,x0)
 
  x  = Simulated Markov chain.
 
  P  = Transition matrix.    [rxr]
  T  = Length of simulation.
  x0 = Initial state. (Default:  x0 will be from the stationary
       distribution of the Markov chain.) 
  
  Simulates a Markov chain with state space {1,2,...,r} 
  
  Example: 
    P=[0.8 0.2; 0.1 0.9]; 
    x=mcsim(P,100); plot(x) 
  
  See also  mc2stat, smcsim

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function x=mcsim(P,T,x0)
002 
003 %MCSIM   Simulates a Markov chain.
004 %
005 % CALL: x = mcsim(P,T)
006 %       x = mcsim(P,T,x0)
007 %
008 % x  = Simulated Markov chain.
009 %
010 % P  = Transition matrix.    [rxr]
011 % T  = Length of simulation.
012 % x0 = Initial state. (Default:  x0 will be from the stationary
013 %      distribution of the Markov chain.) 
014 % 
015 % Simulates a Markov chain with state space {1,2,...,r} 
016 % 
017 % Example: 
018 %   P=[0.8 0.2; 0.1 0.9]; 
019 %   x=mcsim(P,100); plot(x) 
020 % 
021 % See also  mc2stat, smcsim
022  
023 % Tested  on Matlab  5.3 
024 % 
025 % History: 
026 % Revised by PJ 26-Jul-2000 
027 %   updated help text. 
028 % Created by PJ (Pär Johannesson) 1997 
029 %   Copyright (c) 1997 by Pär Johannesson
030 %   Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997
031 
032 if nargin < 3
033   x0 = [];
034 end
035 
036 % Check that the rowsums of P are equal to 1
037 
038 sumP = sum(P');
039 if sum(sumP == 1) ~= length(P)
040   disp(['Warning: Rowsums of P not equal to 1. Renormalizing.']);
041   for i = 1:length(P)
042     P(i,:) = P(i,:)/sumP(i);
043   end
044 end
045 
046 x=ones(T,1);
047 e=rand(T,1);
048 cumsumP = cumsum(P')';
049 
050 % Initial state
051 
052 if isempty(x0)    % From stationary distribution
053 
054   ro=mc2stat(P);
055 
056   x(1) = min(find( e(1)<=cumsum(ro) ));
057 %  x(1) = sum( cumsum(ro) <= e(1) ) + 1;
058 
059 else
060 
061   x(1) = x0;   % Given
062 
063 end
064 
065 %---------------------
066 
067 for k=2:T
068 
069   x(k) = min(find( e(k)<=cumsumP(x(k-1),:) ));
070 
071 end
072  
073 %  x(k) = sum( cumsumP(x(k-1),:) <= e(k) ) + 1;
074

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003