function [Pexact,Pnumber,Pest]=craps(number,trials)
% format [Pexact,Pnumber,Pest]=craps(number,trials)
% input  number=value of 1st roll
%        trials=Number of rolls for a Monte Carlo simulation of probabilities
%               1000 would be an appropriate number, 10000 for 3 signicant figures
%              [optional]
% output: Pexact= Exact probability of winning, given 1st roll equal to number
%        Pnumber=Probability of rolling that number by rolling 2 dice.
%        Pest   =Estimated probability from Monte Carlo simulations
%
% calls absorb.m, an absorbing Markov chain solver
% This problem is solved more simply by crapsmof65.m, based on Mosteller's
% elegant simple solution, see crapsmof65.m
% absorb.m written by Eugene.Gallagher@umb.edu
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% written 6/2/03 by E. Gallagher for EEOS601, revised 1/12/11

die=1:6; % create a sample space of two dice rolls
die1=repmat(die,6,1); % create a matrix, with each row equal to 1 to 6;
die=die';  % transpose of die, now a column vector
die2=repmat(die,1,6); % create a matrix, with each column equal to 1 to 6;
sumdice=die1+die2;  % This gives the 36 sums of two die.

i=find(sumdice==number);  % find the indices for each value equal to number;
lengthi=length(i);     % find how many sums of 2 die equal the number in input;
if lengthi==0          % only true if the number isn't in the sample space.
    disp('Your number isn''t possible, not in the sample space')
    number
    % The following odd syntax will print the number
    fprintf('Your number, %6.0f, isn''t possible, not in the sample space\n',number)
    return % This stops the program after displaying the message;
end

[r,c]=size(sumdice); % find size of sample space
Pnumber=lengthi/(r*c); % What is the probability of rolling your number with 1 roll
                       % of two dice?
                       
% Is the craps shooter an immediate winner?
if (number==7 | number==11) % 2 logical statements using the | OR command.
    Pexact=1;              % Matlab performs statements in this if section only
                           % if the logical statement number equals 7 or equals 11
                           % is true.
    if nargout>1 % did the user request more than 1 output variables? (nargout stands for
        Pest=1;  % number of OUTPUT arguments requested in the call to the function.
    end
    return % exit the function and return to the command window.
end

% Is the craps shooter an immediate loser?, use cut and paste the previous code & modify
if (number==2 | number==3 | number==12)
    Pexact=0;
    if nargout>2 % did the user request 3 outputs (nargout stands for
        Pest=0;  % number of OUTPUT arguments requested in the call to the function.
    end
    return  % exit the function and return to the command window.
end

%                   Making a point calculations
               
% Exact probability of making a point
%  First find probability of rolling a 7, which results in a loss.
j=find(sumdice==7);P7=length(j)/(r*c); % P7 is the probability of rolling a 7;

% Advanced topic: solve for the exact P using a Markov chain model:
% Note that Mosteller's 50 problem book, problem 9, gives a simpler solution.
% (Gallagher programmed as crapmof65.m
% The exact probability can be modeled as a 3-state Markov chain.
% Markov chains are used to solve 'the Gambler's ruin problem' also called
% the drunkard's walk problem. There are 3 states to the problem:
%   Winning
%   Losing
%   Roll again
% In a gambler's ruin problem, when you've won, the game is over
%                              when you've lost, the game is over
%                              in some outcomes, you play again.
% The probabilities of these events are placed in a square matrix, which you
% read from rows to columns
% The three states in words are
% Probability of outcomes on a 2nd roll of the dice:
% Poutcome=
%                           To:
%         F              _ Won    Lost       Roll Again _
%         R  Won        |  1.0      0            0       |
%         O  Lost       |   0      1.0           0       |
%         M  Roll Again |_Pnumber   P7    1-(Pnumber+P7)_|
% Each row in the matrix expresses the probability of an event after the event
% in the row has happened. 
% If you made your point in the previous through, you've won, game over
%                P (Winning, given that you've won) is 1.0
% If you rolled a 7 on your previous toss of the dice, you've already lost
%                P (Losing, given that you rolled a 7 on your previous toss) is 1.0
% All of the information in this matrix is in the final row. After you've rolled 
% a 4, 5, 6, 8, 9, or 10, you must make your point before rolling a 7
%         The probability of rolling your number with one roll of the dice has
%         already been calculated. and it is Pnumber
%         The probability of losing when rolling a 7 is the probability of rolling
%         a 7 or P7
%         Using the idea of the complement of an event, the probability of rolling
%         again is 1 - (Probability of winning on one roll + Probability of losing on one roll)
% absorb.m, based on a remarkable book by John Kemeny & J. Laurie Snell will
% calculate the full set of probabilities for this model.
Ppoint=[1 0 0;0 1 0;Pnumber P7 1-(Pnumber+P7)];
[tau,N,B]=absorb(Ppoint);
% the B matrix gives the probability of winning and losing over the long run
% This calculates the exact probability for a potentially countably infinite sequence.
Pexact=B(1);  % This gives the exact probability of making your point.
              % The tau matrix says how long you will play on average before
              % making your point or losing (we don't need that information here).
 % Nahin in 'Duelling idiots' solves two similar problems
 %           by solving for the sum of an infinite geometric 
 %           series (see index in Nahin, sum geometric series for method)

% We can estimate the exact probability using a Monte Carlo simulation
% Now let's do a Monte Carlo simulation of this event if the user requests it
% by entering the number of trials
tallywins=0;    % This variable will keep track of the number of successes.
if nargin > 1  | nargout>2 % Only do this if the user has requested a number of trials
    if nargin<2
        % Theorem 5.3.2 can be used to calculate how many Monte Carlo trials are
        % needed to produce a result such that half the 95% CI is equal to d.
        % For this Craps problem, we'll just calculate d so that it is within
        % 0.01 of the exact value 95% of the time.
        d=0.01;trials=LMtheorem050302(0.05,0.01,Pexact)
    end
    for i=1:trials  % do this for as many times as requested.
        stillplaying=1;   % you continue playing, stillplaying=1, 
                          % until you lose or make your point.
        while stillplaying==1
            % generate a single random number
            r=rand(1);    % generates a uniformly distributed random number between
                          % 0 and 1
            if r<=Pnumber % if the random number is less than or equal to your
                          % probability of making your point, you win
                tallywins=tallywins+1;  % tally your win in the win column.
                stillplaying=0; % you've won, play again
            elseif r<=(Pnumber+P7) % if you haven't won, you could lose by rolling a 7
                stillplaying=0;
            end
        end
    end
    Pest=tallywins/trials;  % calculate your estimated probability of making your point.
end

function  [tau,N,B,CP,Q,R,tau2,N2,H,MR,VR]=absorb(P);
% format: [tau,N,B,CP,Q,R,tau2,N2,H,MR,VR]=absorb(P); 
% input: P an absorbing Markov transition matrix in 'from rows
%        to column' form.  Rows are probability vectors and
%        must sum to 1.0.  At least one main diagonal element=1
% output: tau: Time before absorption from a transient state.
%        B: P(ending up in absorbing state from a transient state)
%        tau2: variance matrix for tau;N: fundamental matrix
%        N2: covariance matrix for N; 
%        CP: transition matrix in canonical form.
%        Q: Submatrix of transient states, R: from transient to 
%        absorbing state submatrix,
%        H: hij the probability process will ever go to transient 
%           state j starting in transient state i (not counting the
%           initial state (K & S, p. 61))
%        MR: expected number of times that a Markov process remains
%            in a transient state once entered (including entering step)
%        VR: Variance for MR (K & S, Theorem 3.5.6, p. 61)
% written by E. Gallagher, ECOS, now EEOS
% UMASS/Boston email: Eugene.Gallagher@umb.edu
% written 9/26/93, revised: 10/26/94
% refs: 
% Kemeny, J. G. and J. L. Snell.  1976.  Finite Markov chains.  
%        Springer-Verlag, New York, New York, U.S.A. 
% Roberts, F. 1976.  Discrete Mathematical Models with applications
%        to social, biological, and environmental problems. Prentice-Hall
%
[pdim,pdim]=size(P);
if sum([sum(P')-ones(1,pdim)].^2)>0.001
      error('Rows of P must sum to 1.0')
end
dP=diag(P);
ri=find(dP==1.0);
if isempty(ri)
       error('No absorbing states (reenter P or use ergodic.m)')
end
rdim=length(ri);
qi=find(dP~=1.0);qdim=length(qi);
I=eye(qdim,qdim);
Q=P(qi,qi);
N=inv(I-Q);   % the fundamental matrix
tau=sum(N')';
CP=P([ri' qi'],[ri' qi']);  % the canonical form of P
R=CP(rdim+1:pdim,1:rdim);
B=N*R;
if nargout>6          % only if variances requested 
     Ndg=diag(diag(N));
     N2=N*(2*Ndg-I)-N.^2;
     tau2=(2*N-I)*tau-tau.^2;
     H=(N-I)/Ndg;
     dQ=diag(Q);oneq=ones(qdim,1);
     MR=oneq./(oneq-dQ);
     VR=dQ./(oneq-dQ).^2;
end