function  [L,chisq] = binomial_logL(N,Po,Pe)

%BINOMIAL_LOGL -- Binomial log-likelihood goodness-of-fit measure.
%
% [L,chi2] = binomial_logL(N,Po,Pe)
%
% N is an non-negative integer array of arbitrary shape. It specifies the
% number of replications in each "bin" I.  PO is an array of observed
% probabilities (=proportions); PE is an array of expected probabilities.
% When not supplied, PE defaults to PO. It must have the same shape as N.
%
% BINOMIAL_LOGL calculates the log-likelihood to observe PO if the data
% come from a (family of) binomial distribution(s) governed by N and PE:
%   L = sum(N.*Po.*log(Pe) + N.*(1-Po)*log(1-Pe))
% This differs from the true binomial log-likelihood by the additive
% constant sum(nchoosek(N,N.*Po)). Note also that lim[p*log(p)]=0, p->0.
%
% If a second output argument is requested, the chi-square statistic is
% calculated too. Note that chi2<>0 only when Po<>Pe. Specifically:
%   chi2=-2*log[L1/L0]=-2*(binomial_logL(N,Po,Pe)-binomial_logL(N,Po,Po))
% It is asymptotically chi-square with degrees of freedom equal to the
% difference between the number of bins #(N) and the number of parameters
% specifying the expected probabilities Pe:  df = #(N)-#(par(Pe))
%
% References:
%  Watson, A. B. (1979). Probability summation over time.
%    Vision Research 19, 515-522. [Especially the Appendix]
%  Hoel, P. G., Port, S. C., & Stone, C. J. (1971).
%    Introduction to Statistical Theory. Houghton Miffin, Boston.
%
% Example 1:  Suppose we have an experiment with three conditions:
%  1=control, 2=down, and 3=up.  We have 50 replications in the
%  control condition and 25 reps in each of the other two. Suppose the
%  observed frequencies are Po=[.48 .24 .72]. Finally, let our theory
%  predicts that Pe(1)=.5, Pe(2)=.5-p, Pe(3)=.5+p.  Then for p=0.3:
%    [L,chi2] = binomial_logL([50 25 25],[.48 .24 .72],[.5 .2 .3])
%    L = -72.7219 , chi2 = 19.0075
%  We can search for the parameter p that minimizes L (see FMINSEARCH).
%
% See also ENTROPY, CHI2, WEIBULLFIT, BINOFIT, NCHOOSEK, XTAB.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2004/03/22 11:00 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

%-- Expand data, unifying P and 1-P into a long column-vector.
N2  = [N(:) ; N(:)] ;       % assume integers >=0, no error checking
Po2 = [Po(:) ; 1-Po(:)] ;   % assume >=0, <=1
if (nargin==3)   % Expected Pe supplied explicitly
  Pe2 = [Pe(:) ; 1-Pe(:)] ;
  chisq = -1 ;   % recompute later, if requested
else  % Pe defaults to Po, for a null hypothesis
  Pe2 = [Po(:) ; 1-Po(:)] ;
  chisq = 0 ;    % since L1==L0 when Pe==Po
end
%disp([N2 Po2 Pe2 log(Po2) log(Pe2)]) ;

%-- Log-likelihood
%Note: As lim[p*log(p)]=0 when p-->0, zero entries don't count.
idx = find(Po2>0 & Pe2>0) ;   % avoid log(0)
L = sum(N2(idx).*Po2(idx).*log(Pe2(idx))) ;  % [1x1]

%-- Chi-square, if requested
if (nargout==2 & chisq<0)
  L0 = sum(N2(idx).*Po2(idx).*log(Po2(idx))) ;
  chisq = 2*(L0-L) ;
end

%%*   Return L and possibly CHISQ
%%%%% End of BINOMIAL_LOGL.M
