function exercise_C_3_3(nTrials)
% EXERCISE_C_3_3 - Generates the probabilities for exercise C-3.3 via Monte-Carlo trials
% 
% Written by:
% -- 
% John L. Weatherwax                2007-05-04
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% unif_discrete_sample.m 
BayesianNetworks;
  
% calcuate the analytical formulas: 
% 
format rat; 
p0 = (5/6)^3;
p1 = 3*(1/6)*((5/6)^2); 
p2 = 3*((1/6)^2)*(5/6);
p3 = 1/(6^3); 
p0 + p1 + p2 + p3 

% the analytical expectation: 
% 
e = -p0 + p1 + 2*p2 + 3*p3;

% the Monte-Carlo calculation (begins here): 
% 
format short; 

%nTrials = 5; 

% select a choice for the special number (the one picked by the player)
myDieChoice = 3; 

% simulate nTrials number of three die rolls: 
d1 = unif_discrete_sample(6,1,nTrials); 
d2 = unif_discrete_sample(6,1,nTrials); 
d3 = unif_discrete_sample(6,1,nTrials); 

%--
% do statistics on the payoff for each trial:
%--
fprintf('Printing the true probabilities and the Monte-Carlo approximations\n'); 

% the number of times we have NO matching die: 
payMO = length( find( (d1~=myDieChoice) & (d2~=myDieChoice) & (d3~=myDieChoice) ) ); 
fprintf('p0=%6.4f; payMO/nTrials=%6.4f; diff=%10.8e\n',p0,payMO/nTrials,p0-(payMO/nTrials)); %payMO/nTrials

% the number of times we have ONE matching die: 
r1 = (d1==myDieChoice) & (d2~=myDieChoice) & (d3~=myDieChoice); % <- we match roll number 1
r2 = (d1~=myDieChoice) & (d2==myDieChoice) & (d3~=myDieChoice); % <- we match roll number 2
r3 = (d1~=myDieChoice) & (d2~=myDieChoice) & (d3==myDieChoice); % <- we match roll number 3
payPO = length( find( r1 | r2 | r3 ) );
fprintf('p1=%6.4f; payPO/nTrials=%6.4f; diff=%10.8e\n',p1,payPO/nTrials,p1-(payPO/nTrials)); %payPO/nTrials

% the number of times we have TWO matching die: 
r12 = (d1==myDieChoice) & (d2==myDieChoice) & (d3~=myDieChoice); % <- we match roll number 1 & 2 
r13 = (d1==myDieChoice) & (d2~=myDieChoice) & (d3==myDieChoice); % <- we match roll number 1 & 3 
r23 = (d1~=myDieChoice) & (d2==myDieChoice) & (d3==myDieChoice); % <- we match roll number 2 & 3 
payPT = length( find( r12 | r13 | r23 ) );
fprintf('p2=%6.4f; payPT/nTrials=%6.4f; diff=%10.8e\n',p2,payPT/nTrials,p2-(payPT/nTrials)); %payPT/nTrials

% the number of times we have THREE matching die:
r123 = (d1==myDieChoice) & (d2==myDieChoice) & (d3==myDieChoice); % <- we match roll number 1, 2 & 3 
payP3 = length( find( r123 ) );
fprintf('p3=%10.8f; payP3/nTrials=%10.8f; diff=%10.8e\n',p3,payP3/nTrials,p3-(payP3/nTrials)); %payP3/nTrials

mce_approx = -(payMO/nTrials) + (payPO/nTrials) + 2*(payPT/nTrials) + 3*(payP3/nTrials);
fprintf('e=%10.8f; MC_approx=%10.8f; diff=%10.8e\n',e,mce_approx,e-mce_approx);