%-----------------------------------------------------------------------------
%
% $Id: testlr.m,v 1.4 1999/09/22 03:29:26 philipc Exp $
%
% Program to test the likelihood ratio code. Requires the Statistics Toolbox.
% Parameters specify which test case to run, and number
% of generated points to use when needed.
%
%-----------------------------------------------------------------------------

function testlr(testcase, numPoints)

if (~exist('testcase'))
    testcase = 1;
end;

if (~exist('numPoints'))
    numPoints = 100;
end;

if (testcase == 1 | testcase == 0)
%
fprintf('Test case %d: ', testcase);
fprintf('Gaussian distribution\n\n');
%

mu = 0.0;
sigma = 1.0;

y = normrnd(mu, sigma, 1, numPoints);
mur=mean(y)
sdr=std(y)
x = [-25:25]/5;
n = hist(y, x);

lr = lrtest(x, n, mu, sigma);
%lr = lrtest(x, n, mur, sdr);

fprintf('Likelihood-ratio = %f\n\n', lr);

return;
end;
%
% End of test case 1
%
%----------------------------------------------------------------------------

if (testcase == 2 | testcase == 0)
%
fprintf('Test case %d: ', testcase);
fprintf('Gaussian distribution\n\n');
%

mu = 5.0;
sigma = 3.0;

y = normrnd(mu, sigma, 1, numPoints);

x = [-12.5:37.5]/2.5;
n = hist(y, x);

lr = lrtest(x, n, mu, sigma);

fprintf('Likelihood-ratio = %f\n\n', lr);

return;
end;
%
% End of test case 1
%
%----------------------------------------------------------------------------

if (testcase == 3 | testcase == 0)
%
fprintf('Test case %d: ', testcase);
fprintf('Real data\n\n');
%

load results numFreq yy numBins nnr nni;

% Add up how many points there are in each histogram, should be the same
% for all histograms

N = sum(nnr(1,:));

% Calculate mean mr(nf) for each frequency nf.

mr = sum(yy.*nnr, 2)/N;

% Calculate the variance and standard deviation sdr(nf) for each frequency nf.
yy_new = yy - mr*ones(1,numBins);

vr = sum(yy_new.*yy_new.*nnr,2)/(N-1);
sdr = sqrt(vr);

nf = 8;

fprintf('mu = %f, sigma = %f\n',  mr(nf), sdr(nf));

lr = lrtest(yy(nf,:), nnr(nf,:), mr(nf), sdr(nf));

fprintf('Likelihood-ratio = %f\n\n', lr);

return;
end;
%
% End of test case 3
%
%----------------------------------------------------------------------------

if (testcase == 4 | testcase == 0)
%
fprintf('Test case %d: ', testcase);
fprintf('Perfect Gaussian data\n\n');
%

mu = 0.0;
sigma = 1.0;

x = -5:5;

numBins = size(x,2);

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

prob = zeros(1, numBins);

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

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

n = round(numPoints*prob);

lr = lrtest(x, n, mu, sigma);

fprintf('Likelihood-ratio = %f\n\n', lr);

return;
end;
%
% End of test case 4
%
%----------------------------------------------------------------------------