%-----------------------------------------------------------------------------
%
% $Id: lrtest.m,v 1.6 1999/09/22 03:29:26 philipc Exp $
%
% Program to calculate the goodness-of-fit of a sample to a normal
% distribution with the same mean and standard deviation
%
% Inputs:
%     x - array of length k giving center of each bin
%     n - number of points in each bin
%     mu, sigma - standard deviation and mean to test against
% 
%-----------------------------------------------------------------------------
function lr=lrtest(x, n, mu, sigma)

numBins = size(x, 2)

% Add up how many points there are in the histogram

N = sum(n);

%
% For each frequency, we need the probability that a r.v. falls into a
% given bin. We obtain this by integrating the Gaussian d.f. over the
% range of the bin. The centers of the bins are given by x(1) to x(numBins),
% so the ranges for the bins are given by
%
% (-Inf, (x(1)+x(2))/2] ((x(1)+x(2))/2, (x(2)+x(3))/2] . . .
%        . . . ((x(k-1) + x(k))/2, +Inf)
%
% noting that the two end bins need special treatment because they extend to
% infinity.
%
% The integral is found by integrating the Gaussian d.f. with zero mean and
% unit s.d., and then scaling appropriately, using the Matlab error function
% erf().
%

% First calculate the upper and lower ends of the bins, and translate them
% to zero mean and unit s.d.

a = zeros(1, numBins-1);
for i = 1:numBins-1
    a(i) = (0.5*(x(i) + x(i+1)) - mu)/(sigma*sqrt(2));
end;

% Now calculate the areas (probabilities) for interior bins

prob = zeros(1, numBins);

for i = 2:numBins-1
    prob(i) = 0.5*(erf(a(i)) - erf(a(i-1)));
end;

% Now the last bin from a(end) to +Inf

prob(end) = 0.5*erfc(a(end));

% Now the first bin from a(1) to -Inf (actually, from -a(1) to +Inf)

prob(1) = 0.5*erfc(-a(1));

mur=sum(x.*n)/N
sdr=sqrt(sum((x-mur).^2.*n)/N)
Nprob=gausshist(mu,sigma,N,x);
prob-Nprob/N
Nprob=round(Nprob);
toplot=n;
%toplot=Nprob;
%toplot=round(Nprob);
doublehist(x,toplot,Nprob)

% Now calculate the likelihood ratio

% One method, just take the product, doesn't work well when N is large
%
%lr = 1;
%
%for i = 1:numBins
%    if (n(i) > 0)
%        lr = lr*(N*prob(i)/n(i))^n(i);
%    end;
%end;

% Alternative lnlr formula
%
%lnlr   = N*log(N)   + sum(n.*log(prob)) ...
%             - sum(n.*log(n));

lrcomp = zeros(1, numBins);

for i = 1:numBins
    if (n(i) ~= 0)
        lrcomp(i) = n(i)*log(N*prob(i)/n(i));
    end;
end;

lnlr = sum(lrcomp);

lr = exp(lnlr);


lnlr=gammaln(sum(toplot))-gammaln(N);
for i = 1:numBins
    if (prob(i) > 0)
        term1=gammaln((toplot(i))+1)-gammaln(N*prob(i)+1);
        term2=((toplot(i))-N*prob(i))*log(prob(i));
        lnlr = lnlr-term1+term2;    
    end;
end;

lr =exp(lnlr/N);