%-- SOLUTIONS to Day 4 Exercises
%-- Peter S. Shawhan

%-------------------------------
% Section 1

%-- False negative probability:
p_neg_given_flu = 0.08;
%-- False positive:
p_pos_given_noflu = 0.02;

%-- True probabilities:
p_neg_given_noflu = 1 - p_pos_given_noflu;
p_pos_given_flu = 1 - p_neg_given_flu;

%----- Part 1a
disp('(1a) --- Maryland in wintertime ---');
%-- The prior is p_flu
p_flu = 0.03;
p_noflu = 1 - p_flu;
%-- Calculate the evidence, p_pos :
p_pos = p_pos_given_flu * p_flu + p_pos_given_noflu * p_noflu ;
p_flu_given_pos = p_pos_given_flu * p_flu / p_pos
%-- Now do the same for p_neg :
p_neg = p_neg_given_flu * p_flu + p_neg_given_noflu * p_noflu ;
p_flu_given_neg = p_neg_given_flu * p_flu / p_neg

%----- Part 1b
disp('(1b) --- People sick enough to come to doctor''s office ---');
%-- The prior is p_flu
p_flu = 0.25;
p_noflu = 1 - p_flu;
%-- Calculate the evidence, p_pos :
p_pos = p_pos_given_flu * p_flu + p_pos_given_noflu * p_noflu ;
p_flu_given_pos = p_pos_given_flu * p_flu / p_pos
%-- Now do the same for p_neg :
p_neg = p_neg_given_flu * p_flu + p_neg_given_noflu * p_noflu ;
p_flu_given_neg = p_neg_given_flu * p_flu / p_neg

%----- Part 1c
disp('(1c) --- South Padre Island in summertime ---');
%-- The prior is p_flu
p_flu = 0.001;
p_noflu = 1 - p_flu;
%-- Calculate the evidence, p_pos :
p_pos = p_pos_given_flu * p_flu + p_pos_given_noflu * p_noflu ;
p_flu_given_pos = p_pos_given_flu * p_flu / p_pos
%-- Now do the same for p_neg :
p_neg = p_neg_given_flu * p_flu + p_neg_given_noflu * p_noflu ;
p_flu_given_neg = p_neg_given_flu * p_flu / p_neg

input('End of part 1; Press RETURN to continue');

%-------------------------------
% Section 2

%----- Part 2a
%-- 90% probability of MORE THAN zero events implies 10% probability of
%-- getting EXACTLY ZERO.
%-- We want to find mu for which P(0|mu)=0.10
%--  ==> exp(-mu) = 0.10
disp('(2a) --- 90% UL for N=0');
UL90 = -log(0.10)

%----- Part 2b
%-- Want to find mu which gives 10% probability of getting either 0 or 1:
%-- exp(-mu)+mu*exp(-mu)/1! = 0.10
%--   ==> exp(-mu)*(1+mu) = 0.10
%-- This can't be solved in closed form, so use trial-and-error:
disp('(2b) --- Using trial-and-error:');
mu=3;    c=exp(-mu)*(1+mu); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=3.5;  c=exp(-mu)*(1+mu); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=4;    c=exp(-mu)*(1+mu); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=3.9;  c=exp(-mu)*(1+mu); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=3.88; c=exp(-mu)*(1+mu); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=3.89; c=exp(-mu)*(1+mu); fprintf('  mu=%.3f  --> %.4f\n',mu,c);

%-- Another way: use Matlab's fzero function, which finds a zero of a function.
%-- Do 'help fzero' to read about it.  Here's how to use it in this case:
disp('(2b) --- Using fzero:');
fzero(@(mu) exp(-mu)*(1+mu)-0.1 , 2 )

%----- Part 2c
%-- Want to find mu which gives 10% probability of getting either 0, 1 or 2:
%-- exp(-mu)+mu*exp(-mu)/1!+mu^2*exp(-mu)/2! = 0.10
%--   ==> exp(-mu)*(1+mu+mu^2/2) = 0.10

disp('(2c) --- Using trial-and-error:');
mu=4;    c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=5;    c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=5.2;  c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=5.4;  c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=5.3;  c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=5.31; c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);
mu=5.32; c=exp(-mu)*(1+mu+mu^2/2); fprintf('  mu=%.3f  --> %.4f\n',mu,c);

disp('(2c) --- Using fzero:');
fzero(@(mu) exp(-mu)*(1+mu+mu^2/2)-0.1 , 3 )

%----- Part 2d
%-- 90% probability of LESS THAN one events implies 90% probability of
%-- getting EXACTLY ZERO.
%-- We want to find mu for which P(0|mu)=0.90
%--  ==> exp(-mu) = 0.90
disp('(2d) --- 90% LL for N=1');
LL90 = -log(0.90)

%----- Part 2e
%-- Want to find mu which gives 90% probability of getting either 0 or 1:
disp('(2e) --- Using fzero:');
fzero(@(mu) exp(-mu)*(1+mu)-0.9 , 2 )

input('End of part 2; Press RETURN to continue');

%-------------------------------
% Section 3

%----- Part 3a
ndata = 100;
trueval = 0.3;
xj = trueval + randn(1,ndata);

%----- Part 3b
%-- The 'frequentist point estimate' of the true mean is simply the mean
%-- of the 100 samples.  The distribution of the sample mean is Gaussian
%-- with a sigma of stdev/sqrt(N), where stdev is the standard deviation
%-- for the 100 samples.  The 90% confidence interval is
%-- plus-or-minus 1.65 sigmas around that
sampmean = mean(xj);
sigmean = std(xj,1) / sqrt(ndata);
lowerlimit = sampmean - 1.65*sigmean;
upperlimit = sampmean + 1.65*sigmean;
fprintf('(3b) Mean=%7.4f, confidence interval= [%7.4f,%7.4f]\n', ...
	sampmean, lowerlimit, upperlimit );

%----- Part 3c
%-- Lots more realizations of the noise...
nloop = 10000;
%-- Keep track of how many contain the true value
ncontain = 0;
for iloop=1:nloop
  xj = trueval + randn(1,ndata);
  sampmean = mean(xj);
  sigmean = std(xj,1) / sqrt(ndata);
  lowerlimit = sampmean - 1.65*sigmean;
  upperlimit = sampmean + 1.65*sigmean;
  if lowerlimit <= trueval & upperlimit >= trueval
    ncontain = ncontain + 1;
  end
end
fprintf('(3c) Out of %d trials, %d of the intervals (%.1f%%) contain the true value\n',nloop,ncontain,100*ncontain/nloop);

%----- Part 3d
ndata = 1000;
nloop = 10000;
%-- Keep track of how many contain the true value
ncontain = 0;
for iloop=1:nloop
  xj = trueval + randn(1,ndata);
  sampmean = mean(xj);
  sigmean = std(xj,1) / sqrt(ndata);
  lowerlimit = sampmean - 1.65*sigmean;
  upperlimit = sampmean + 1.65*sigmean;
  if lowerlimit <= trueval & upperlimit >= trueval
    ncontain = ncontain + 1;
  end
end
fprintf('(3d) For n=%d, out of %d trials, %d of the intervals (%.1f%%) contain the true value\n',ndata,nloop,ncontain,100*ncontain/nloop);

input('End of part 3; Press RETURN to continue');


%-------------------------------
% Section 4

%----- Part 4a
ndata = 100;
sigma = 1.0;
xj = trueval + sigma*randn(1,ndata);

%----- Part 4b
%-- Consider different possible values of s
sarray = (-1.0:0.01:1.0);

%-- Calculate the likelihood
L = zeros(size(sarray));
for is=1:length(sarray)
  s = sarray(is);
  L(is) = (2*pi)^(-ndata/2) * sigma^(-ndata) * exp(-sum((xj-s).^2/sigma^2)/2);
end
plot(sarray,L);
input('(4b) Likelihood plot; press RETURN to continue');

%----- Part 4c
prior = zeros(size(sarray));
prior(find(sarray>=0)) = 1;
evidence = sum(L.*prior);
posterior = L.*prior/evidence;
%-- Check that the integral of the posterior is 1
sum_posterior = sum(posterior)
plot(sarray,posterior);
input('(4c) Posterior plot; press RETURN to continue');

%----- Part 4d
%-- Find the mode
imax = find(posterior==max(posterior))
%-- Extend symmetrically until we enclose 90% of the area
addbins = 0;
while sum(posterior(imax-addbins:imax+addbins))<0.90
  addbins = addbins+1;
end
fprintf( '(4d) Mode=%7.4f , interval = [%7.4f,%7.4f]\n', ...
        sarray(imax), max(sarray(imax-addbins),0), sarray(imax+addbins) );

%----- Part 4e
fprintf('(4e) Repeat those steps with different realizations of the data...\n');
for iloop=1:10
  %-- Generate the data
  xj = trueval + sigma*randn(1,ndata);
  %-- Calculate the likelihood as function of s
  for is=1:length(sarray)
    s = sarray(is);
    L(is) = (2*pi)^(-ndata/2) * sigma^(-ndata) * exp(-sum((xj-s).^2/sigma^2)/2);
  end
  %-- Calculate the posterior
  evidence = sum(L.*prior);
  posterior = L.*prior/evidence;
  %-- Find the 90% interval
  imax = find(posterior==max(posterior));
  addbins = 0;
  while sum(posterior(imax-addbins:imax+addbins))<0.90
    addbins = addbins+1;
  end
  fprintf( '     Mode=%7.4f , interval = [%7.4f,%7.4f]\n', ...
        sarray(imax), max(sarray(imax-addbins),0), sarray(imax+addbins) );
end
