% Chutes and Ladders(TM), Markov Chain Problem
% Copyright (C) 2002 Jeffrey Humpherys.
%

SPIN = 6;       % The number of sides of the die.
N = 100;        % The number of squares.

%
% Fills the transition matrix as if there were no chutes or ladders.
%

A = zeros(N+1,N+1);
for i = 1:SPIN
    A = A + diag(ones(N+1-i,1),-i);
end;

%
% Now adjust for the chutes and ladders
%

B = A(38+1,:);A(38+1,:) = A(1+1,:) + B;A(1+1,:) = zeros(1,N+1);B=[];        % Ladder from 1->38
B = A(14+1,:);A(14+1,:) = A(4+1,:) + B;A(4+1,:) = zeros(1,N+1);B=[];        % Ladder from 4->14
B = A(31+1,:);A(31+1,:) = A(9+1,:) + B;A(9+1,:) = zeros(1,N+1);B=[];        % Ladder from 9->31
B = A(42+1,:);A(42+1,:) = A(21+1,:) + B;A(21+1,:) = zeros(1,N+1);B=[];      % Ladder from 21->42
B = A(84+1,:);A(84+1,:) = A(28+1,:) + B;A(28+1,:) = zeros(1,N+1);B=[];      % Ladder from 28->84
B = A(44+1,:);A(44+1,:) = A(36+1,:) + B;A(36+1,:) = zeros(1,N+1);B=[];      % Ladder from 36->44
B = A(67+1,:);A(67+1,:) = A(51+1,:) + B;A(51+1,:) = zeros(1,N+1);B=[];      % Ladder from 51->67
B = A(91+1,:);A(91+1,:) = A(71+1,:) + B;A(71+1,:) = zeros(1,N+1);B=[];      % Ladder from 71->91
B = A(100+1,:);A(100+1,:) = A(80+1,:) + B;A(80+1,:) = zeros(1,N+1);B=[];    % Ladder from 80->100

B = A(78+1,:);A(78+1,:) = A(98+1,:) + B;A(98+1,:) = zeros(1,N+1);B=[];      % Chute from 98->78
B = A(75+1,:);A(75+1,:) = A(95+1,:) + B;A(95+1,:) = zeros(1,N+1);B=[];      % Chute from 95->75
B = A(73+1,:);A(73+1,:) = A(93+1,:) + B;A(93+1,:) = zeros(1,N+1);B=[];      % Chute from 93->73
B = A(24+1,:);A(24+1,:) = A(87+1,:) + B;A(87+1,:) = zeros(1,N+1);B=[];      % Chute from 87->24
B = A(19+1,:);A(19+1,:) = A(62+1,:) + B;A(62+1,:) = zeros(1,N+1);B=[];      % Chute from 62->19
B = A(60+1,:);A(60+1,:) = A(64+1,:) + B;A(64+1,:) = zeros(1,N+1);B=[];      % Chute from 64->60
B = A(53+1,:);A(53+1,:) = A(56+1,:) + B;A(56+1,:) = zeros(1,N+1);B=[];      % Chute from 56->53
B = A(11+1,:);A(11+1,:) = A(49+1,:) + B;A(49+1,:) = zeros(1,N+1);B=[];      % Chute from 49->11
B = A(26+1,:);A(26+1,:) = A(48+1,:) + B;A(48+1,:) = zeros(1,N+1);B=[];      % Chute from 48->26
B = A(6+1,:);A(6+1,:) = A(16+1,:) + B;A(16+1,:) = zeros(1,N+1);B=[];        % Chute from 16->6

%
% Now change the diagonal when close to the end so that the player loses his turn if (s)he overshoots.
%

for i = 0:SPIN-1
    A(N+1-i,N+1-i) = SPIN-i;
end;

%
% Now normalize A to get the Markov chain.
%

B = A;
A = (1/SPIN)*B;
B = [];

%
% The probability distribution (d), i.e., the probability that game is over by the ith move,
% will approach 1 in the limit as i goes to infinity (ITERATIONS).
%

ITERATIONS = 500;

d=zeros(1,ITERATIONS);
for i = 1:ITERATIONS
    B = A^i;
    d(i) = B(N+1,1);
end;
B=[];

%
% The probability density function is then the difference function on d.
%

prob_density = [0 diff(d)];

figure;
bar(prob_density(1:200));   % graph the density function

figure;
bar(d(1:200))               % graph the distribution function

%
% The expected value (mu) is the inner product of the probability density and the values of the random variable (X),
% which is just the number of turns taken.
%

X = [1:ITERATIONS]';
mu = prob_density * X