function PoissHTsim(A0,A1,p1,T,Nbins);
%function PoissHTsim(A0,A1,p1,T,Nbins);
%Demonstrates the testing of two Poiss Hypth
% H0:Poiss(A0), H1:Poiss(A1) with prob p1
% T iid samples for Nbins histogram and testing
% A typical call is
% PoissHTsim(A0,A1,p1,T,Nbins);

lgA0 = log(A0);lgA1 = log(A1);
Kmax = round( 2*max(A0,A1) );% A choice of plotting range
K = [1:Kmax];lgK = log(K);%Row
%These could be computed simultaneously ...

pmf0 = exp( K*lgA0 -A0 -cumsum( lgK ));
sumpmf0 = sum(pmf0);pmf0=pmf0/sumpmf0;%This is an approx!
pmf1 = exp( K*lgA1 -A1 -cumsum( lgK ));
sumpmf1 = sum(pmf1);pmf1 = pmf1/sumpmf1;%This is an approx!

K = K(:);PMF = [ pmf0; pmf1]';%Convert to Column Format
% Let's compute the mean values from theory
MN = K'*PMF;% Sample means computed from dataset

% Hypotesis Testing Quantities
% thb is derived from the Likelihood ratio test with Bayes Threshold.
% Thresh_bayes depnds on Poisson parameters A0,A1 and p1=1-p0, the prior prob
p0 = 1 -p1;eta = log( p0/p1); beta= (A1/A0);
thb = (eta+A0*(beta-1) )/log(beta);
thb = round(thb);
disp(' ')
disp('Threshold and Estimated Mean values')
disp(' ')

disp(' Thb Mn0 Mn1 ')
disp( [ thb MN ] )
disp(' ')

figure(1)
K = K(:);PMF = [ pmf0; pmf1]';%Convert to Column Format
plot( K, PMF, '*');grid
xlabel('K Counts');ylabel('PMF');
title('Theory: Poisson PMF under each Hm');
PARaxis = axis;
hold on
h1= plot( [thb thb], [0 PARaxis(4)],'r' );
% Draw the threshold line in red
set(h1,'LineWidth',2);
hold off

%Compute Error Rates Qth11 Qth01
%Typically Qth11 = Qth01 only when p1=p0 1/2
Qth11 = sum(pmf1(find(K >= thb)));%Notethe >= choice
Qth01 = sum(pmf0(find(K >= thb)));% These are computed from the PMF

%Qth01th=gammainc(A0,thb);% alternative computation in terms of the incomplete gamma
%Qth11th=gammainc(A1,thb);% These check
disp(' ')
disp('Qthij := Prob( Poiss(H1) > Thb');
disp(' ')
disp(' Qth11 Qth01 ')
disp( [ Qth11 Qth01 ] )
disp(' ')

%Now lets run a numerical simuation experiment
% based on Monte Carlo with T IID trials

% Generate some random data according to the Hypoth
AlldataK = poissrnd(repmat([A0 A1],T,1));%does both hypoth
%Column 1 is for H0, Column 2 for H1

%our data sample whould have a proportion p1 of Poiss(A1) and p0 of Poiss(A0)
%Get a random proportion
H1indx = (rand(T,1) < p1);% Indicates Ho or H1 - Bern(p1)
Lh1 = sum(find(H1indx==1));Lh0 = T-Lh1;% Lengths

%Sindx = zeros(T,1);
%whos
%Sindx(T1rnd) = ones( LT1,1);

%T1rnd is an indicator to choose whether sample t of t=1:T is drawn from H1
%LT0 = length(T0rnd);

dataK0 = AlldataK([1:T0rnd], 1);% A random selection from H0
dataK1 = AlldataK([1:T1rnd], 2);% A random selection from H1
% We can now randomize the order

[sx,indx ]= sort(rand(T,1));
%dataK = dataK(indx);

figure(2);
plot( T1rnd,'r*');grid

figure(1)
break

BinC = Kmax*[1:Nbins]/Nbins;
[Nknt1] = hist(dataK1,BinC);
[Nknt0] = hist(dataK0,BinC);

figure(3)
plot( BinC,Nknt1,'b*',BinC,Nknt0,'g*');grid;
xlabel('K Counts');ylabel('PMF');
title('Experiment: Poisson Histogram under each Hm');

PARaxis = axis;
hold on
h2= plot( [thb thb], [0 PARaxis(4)],'r' );
set(h2,'LineWidth',2);
hold off

%hold on
%plot( ones(T1rnd,1),dataK0,'g*');grid
%hold off

figure(1)