% Written by Min-Hua Huang
% June 9, 2011
% Testing eight cumulative distribution functions related to the error function with the "h" factorization method.
% Section 4.4 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=3.5; mu=2; sigma=1.7;
fprintf('Testing eight distribution functions associated with the Error function:\n\n');
% ************************************************************************
% Normal distribution
answer1=normcdf(x,mu,sigma); 
check1=(1+(x-mu)/(sqrt(2*pi)*sigma)*exp(-(x-mu)^2/(2*sigma^2))*h(-1/2,-(x-mu)^2/(2*sigma^2),n))/2;
fprintf('Normal=%8.5f, h method=%8.5f,\n',answer1,check1);
% ************************************************************************
% Inverse Gaussian distribution
lamda=1;
answer2=normcdf(sqrt(lamda/x)*(x/mu-1),0,1)+exp(2*lamda/mu)*normcdf(-sqrt(lamda/x)*(x/mu+1),0,1);
check21=1/2*(1+sqrt(lamda/(2*pi*x))*(x/mu-1)*exp(-lamda/(2*x)*(x/mu-1)^2)*h(-1/2,-(x/mu-1)^2*lamda/(2*x),n));
check22=1/2*exp(2*lamda/mu)*(1-sqrt(lamda/(2*pi*x))*(x/mu+1)*exp(-lamda/(2*x)*(x/mu+1)^2)*h(-1/2,-(x/mu+1)^2*lamda/(2*x),n));
check2=check21+check22;
fprintf('Inverse Gaussian=%8.5f, h method=%8.5f,\n',answer2,check2);
% ************************************************************************
% Log-normal distribution
answer3=1/2+1/2*erf((log(x)-mu)/sqrt(2*sigma^2));
check3=1/2+(log(x)-mu)/(2*sqrt(2*pi)*sigma)*exp(-(log(x)-mu)^2/(2*sigma^2))*h(-1/2,-(log(x)-mu)^2/(2*sigma^2),n);
fprintf('Log-normal=%8.5f, h method=%8.5f,\n',answer3,check3);
% ************************************************************************
% Logit-normal distribution
p=0.3; 
logit=log(p/(1-p));
answer4=1/2*(1+erf((logit-mu)/sqrt(2*sigma^2)));
check4=1/2*(1+(logit-mu)/(sqrt(2*pi)*sigma)*exp(-(logit-mu)^2/(2*sigma^2))*h(-1/2,-(logit-mu)^2/(2*sigma^2),n));
fprintf('Logit-normal=%8.5f, h method=%8.5f,\n',answer4,check4);
% ************************************************************************
% Half-normal distribution
answer5=erf(x/(sqrt(2)*sigma));
check5=x/(sqrt(2*pi)*sigma)*exp(-x^2/(2*sigma^2))*h(-1/2,-x^2/(2*sigma^2),n);
fprintf('Half-normal=%8.5f, h method=%8.5f,\n',answer5,check5);
% ************************************************************************
% Folded normal distribution
answer6=1/2*(erf((x+mu)/(sqrt(2)*sigma))+erf((x-mu)/(sqrt(2)*sigma)));
check61=1/2*((x+mu)/(sqrt(2*pi)*sigma)*exp(-(x+mu)^2/(2*sigma^2))*h(-1/2,-(x+mu)^2/(2*sigma^2),n));
check62=1/2*((x-mu)/(sqrt(2*pi)*sigma)*exp(-(x-mu)^2/(2*sigma^2))*h(-1/2,-(x-mu)^2/(2*sigma^2),n));
check6=check61+check62;
fprintf('Folded normal=%8.5f, h method=%8.5f,\n',answer6,check6);
% ************************************************************************
% Maxwell-Boltzmann distribution
a=2.5;
answer7=erf(x/(sqrt(2)*a))-sqrt(2/pi)*x*exp(-x^2/(2*a^2))/a;
check7=x/(sqrt(2*pi)*a)*exp(-x^2/(2*a^2))*h(-1/2,-x^2/(2*a^2),n)-sqrt(2/pi)*x*exp(-x^2/(2*a^2))/a;
fprintf('Maxwell-Boltzmann=%8.5f, h method=%8.5f,\n',answer7,check7);
% ************************************************************************
% Levy distribution
c=3;
answer8=erfc(sqrt(c/(2*(x-mu))));
check8=1-sqrt(c/(2*pi*(x-mu)))*exp(-c/(2*(x-mu)))*h(-1/2,-c/(2*(x-mu)),n);
fprintf('Levy=%8.5f, h method=%8.5f.\n\n',answer8,check8);

