% Written by Min-Hua Huang
% June 9, 2011
% Testing eight cumulative distribution functions related to the beta function with the "h" factorization method.
% Section 4.5 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=0.4; c=20; a=1.6; b=2.3;
fprintf('Testing eight distribution functions associated with the Beta function:\n\n');
% ************************************************************************
% Beta distribution
answer1=betainc(x,a,b); 
check1=double(hbetainc(x,a,b,100,700)*h(a+b-1,-c,n)/(exp(-c)*h(a-1,-c,n)*h(b-1,-c,n)));
fprintf('Beta=%8.5f, h method=%8.5f,\n',answer1,check1);
% ************************************************************************
% Binomial distribution
p=0.67; k=7; N=10; c=50;
answer3=betainc(1-p,N-k,1+k);
check3=double(hbetainc(1-p,N-k,1+k,150,1000)*h1(N,-c,100)/(exp(-c)*h1(N-k-1,-c,100)*h1(k,-c,100)));
fprintf('Binomial=%8.5f, h method=%8.5f,\n',answer3,check3);
% ************************************************************************
% F distribution
d1=2; d2=3; 
answer6=betainc(d1*x/(d1*x+d2),d1/2,d2/2);
check6=double(hbetainc(d1*x/(d1*x+d2),d1/2,d2/2,100,700)*h((d1+d2)/2-1,-c,n)/(exp(-c)*h(d1/2-1,-c,n)*h(d2/2-1,-c,n)));
fprintf('F=%8.5f, h method=%8.5f,\n',answer6,check6);
% ************************************************************************
% Beta prime distribution
answer2=betainc(x/(1+x),a,b);
check2=double(hbetainc(x/(1+x),a,b,100,700)*h(a+b-1,-c,n)/(exp(-c)*h(a-1,-c,n)*h(b-1,-c,n)));
fprintf('Beta prime=%8.5f, h method=%8.5f,\n',answer2,check2);
% ************************************************************************
% Negative binomial distribution
p=0.67; k=5; r=3; c=100;
answer4=1-betainc(p,k+1,r);
check4=double(1-hbetainc(p,k+1,r,150,1000)*h1(k+r,-c,700)/(exp(-c)*h1(k,-c,700)*h1(r-1,-c,700)));
fprintf('Negative binomial=%8.5f, h method=%8.5f,\n',answer4,check4);
% ************************************************************************
% Yule-Simon distribution
k=2; rho=1.5; c=20;
answer5=1-k*beta(k,rho+1);
check5=1-k*exp(-c)*h(k-1,-c,n)*h(rho,-c,n)/h(k+rho,-c,n); 
fprintf('Yule-Simon=%8.5f, h method=%8.5f,\n',answer5,check5);
% ************************************************************************
% Noncentral F distribution
d1=2; d2=3; lamda=0.5; 
answer7=ncfcdf(x,2,3,0.5);
check7=0;
for i=0:1:15
    t=double(hbetainc(d1*x/(d1*x+d2),d1/2+i,d2/2,30,200)*h((d1+d2)/2+i-1,-c,n)/(exp(-c)*h(d1/2+i-1,-c,n)*h(d2/2-1,-c,n)));
    check7=check7+(lamda/2)^i*exp(-lamda/2)/factorial(i)*t; 
end    
fprintf('Noncentral F=%8.5f, h method=%8.5f,\n',answer7,check7);
% ************************************************************************
% Noncentral t-distribution
v=3; lamda=2;
answer8=nctcdf(x,3,2);
if x>=0
   check8=normcdf(-lamda);
   for i=0:1:30 
       t1=(lamda^2/2)^(i/2)/(c^(i/2+1)*exp(-c)*h(i/2,-c,n));
       t2=double(hbetainc(x^2/(v+x^2),(i+1)/2,v/2,50,700)*h((i+v-1)/2,-c,n)/(exp(-c)*h((i-1)/2,-c,n)*h(v/2-1,-c,n)));
       check8=check8+t1*t2*exp(-lamda^2/2)/2;
   end
else
   check8=normcdf(-lamda);
   for i=0:1:20 
       t1=(lamda^2/2)^(i/2)/(c^(i/2+1)*exp(-c)*h(i/2,-c,n));
       t2=double(hbetainc(x^2/(v+x^2),(i+1)/2,v/2,50,700)*h((i+v-1)/2,-c,n)/(exp(-c)*h((i-1)/2,-c,n)*h(v/2-1,-c,n)));
       t3=(lamda^2/2)^i/factorial(i);
       t4=double(hbetainc(x^2/(v+x^2),i+1/2,v/2,50,700)*h(i+(v-1)/2,-c,n)/(exp(-c)*h(i-1/2,-c,n)*h(v/2-1,-c,n)));
       check8=check8+t1*t2*exp(-lamda^2/2)/2-t3*t4*exp(-lamda^2/2);
   end
end
fprintf('Noncentral t=%8.5f, h method=%8.5f.\n\n',answer8,check8);
