关于求解马尔科夫链的平稳分布，大牛进！！！

yangfan715

知道马尔科夫链的转移流量图，如何求解马尔科夫链的平稳分布？或者说，知道N*N矩阵中各元素之间的关系，如何求各元素的值？在线等。。。。。。
例如：状态（i，j），0<i,j<N，转移流程图如下：

如何求状态（i，j）的平稳分布？

zsn828421

happy

一个马尔科夫链的matlab程序

1. function [chain,state]=markov(T,n,s0,V);
   
2. %function [chain,state]=markov(T,n,s0,V);
   
3. %  chain generates a simulation from a Markov chain of dimension
   
4. %  the size of T
   
5. %
   
6. %  T is transition matrix
   
7. %  n is number of periods to simulate
   
8. %  s0 is initial state (initial probabilities)
   
9. %  V is the quantity corresponding to each state
   
10. %  state is a matrix recording the number of the realized state at time t
    
11. %
    
12. % Original author: Tom Sargent
    
13. % Comments added by Qiang Chen
    
14. 
    
15. 
    
16. [r c]=size(T);   % r is # of rows, c is # of columns of T
    
17. if nargin == 1;  % "nargin" refers to "number of arguments in". So only T is provided in this case
    
18.   V=[1:r];
    
19.   s0=1;
    
20.   n=100;
    
21. end;
    
22. if nargin == 2;  % both T and n are provided
    
23.   V=[1:r];
    
24.   s0=1;
    
25. end;
    
26. if nargin == 3;  % T, n and S0 are provided
    
27.   V=[1:r];
    
28. end;
    
29. % check if the transition matrix T is square
    
30. if r ~= c;  
    
31.   disp('error using markov function');
    
32.   disp('transition matrix must be square');
    
33.   return;  % break the program and return
    
34. end;
    
35. % check if each row of T sums up to 1
    
36. for k=1:r;
    
37.   if sum(T(k,:)) ~= 1;
    
38.    disp('error using markov function')
    
39.     disp(['row ',num2str(k),' does not sum to one']);  % "num2str" converts numbers to a string.
    
40.     disp(' it sums to :');
    
41.     disp([ sum(T(k,:)) ]);
    
42.     disp(['normalizing row ',num2str(k),'']);
    
43.     T(k,:)=T(k,:)/sum(T(k,:));
    
44.   end;
    
45. end;
    
46. [v1 v2]=size(V);
    
47. if v1 ~= 1 | v2 ~=r   % "|" means "or"
    
48.   disp('error using markov function');
    
49.   disp(['state value vector V must be 1 x ',num2str(r),''])
    
50.   if v2 == 1 &v2 == r;
    
51.     disp('transposing state valuation vector');
    
52.     V=V';  % change it to a column vector
    
53.   else;
    
54.     return;
    
55.   end;  
    
56. end
    
57. if s0 < 1 |s0 > r;
    
58.   disp(['initial state ',num2str(s0),' is out of range']);
    
59.   disp(['initial state defaulting to 1']);
    
60.   s0=1;
    
61. end;
    
62. 
    
63. % The simulation of Markov chain formally starts from here
    
64. %T
    
65. %rand('uniform');
    
66. 
    
67. X=rand(n-1,1);  % generate (n-1) random numbers drawn from uniform distribution on [0,1], each number to be used in one simulation.
    
68. 
    
69. s=zeros(r,1);   % initiate the state vector "s" to be a rx1 zero vector
    
70. 
    
71. s(s0)=1;        % change the "s0"th element of "s" to 1
    
72. 
    
73. cum=T*triu(ones(size(T)));
    
74. % "triu(ones(size(T)))" generates an upper triangular matrix with all elements equal to 1
    
75. % cum is a rxr matrix whose ith column is the cumulative sum from the 1st column to the ith column
    
76. % the ith row of cum is the cumulative distribution for the next period given the current state.
    
77. 
    
78. for k=1:length(X);  % "length(X)" returns the size of the longest dimension of X. "k" indicates the kth simulation.
    
79.    
    
80.   state(:,k)=s;     % state is a matrix recording the number of the realized state at time k
    
81.   
    
82.   ppi=[0 s'*cum];   % this is the conditional cumulative distribution for the next period given the current state s
    
83.   
    
84.   s=((X(k)<=ppi(2:r+1)).*(X(k)>ppi(1:r)))';   
    
85.   % compares each element of ppi(2:r+1) or ppi(1:r) with a scalar X(k), and
    
86.   % returns 1 if the inequality holds and 0 otherwise
    
87.   % this formula assigns 1 when both inequalities hold, and 0 otherwise
    
88.   
    
89. end;
    
90. chain=V*state;

复制代码

hanwangwan

我想问一下这个转移流量图是怎么建立的啊?