Markov Chains and partial views/Code

From GusWiki

M = [1,0,1,0,0,0,0,0,-1,0,-1,0,0,0,0,0,0; 0,1,0,1,0,0,0,0,0,-1,0,- 1,0,0,0,0,0;0,0,0,0,1,0,1,0,0,0,0,0,-1,0,-1,0,0;0,0,0,0,0,1,0,1,0,0,0,0,0,-1,0,- 1,0;1,1,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0;0,0,1,1,0,0,-1,- 1,0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,1,1,0,0,-1,-1,0,0,0;0,0,0,0,0,0,0,0,0,0,1,1,0,0,-1,- 1,0;1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1;0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,1;0,0,0,0,0,0,0,0,1,1, 1,1,0,0,0,0,1;0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1;1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0.3;0,1,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0.7;0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0.7;0,0,0,0,0,1,0,1,0,0,0,0,0,0, 0,0,0.3;1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5;0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0.5;0,0,0,0,0,0, 0,0,1,1,0,0,0,0,0,0,0.75;0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0.25] function ret = lastcolumn(mat) ret = mat(:,columns(mat)) endfunction function ret = matricize (vec) r(1,:) = vec(1:4) r(2,:) = vec(5:8) r(3,:) = vec(9:12) r(4,:) = vec(13:16) ret = r endfunction # given the parameters, joint matrix function ret = jointmatrix(M,parms) L = rref(M) J = L(1:12,:) #remove the zero rows J(13,:) = [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,parms(1)] J(14,:) = [0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,parms(2)] J(15,:) = [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,parms(3)] J(16,:) = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,parms(4)] I = rref(J) ret = matricize(lastcolumn(I)) endfunction function ret = entropy(vec) sum=0 for i=1:length(vec) sum+=vec(i)*log(1/vec(i)) endfor ret = sum endfunction function ret = normalize(vec) ret = vec/sum(vec) endfunction function ret = matentropy(mat) [a,b] = eig(mat) if (b(1,1)=1) distr = normalize(a(:,1)) ret = entropy(distr) endif endfunction