% Funkce pro testovani Poissonova rozlozeni
%
% Vstupni parametry: Y ... vektor pozorovanych cetnosti
%                    alfa ... hladina vyznamnosti testu (defaultne=0.05)
%                    epsilon ... iteracni rozdil (defaultne=10^(-9))
% Vystupni parametry: Y ... 'modifikovany' vektor pozorovanych cetnosti
%                     p ... vektor empirickych cetnosti
%                     np ... vektor teoretickych cetnosti
%                     ts ... testova statistika
%                     w ... (1-alfa) kvantil Pearsonova rozlozeni
%                           s poctem stupnu volnosti lenght(Y)-2
%                     pval ... p-hodnota pro tento test
%                     iter ... pocet iteraci
%                     lambda ... odhadnuty lambda parametr
%                     H ... vysledek nulove hypotezy (0=nezamitame, 1=zamitame)
% Syntaxe: [Y,p,np,ts,w,pval,iter,lambda,H]=poiss_test(Y,alfa,epsilon)
%
function [Y,p,np,ts,w,pval,iter,lambda,H]=poiss_test(Y,alfa,epsilon);

% defaultni alfa
if (nargin<2) alfa=0.05; end

% defaultni epsilon
if (nargin<3) epsilon=10^(-9); end

% rozsah datoveho souboru
n=sum(Y);

% vazeny prumer
lenY=length(Y);
countY=0:1:(lenY-1);
m=(1/n)*Y*countY';

% vypocet lambda
lambda=m;
diffLambda=inf;
iter=0;
while (diffLambda>=epsilon)
  % vektor p
  p=poisspdf(countY,lambda);
  p(lenY)=1-sum(p)+p(lenY);

  % kontrola podminek npi>5
  np=n*p;
  [min_np,ind]=min(np);
  while ((min_np<=5)&&(lenY>1))
    if (ind==lenY)
      sumY=sum(Y(lenY-1:lenY));
      Y(lenY-1:lenY)=[];
      Y(lenY-1)=sumY;
    else
      sumY=sum(Y(ind:lenY));
      Y(ind:lenY)=[];
      Y(ind)=sumY;
    end  
    lenY=length(Y);   
    countY=0:1:(lenY-1);
    p=poisspdf(countY,lambda);
    p(lenY)=1-sum(p)+p(lenY);
    np=n*p;
    [min_np,ind]=min(np);
  end
  if (min_np<=5)
    warning='Vektor pozorovanych cetnosti nikdy nemuze splnit podminky dobre aproximace, tudiz nasledujici vysledky mohou byt zavadejici'
  end
  
  % citatel podilu dvou sum
  countYY=(lenY-1):1:(1000+lenY-2);
  pp=poisspdf(countYY,lambda);
  numerator=pp*countYY';

  % nova lambda
  tempY=Y;
  tempY(lenY)=0;
  tempY(1)=0;
  newLambda=(1/n)*(tempY*countY'+Y(lenY)*(numerator/p(lenY)));
  
  % rozdil dvou po sobe jdoucich lambd
  diffLambda=abs(lambda-newLambda);
  
  % iterace
  iter=iter+1;
  lambda=newLambda;
end

% vektor teoreticke cetnosti np
np=n*p;

% testova statistika
ts=sum(((Y-np).^2)./np);

% dolni hranice kritickeho oboru
w=chi2inv(1-alfa,lenY-2);

% vysledek hypotezy
if (ts<w) H=0; else H=1; end

% p-value
pval=min(chi2cdf(ts,lenY-2),1-chi2cdf(ts,lenY-2));