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