Page 71
roll=ones(1,11);
roll=cumsum(roll);
roll=roll+1;
prob=[1 2 3 4 5 6 5 4 3 2 1]/36;
bar(roll,prob)

                   Page 72
weightedRoll=prob.*roll;

                   Page 73
m=sum(weightedRoll)
v=(roll-m).^2; % sum of squared deviations
var=sum(v.*prob)

                   Page 76
% use the previous m-file, gaussian.m:
   %   % function y=gaussian(x,m,s);
   %   % m=mean, s=stddev
   %   % note 1/sqrt(2*pi)=.3989422803
   % y=(.3989422803/s)*
   % exp(-0.5*((x-m)./s).^2);
x=[-10:.1:10];
y=gaussian(x,30,sqrt(15)); plot(x,y)


                  Exercises/Experiments
1.
% No built-in combinatorics in matlab but it
% is easy to do factorials and hence
% permutations and combination calculations
% permutations of 6 things taken 3 at a time
n6=1:6; n3=1:3;
perm6t3=prod(n6)/prod(n3)

% combinations of 6 things taken 3 at a time
comb6t3=perm6t3/prod(n3)

2.
% combinations of 5 things take 2 at a time
comb5t2=prod(1:5)/(prod(1:2)*prod(1:3))
comb5t3=prod(1:5)/(prod(1:3)*prod(1:2))
perm5t3=prod(1:5)/prod(1:2)

3.
% rand(1,300) is a random vector with components
% between 0 and up to but not including 1; then
% 6 times this gives numbers from 0 up to 6, add
% 1 and get numbers 1 up to 7, finally fix()
% truncates the fractional part
die=fix(6*rand(1,300)+1);
% now count the number of 3's
count3s=1./(die-3);  % gives infinity at every 3
count3s=isinf(count3s); % 1 for infinity, 0 otherwise
number3s=sum(count3s)

4.
red=fix(6*rand(1,360)+1);
blue=fix(6*rand(1,360)+1);
pairDice=red+blue;
x=2:12;
hist(pairDice,x)
hold on
h=hist(pairDice,x)
mu=dot(x,h)/sum(h)
% wt each int by its frac. of outcomes, add
v=(x-mu).^2; % vector of diffs squared
var = dot(v,h)/sum(h); % variance
sigma=sqrt(var)
t=linspace(2,12);
y=360*exp(-(t-mu).^2/(2*sigma^2))/(sigma*sqrt(2*pi));
plot(t,y)
% the theoretical probability for seeing 2 is
% 1/36, same for 12, for seeing 3 is 2/36,
% same for 11, etc, for seeing 7 is 6/36.
% compare with h above
theory=[10 20 30 40 50 60 50 40 30 20 10];
mu=dot(x,theory)/sum(theory)
v=(x-mu).^2; var=dot(v,theory)/sum(theory);
sigma=sqrt(var)
y=360*exp(-(t-mu).^2/(2*sigma^2))/(sigma*sqrt(2*pi));
hold off
plot(t,y); hold on
hist(pairDice,x)

5a.
% # ways for two blue eyed is C(2,2)
% (2 choose 2)=2!/(2!*0!) so the
% probability is that times (1/4)^2, etc.
blublu=prod(1:2)/(prod(1:2)*1)*(1/4)^2 % blu/blu children
Bwnblu=prod(1:2)/(prod(1:1)*prod(1:1))*(3/4)*(1/4) % Bwn/Blu
BwnBwn=prod(1:2)/(1*prod(1:2))*(3/4)^2 % Bwn/Bwn children
blublu+Bwnblu+BwnBwn

5b
% exactly two are brown eyed is 5 choose 2
exact2=prod(1:5)/(prod(1:2)*prod(1:3))*(3/4)^2*(1/4)^3}

5c.
exact3=prod(1:5)/(prod(1:3)*prod(1:2))*(3/4)^3*(1/4)^2
exact4=prod(1:5)/(prod(1:4)*prod(1:1))*(3/4)^4*(1/4)^1
exact5=prod(1:5)/(prod(1:5)*1)*(3/4)^5
atleast2=exact2+exact3+exact4+exact5