Dear Young people of this universe.

I am in an obsessive effort to make a calculation tool in VB from Scrach... At the moment I am trying to do the Gamma function in VB codes.

Using trapezoidal rule, I managed to do approximation for finite integration... just about anything. I can even do the normal curve without too much problem... cross checked with the minitab table, it is right and quite okay to use.

but when I try to do T distribution, I run into the Gamma function. it was a pain. using 1000 strips on trapesoidal rule, the errors was quite obvious. So I made 1 million strips, the accracy improved and seems to get quite good result for T, X2 and F distribution.

However, it was still an estimation which will not work well on larger numbers. as a result, I looked two ways. one is to find a good estimation for the function: at the moment I am looking at Lanczos approximation, I am still trying to do the matrix bit in code.

another approach is to look if the function can actually be integrated... silly me, if it is possible then the person who make the Gamma function would not have the intgral in the first place. Anyhow, I tried to integrate it... using "integration by part", the result turned out 0=0.. not very cool.

so my question is that is there any way that I can integrate the Gamma function (I guess not)? and, is Lanczos the way to go? Microsoft Excel seems to have a very fast equation for it... anyone know how they did it?

my many thanks in advance...

JL

There are several known approximations to the gamma function. Stirling's formula is asymptotic to the gamma function, but is not a particularly good approximation (even for very large numbers, the second significant figure is often incorrect). An excellent approximation can be obtained from the following formula (due to Gosper I believe):

G(x+1) » ((2x+1/3)p)^1/2x^x e^-x
It appears to give around four significant figure accuracy. I haven't heard of the Lanczos approximation so I cannot tell you how good it is.

I hope this helps.

Paul

Thanks for the reply, at the moment trapezoidal rule 1 million strips provide a very very good result on for small number inputs... between 2.2 and 3.2 the error (compare to MicroSoft Excel Exp(GAMMALN())) ranged between 1.55442E-09 to
1.25174E-08. bigger number like 50 has a bigger error term of -1.9599E-06.

I am on a P4 1.7GHz 256 ram (at work, my home PC is more ... humble)it take about 5 sec to calculate.

I am thinking this error is quite okay and may be I can use the recurrence property of Gamma to calculate the rest of the values. however if the number get bigger the error mentioned above will get manified... which might result as un accpetable. that's why I want to get it right... well or very close to right.

saying so, I think T distribution and F and X2 tend to be affected by the errors in the Gamma function when the gamma variable is small. for large variables the existing Gamma is not right but the error seems to be "hidden" in the T F and X2 distributions...

The improvement that I am require of is to, above most, improve the accuracy of the Gamma and secondly if possible improve the run speed.

I find the info on Lanzcos via this website http://www.rskey.org/gamma.htm the dude did it in C++ manage to have 158 S.F. for 100!... very impressive... but it also said it took a while to run...

I just tried the Stirlings approx. since the errors will be manified in the distributions I really dun think it fits for me..

anyway thanks for the reply and may be we can evaluate various approx for Gamma here.

JL

for got to say the errors terms are percentages... and the upper limit of Gamma (the bit that suppose to be infinity) is in fact just 99. for some reason it seems to work fine.

John,

You mention integrating the normal curve ... There is a very good approximation for the error function, i.e.

2/

Ö

p

ò0xexp(-x2) dx

in ''Handbook of mathematical functions'' (Abramowitz and Stegun) of the form
erf(x)=1-(a₁ t+...+a₅ t⁵ )exp(-x² )

where t=1/(1+px)

for constants p, t_i .

The error is < =1.5*10^-7 for x in [0,infinity ]

This is made use of extensively in finance, in computing option prices using the famous/infamous Black Scholes formula (which requires 2 normal curve integrations). Numerical integration of the normal curve is just too slow. Indeed, there are some options ('double knockouts') whose value can be expressed as an infinite sum where each term involves an error function - so you can see why speed is important!

As for the gamma function, Abramowitz & Stegun suggests a few approximations:

i) The Stirling approximation complete with the (1+1/(12x)+1/(288*x² ) etc) factor

ii) A power series expansion of 1/gamma(x)

iii) A polynomial expansion for gamma(x+1) with an error of less than 5*10^-5 for x in [0,1]

Andre