%
% Plots a "binomialness" plot
% 
% Problem Epage 201 
% Example Epage 142
%
% 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 biology;   % <- freqs/numpaps

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

%--
% 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,n);
for ti=1:n,
  nck      = nchoosek(n,ti);
  phik(ti) = log( (n_k_star(ti)+eps)/(N*nck) ); 
end
figure; plot( k, phik, 'x' ); grid on; 
xlabel( 'k' ); ylabel( '\phi(n_k)' ); 
title( 'binomialness plot for the biology data' ); 
saveas( gcf, 'biology_binomialness_plt', 'epsc' );