hw.m

% 7.1 4
if false
px = [1/4 1/4 1/2];
n = 10;
d = nconv(px, n);

x = 1:size(d,2);
disp(sprintf('The distribution for %1.0f games is:', n))
plot(x, d, 'o')
disp(d')
end

% 7.1.8
if false
p = [.4 .2 .2 .1 .1];

% part 1 
% distributions for different series
[d4, x4] = nconv(p,4);
[d6, x6] = nconv(p,6);
[d7, x7] = nconv(p,7);

clf;
subplot(1,3,1);
plot(x4,d4, 'o')
title('4 game series');

subplot(1,3,2);
plot(x6,d6, 'o')
title('6 game series');

subplot(1,3,3);
plot(x7,d7, 'o')
title('7 game series');

% part 2
% p(avg > .4)

num_games = 4;
num_hits_per_game = 4;
batting_avg = .4;
num_hits = ceil(num_games*num_hits_per_game*batting_avg);
prob = sum(d4(num_hits:end));
disp(sprintf('%2.2f%% of the time the player hits above a %2.0f in a %1.0f game series', prob, batting_avg*1000, num_games))

% part c
% long term batting avg
avg = sum([0 1 2 3 4] .* p)/5;
disp(sprintf('His longterm batting average is %0.3f', avg))
end

% 7.2.15
if true

reps = 10000;
n = 1000;

% count number of heads & tails
heads = @(p) sum(rand(n, reps)<p);
fair_heads = heads(1/2);
weighted_heads = heads(1/3);

% calculate Zs
z = @(x) (x - n*p).^2 ./ (n*p);
fair_z = z(fair_heads) + z(n-fair_heads);
weighted_z = z(weighted_heads) + z(n-weighted_heads);

% plot results
clf;
subplot(2,1,1);
hist(fair_z,100);
title('fair coin');

subplot(2,1,2);
hist(weighted_z,100);
title('weighted coin');

end

% 9.1 1
if false
  
n = 100;
p = .5;
z = @(x) (x-.5-n*p)/sqrt(n*p*p);

% P(Sn <= 45)
x1 = 45;
z1 = z(x1);
p1 = normcdf(z1);
disp(sprintf('P(Sn<= %2.0f) = %0.3f', x1, p1))

% P(45 <= Sn <= 55)
x2 = 55;
z2 = z(x2);
p2 = (normcdf(z2)-.5)*2;
disp(sprintf('P(%2.0f <= Sn <= %2.0f) = %0.3f',45, x2,p2))

% P(Sn >= 63)
x3 = 63;
z3 = z(x3);
p3 = 1 - normcdf(z3);
disp(sprintf('P(Sn >= %2.0f) = %0.3f',x3,p3))

% P(Sn <= 57)
x4 = 57;
z4 = z(x4);
p4 = normcdf(z4);
disp(sprintf('P(Sn <= %2.0f) = %0.3f',x4,p4))

  
end

% 9.1 3
if false
n = 48;
x = 30;
z = @(x, p) (x-.5-n*p)/sqrt(n*p*(1-p));

pj = .75;
zj = z(x, pj);
pj = 1 - normcdf(zj);

pa = .5;
za = z(x, pa);
pa = 1 - normcdf(za);

disp(sprintf('June passes %2.2f%% of the time', pj*100))
disp(sprintf('April passes %2.2f%% of the time', pa*100))
end

% 9.1.5
if false
n = 300;
p = .3;
pact = .267;
x = ceil(pact*n);
z = (x-.5-n*p)/sqrt(n*p*(1-p));
p1 = normcdf(z);  
disp(sprintf('The rookie should miss %2.0f or fewer balls %2.2f%% of the time.\nThus he probably is worse than a 300 hitter.', x,p1*100 ))
end

% 9.1.7
if false
  
n = 1750;
p = .6;
x = 1061;
z = (x-.5-n*p)/sqrt(n*p*(1-p));

p1 = 1 - normcdf(z);
disp(sprintf('Theres a %2.2f%% chance that %4.0f or more students accept', p1*100, x))
  
end

% 9.2 1
if false
px = ones(1,6)*1/6;
n = 24;
x = 84;

mu = 7/2;
s2 = 35/12;
z = (x+1-.5 - n*mu) / sqrt(n*s2);
p1 = 1 - normcdf(z);
disp(sprintf('P(S>84)= %0.5f', p1))

% doesn't match
[d, x] = nconv(px, n);
disp(sprintf('P(S=84)= %0.5f', d(85)))

end

% 9.2.4
if false
n = 1000;
mu = 4.5;
s2 = sum([[0:9] - mu] .^ 2)/10;
std3 = sqrt(s2)*3/sqrt(n);

% one set of a 1000 rand digits
x = randint(0,9,n);
avg = sum(x)/n;
disp(sprintf('The average is %1.3f, which is %0.3f away from the mean and falls within 3 standard deviations', avg, mu-avg))

% m sets of 1000 rand digits
m = 10000;
x2 = randint(0,9,n,m);
avg2 = sum(x2) ./ n;
p = sum((mu-avg2)<std3) / m ;

disp(sprintf('Yes. %2.2f%% of the averages fell within 3 standard deviations.', p*100))

end

nconv.m

function [d, x] = nconv(p,n)
  d = p;
  for i = 2:n
    d = conv(d,p);
  end

  x = 1:size(d,2);

randint.m

function x = randint(a,b,n,m)
  
  if nargin < 4
      m = 1;
  end

  x = floor(rand(n,m)*abs((b-a+1))) + a;