% Lecture 7 of PHYS 399
%% part 1
nu = 0:10;
mu = 3;
N = 11;
[Ps] = poisson(nu, mu);
% [PsML] = poisspdf(nu,mu); % using matlab stat TB builtin
 
bar(nu,100.*Ps,0.2)
% bar(nu,100.*PsML,0.2,'r')
grid on
xlabel(' \nu--> ','FontSize',16)
yl = [' P_{',num2str(mu),'}(\nu) (%)']
ylabel(yl,'FontSize',16)
title('Poisson distribution for \mu=3','FontSize',16)
hold on
plot(nu,100.*Ps,'b-')


%% part 2
figure
nu = 0:10;
mu = 5;
N = 11;
[Ps] = poisson(nu, mu)
 
bar(nu,100.*Ps,0.2)
grid on
xlabel(' \nu--> ','FontSize',16)
yl = [' P_{',num2str(mu),'}(\nu) (%)']
ylabel(yl,'FontSize',16)
title(['Poisson distribution for \mu=',num2str(mu)],'FontSize',16)
hold on
plot(nu,100.*Ps,'r-')


%% part 3
% QC 11.1  18 hens, each lays 1 egg / day on average.  Check
% henhouse every hour, want to know the probabilities of finding
% 0, 1, 2, 3, and 4 eggs in that hour.

% 18 eggs / 24 hours, means rate is 18/24 = 0.75 (per hour).  
% This is my mu. so calculate Ps for mu = 0.75, plot
figure
[Ps] = poisson(nu, 0.75)
 
bar(nu,100.*Ps,0.4)
grid on
xlabel(' \nu--> ','FontSize',16)
ylabel(' P_{\mu}(\nu) (%)','FontSize',16)
title('Poisson distribution for \mu=0.75','FontSize',16)

% about a 47% chance I find no eggs, 35% chance of finding 1, etc.
% (see bar plot).  Note, this assumes egg-laying rate is constant through
% the day, it is NOT!


%% part 4
% *** fit of Poisson dist. to observed counting data:
% data is in cnt (5 second intervals)
% first make histogram of cnt
load ~/physics/Lab_Courses/PHYS391/MATLAB/cnt

hcnt = hist(cnt);  % *** make histogram of count data
nul = min(cnt); % lowest value of counts observed
nuh = max(cnt); % highest "     "    "     "
nu = [nul:nuh]; % range of nu's to use for this
hcnt = hist(cnt,nu)  % get # observations (each per 5 seconds) per nu 
bar(nu,hcnt)  % plot results using bar instead of hist plot, control!
grid on
xlabel(' \nu--> ')
ylabel(' counts')
fid = gcf;  %  this returns a pointer to the current figure, to be used later

%% OK, have our histogram, now want to fit it using 'brute-force' method
% *** eyeball that mu is somewhere around 15.

% *** step from nul to nuh in 0.1 steps for mu, run
%     poisson with this mu (and standard nu), then
%     calculate un-normalized chi (chi squared) between
%     poisson prediction for this mu and data (hcnt).
i = 1;
N = sum(hcnt);  % *** for later 'anti-normalization' ***
% *** what does this next for loop do? ***
for mu=nul:0.1:nuh
    [Ps] = poisson(nu,mu);
    mus(i) = mu;
    NPs = N.*Ps;  % normalized Ps's 
    chi(i) = sum((NPs-hcnt).^2);
    i = i + 1;
end

%% after stepping through trial mu's, look for minimum misfit (chi-squared)
% for a given mu
% *** plot chi (squared) vs. mus, look for minimum***
figure  % *** new call to 'figure' opens new figure window for plotting
plot(mus,chi,'b*')
grid on
xlabel(' \mu --> ')
ylabel(' chi squared misfit (unnormalized)')
hold on
fig2 = gcf;

% *** see a min at about 18 or so.

%% rather than eye-balling mu with minimum misfit...
% a better way than looking!, 'find' where chi == min(chi)
ipt = find(chi==min(chi)); % *** here's our all-powerful 'find' command
mu = mus(ipt)  % *** how is best-fit mu found in this way ***

plot(mu,chi(ipt),'r*','MarkerSize',12)

%% OK, our 'better way' gives the mu that best fits our observed counts
% let's replot the Poisson dist. with best fit over histogram
figure(fid); % *** this recalls the previous bar/hist plot figure
% we now will plot on our original figure ***
% mu = 8.9;
[Ps] = poisson(nu,mu);  % *** what does this call do? ***
hold on
plot(nu,N.*Ps,'r-','LineWidth',2)
yl = [' P_{',num2str(mu),'}(\nu) (red line)']
ylabel(yl)





