% Written by Min-Hua Huang
% June 9, 2011
% Evaluating the incomplete beta function by using the h functions.
% input: hbetainc(x,a,b,nbeta,pre); x is the random variable; a is parameter 
% alpha; b is parameter beta; nbeta is the number of expansions for 
% calcutaing the incomplete beta function;  pre:level of precision (number
% of singificant digits in matlab's vpa function).
function fv=hbetainc(x,a,b,nbeta,pre)
format long; digits(pre);
record=vpa(ones(1,3)); d=vpa(x^a/a); sum=vpa(d);
record(1,1)=vpa(d);
record(1,2)=vpa(sum);
for i=1:1:nbeta-1
    d=vpa(d*(-x)*(a+i-1)/(i*(a+i)));  
    sum=vpa(sum+d);    
    record(i+1,1)=vpa(d);
    record(i+1,2)=vpa(sum);
end
hv=vpa(sum); hanswer=vpa(sum);
for i=1:1:nbeta
product=1;
    for k=1:1:i+1
        if k==i
           product=vpa(product*1);
        else
           product=vpa(product*(b-k));
        end   
    end
hv=vpa(hv-record(i,1));    
hanswer=vpa(hanswer+product*hv);
record(i,3)=vpa(hv);
record(i,4)=vpa(product);
record(i,5)=vpa(hanswer);
end    
fv=hanswer;

