%
% Plots a "poissonness" plot
% 
% Problem Epage 201 
% Example Epage 133
%
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');
load counting;   % <- freqs/numpaps

k   = count;             % <- k is an array of the number of counts 0:L
L   = max(k);            % <- the maximum number of counts ever recieved
n_k = freqs;             % <- the number of trials that resulted in "k" successes
N   = sum(n_k);          % <- the total number of successes ... 

phik = zeros(1,L+1);
for ti=1:(L+1),
  % don't do the fancy phik adjustment of Hoaglin and Tukey [1985] ... take n_k as is
  phik(ti) = log( factorial(ti-1)*(n_k(ti)+eps)/N ); 
end

figure; plot( k, phik, 'x' ); grid on; 
xlabel( 'k' ); ylabel( '\phi(n_k)' ); 
title( 'poissonness plot for the counting data' ); 

%--
% adjust n_k as suggested by Hoaglin and Tukey: 
%--
n_k_star = n_k; 
indeq1 = find( n_k==1 ); 
n_k_star(indeq1) = 1/exp(1); 
indgt2 = find( n_k>=2 );
n_k_star(indgt2) = n_k(indgt2)-0.67-0.8*n_k(indgt2)/N;

phik = zeros(1,L+1);
for ti=1:(L+1),
  % please DO the fancy phik adjustment of Hoaglin and Tukey [1985] ... DON'T take n_k as is
  phik(ti) = log( factorial(ti-1)*(n_k_star(ti)+eps)/N ); 
end

figure; plot( k, phik, 'x' ); grid on; 
xlabel( 'k' ); ylabel( '\phi(n_k)' ); 
title( 'poissonness plot for the counting data' ); 
saveas( gcf, 'counting_poissonness_plt', 'epsc' );