% Written by Min-Hua Huang
% June 9, 2011
% Testing ten cumulative distribution functions related to the Gamma function with the "h" factorization method.
% Section 4.2 in "Indefinite Integration of the Gamma Integral".
% c is an considerably larger number, as a proxy of the infinite upper limit. 
clear; clc; format long; n=10^5;
x=1.5; c=20; a=2; b=1.5;
fprintf('Testing ten distribution functions associated with the Gamma function:\n\n');
% ************************************************************************
% Gamma distribution
answer1=gamcdf(x,a,b); 
check1=(x/(b*c))^a*exp(c-x/b)*h(a-1,-x/b,n)/h(a-1,-c,n);
fprintf('Gamma=%8.5f, h method=%8.5f,\n',answer1,check1);
% ************************************************************************
% Poisson distribution
k=4; lamda=2.5;
answer4=poisscdf(k,lamda);
check4=1-(lamda/c)^(k+1)*exp(c-lamda)*h(k,-lamda,n)/h(k,-c,n); 
fprintf('Poisson=%8.5f, h method=%8.5f,\n',answer4,check4);
% ************************************************************************
% Chi-square distribution
k=3; 
answer7=chi2cdf(x,k);
check7=(x/(2*c))^(k/2)*exp(c-x/2)*h(k/2-1,-x/2,n)/h(k/2-1,-c,n); 
fprintf('Chi-square=%8.5f, h method=%8.5f,\n',answer7,check7);
% ************************************************************************
% Erlang distribution
answer2=gammainc(x*b,a);
check2=(b*x/c)^a*exp(c-b*x)*h(a-1,-1*b*x,n)/h(a-1,-c,n);
fprintf('Erlang=%8.5f, h method=%8.5f,\n',answer2,check2);
% ************************************************************************
% Inverse-Gamma distribution
answer3=1-gammainc(b/x,a);
check3=1-(b/(c*x))^a*exp(c-b/x)*h(a-1,-b/x,n)/h(a-1,-c,n);
fprintf('Invse-gamma=%8.5f, h method=%8.5f,\n',answer3,check3);
% ************************************************************************
% Chi distribution
k=2.5; 
answer5=gammainc(x^2/2,k/2);
check5=(x/sqrt(2*c))^(k)*exp(c-x^2/2)*h(k/2-1,-x^2/2,n)/h(k/2-1,-c,n); 
fprintf('Chi=%8.5f, h method=%8.5f,\n',answer5,check5);
% ************************************************************************
% Noncentral Chi distribution
k=1; lamda=1; 
answer6=NaN;
check6=exp(-lamda^2/2)/(2^(k/2)*c^(k/2)*exp(-c)*h(k/2-1,-c,n));
sum6=x^(k/2)*exp(-x/2)*h(k/2-1,-x/2,n); product=1;
for i=1:1:150
    product=product*(lamda^2/4)/(i*(k/2-1+i)); 
    sum6=sum6+x^(k/2+i)*exp(-x/2)*h(k/2-1+i,-x/2,n)*product;
end
check6=check6*sum6;
fprintf('Noncentral-chi="not available", h method=%8.5f,\n',check6);
% ************************************************************************
% Inverse-chi-square distribution
v=3; 
answer8=1-gammainc(1/(2*x),v/2);
check8=1-(1/(2*c*x))^(v/2)*exp(c-1/(2*x))*h(v/2-1,-1/(2*x),n)/h(v/2-1,-c,n); 
fprintf('Inverse-chi-square=%8.5f, h method=%8.5f,\n',answer8,check8);
% ************************************************************************
% Scaled-inverse-chi-square distribution
v=3; sigma2=2;
answer9=1-gammainc(sigma2*v/(2*x),v/2);
check9=1-(sigma2*v/(2*c*x))^(v/2)*exp(c-sigma2*v/(2*x))*h(v/2-1,-sigma2*v/(2*x),n)/h(v/2-1,-c,n); 
fprintf('Scaled-inverse-chi-square=%8.5f, h method=%8.5f,\n',answer9,check9);
% ************************************************************************
% Generalized normal distribution
mu=0;
answer10=1/2+sign(x-mu)*gammainc((abs(x-mu)/a)^b,1/b)/2;
check10=(1+sign(x-mu)*abs(x-mu)/(a*c^(1/b))*exp(c-(abs(x-mu)/a)^b)*h(1/b-1,-(abs(x-mu)/a)^b,n)/h(1/b-1,-c,n))/2; 
fprintf('Generalized normal=%8.5f, h method=%8.5f.\n\n',answer10,check10);
