clear
echo on
%Solution to problem 12
%In this problem we explored the Poisson probability
%distribution.  This distribution governs random
%processes such as radioactive decay.  The distribution
%is governed by:
%    p(k) = (r*t)^k/k!*exp(-r*t)
%where r is the decay rate and t is the time of
%observation.
%
pause
%For part a) we wish to calculate the mean and variance
%of the distribution.  We use the attached program which
%must be saved under the name poisson_prob.m
%So:
rt=5;
k=[0:50];
mean=k*poisson_prob(rt,k)'
var=(k-mean).^2*poisson_prob(rt,k)'
pause
%Note that both the mean and the variance of the Poisson
%distribution are just r*t!  This holds true for the
%larger value of r*t as well:
rt=500;
k=[350:650];
mean=k*poisson_prob(rt,k)'
var=(k-mean).^2*poisson_prob(rt,k)'
%It probably would have been alot faster to use the 
%asymptotic formula for k! for these large values.
pause
%
%Now let's do part b) and compare these distributions
%to the gaussian distribution.  The gaussian distribution
%is given by exp(-(x-mean)^2/(2*sigma^2))/sqrt(2*pi*sigma^2).
%Recalling that the mean and variance of the Poisson
%distribution is just r*t, we get:
figure(1)
rt=5;
k=[0:10];
plot(k,poisson_prob(rt,k),'o')
hold on
x=[0:.1:10];
pgauss=exp(-(x-rt).^2/2/rt)/sqrt(2*pi*rt);
plot(x,pgauss)
xlabel('x')
ylabel('probability density')
title(['Probability Distribution for r*t = ',num2str(rt)])
hold off
pause
%And we have the similar result for the larger r*t:
figure(2)
rt=500;
k=[450:550];
plot(k,poisson_prob(rt,k),'o')
hold on
x=[450:550];
pgauss=exp(-(x-rt).^2/2/rt)/sqrt(2*pi*rt);
plot(x,pgauss)
xlabel('x')
ylabel('probability density')
title(['Probability Distribution for r*t = ',num2str(rt)])
hold off
pause
%Comparison of the two graphs shows that the larger
%value of r*t is well fit by a gaussian distribution.
%
pause
%Now for the final part (c):
%We want to have an accuracy of 10% at the 95% confidence
%level.  This means that two standard deviations should
%equal 10% of the total number of counts.  Thus:
%
%   2*(r*t)^.5 = 0.10*r*t
%
%or,
%   t = 400 / r = 13.3 minutes
%
%which is the desired result.
echo off



_________________________________


function p=poisson_prob(rt,k)
%This function calculates the poisson probability
%for a given decay rate times time rt and observed
%decays k.  It is set up so that k is a vector.
%To avoid overflow/underflow problems we calculate
%the log of p and then return the exponential of
%this value.

%We need to calculate the log of k factorial.  This
%must be done up to the maximum value of k.  We only
%need to do this once.  Thus:

logfac=[0,cumsum(log([1:max(k)]))];

%Actually, in this line logfac(n) is (n-1)! so as to
%avoid indices of zero if we feed it k=0.  The rest
%of the program is now just one line!

p=exp(k.*log(rt)-logfac(k+1)-rt);

%Note that we are just using the elements of the
%vector logfac which correspond to the values of
%k for which we are calculating the probability.
%The others are not used.