Calculates the mean hitting time of a state in a homogenous markov chain

meanHittingTime.m

function K = fixedPoint(M)
  [r c] = size(M);
  assert(r==c);
  I = eye(r);
  K = -rref(I-M)(:,r);
  K(r)=1;
  assert(norm(M*K-K)<.0001*norm(K));
end

function p = projectFromPRn(X)
  assert(X(1)!=0);
  X/=X(1);
  p=X(2:end);
end

function k = meanHittingTime(P,destination)
  [r c] = size(P);
  assert(r==c);
  N = r+1;
  M = [ones(N,1) [zeros(1,r); P]];
  M(:,destination+1)=zeros(1,N);
  k = projectFromPRn(fixedPoint(M));
end