def stationary_distrib(F,residtol = 1.E-10,max_iter=100):
"""
Calculates the eigenvector of the matrix F with eigenvalue 1 (if it exists).
Parameters
----------
F : ndarray
A matrix known to have a single left eigenvector with
eigenvalue 1.
residtol : float or scalar
To improve the accuracy of the computation, the algorithm will
"polish" the final result using several iterations of the power
method, z^T F = z^T. Residtol gives the tolerance for the
associated relative residual to determine convergence.
maxiter : int
Maximum number of iterations to use the power method to reduce
the residual. In practice, should never be reached.
Returns
-------
z : ndarray
The eigenvector of the matrix F with eigenvalue 1. For a Markov
chain stationary distribution, this is the stationary distribution.
Normalization is chosen s.t. entries sum to one.
"""
L = len(F) # Number of States
M = np.eye(L)-F
q,r=qr(M)
z=q[:,-1] # Stationary dist. is last column of QR fact
z/=np.sum(z) # Normalize Trajectory
# Polish solution using power method.
for itr in xrange(max_iter):
znew = np.dot(z,F)
maxresid = np.max(np.abs(znew - z)/z) # Convergence Criterion
if maxresid < residtol:
break
else:
z = znew
return z/np.sum(z) # Return normalized (by convention)
评论列表
文章目录
