محاسبات تابع گاما

[arashkey]

می خواستم برنامه ای بنویسم که با آن بتوانم نمودار هایی برای یکسری توابع آماری از جمله پواسن و گاما و t استیودنت و ... بنویسم
متاسفانه تابع گاما به سادگی دیگر توابع نبود
و چون یکی از مهمترین توابع آماری به شمار می رود و در توابع دیگری مانند t استیودنت و فیشر و ... کاربرد دارد
حال می خواستم از شما کمک بخواهم تا بتوانم محاسبات بروی این تابع را انجام دهم
یعنی دنبال یک Function هستم که من به آن alfa را بدهم ( یک ورودی) و اون به من معادل محاسبه شده گامای آن را بدهد
در اینترنت زیاد راجع به این تابع تحقیق کردم ولی فقط و فقط یک برنامه ++C از آن گیر آوردم و آن را به VB ترجمه کردم
که می توانید آنرا در پیوست ببینید


' Returns the imcomplete gamma functions P(a,x) and Q(a,x).
' Arash kiani Jan/19/2007
Imports System.Math
Public Class clsGamma
Dim cerr As String
Sub New()
End Sub
Function gammp(ByVal a As Double, ByVal x As Double) As Double
Dim gamser, gammcf, gln As Double
If x < 0.0 Or a <= 0.0 Then
cerr = "Invalid argument in routine gammp"
End If
If x < (a + 1.0) Then
Call gser(gamser, a, x, gln)
Return gamser
Else
gcf(gammcf, a, x, gln)
Return 1.0 - gammcf
End If
End Function
Function gammq(ByVal a As Double, ByVal x As Double) As Double
Dim gamser, gammcf, gln As Double
If x < 0.0 Or a <= 0.0 Then
cerr = "Invalid argument in routine gammp"
Exit Function
End If
If x < a + 1.0 Then
gser(gamser, a, x, gln)
Return 1.0 - gamser
Else
gcf(gammcf, a, x, gln)
Return gammcf
End If
End Function
Const ITMAX = 100
Const EPS = 0.0000003
Sub gser(ByRef gamser As Double, ByVal a As Double, ByVal x As Double, ByRef gln As Double)
Dim n As Integer
Dim sum, del, ap As Double
gln = gamma_ln(a)
If x <= 0.0 Then
If x < 0.0 Then
cerr = "x less than 0 in routine gser"
gamser = 0.0
Exit Sub
End If
Else
ap = a
sum = 1.0 / a
del = sum
n = 1
While n <= ITMAX
ap += 1
del *= x / ap
sum += del
If Abs(del) < Abs(sum) * EPS Then
gamser = sum * Exp(-x + a * Log(x) - gln)
Exit Sub
End If
n += 1
End While
cerr = "a too large,ITMAX too small in routine gser"
Exit Sub
End If
End Sub
Const FPMIN = 1.0E-30
Sub gcf(ByRef gammcf As Double, ByVal a As Double, ByVal x As Double, ByRef gln As Double)
Dim i As Integer
Dim an, b, c, d, del, h As Double
gln = gamma_ln(a)
b = x + 1.0 - a
c = 1.0 / FPMIN
d = 1.0 / b
h = d
i = 1
While i <= ITMAX
an = -i * (i - a)
b += 2.0
d = an * d + b
If Abs(d) < FPMIN Then
d = FPMIN
End If
c = b + an / c
If Abs(c) < FPMIN Then
c = FPMIN
End If
d = 1.0 / d
del = d * c
h *= del
If Abs(del - 1.0) < EPS Then
Exit While
End If
i += 1
End While
If i > ITMAX Then
cerr = "a too large ,ITMAX too small in gcf"
End If
gammcf = Exp(-x + a * Log(x) - gln) * h
End Sub
Function gamma_ln(ByVal xx As Double) As Double
Dim x, y, tmp, ser As Double
Dim cof() As Double = {76.180091729471457, -86.505320329416776, _
24.014098240830911, -1.231739572450155, _
0.001208650973866179, -0.000005395239384953}
Dim j As Integer
x = xx
y = x
tmp = x + 5.5
tmp -= (x + 0.5) * Log(tmp)
ser = 1.0000000001900149
j = 0
While j <= 5
y += 1
ser += cof(j) / y
j += 1
End While
Return -tmp + Log(2.5066282746310007 * ser / x)
End Function
End Class

این تابع دو ورودی x و a می گیرد و من نمی دانم چگونه کاری کنم که فقط با دادن یک alfa معادل Gamma آن را به من بدهد. ( تابع گاما حاصل یک انتگرال است توضیحات بیشتر در این رابطه را در
http://mathworld.wolfram.com/GammaFunction.html
ببیند)
با تشکر از دوستان خواهشمندم هر چه زودتر به من کمک کنید
در ضمن ورودی من صحیح یا ضریبی از عدد 1/2 نیست هر عددی می تواند باشد

[arashkey]

این هم معادل ++C این تابع
خواهشا به من کمک کنید


/************************************************** *******************
Returns the imcomplete gamma functions P(a,x) and Q(a,x).
C.A. Bertulani May/15/2000
************************************************** *******************/

typedef double Number;

#include <math.h>
#include<iostream.h>
#include <conio.h>


int main()
{
clrscr();
Number gammp(Number , Number );
Number gammq(Number , Number );
Number a, x;
int j;
j=1;
cout << "\n Enter a and x. To stop, enter a<0\n";
cout << "\n Wanna check? Note that P(a,0) = 0, and P(a, infnity) = 1 \n";
cin >> a >> x;
while(a>=0){
if(j>10){ cout << "\n My patience is over. Stop, please!\n";
break;
}
if(j!=1){
cout << "\n\n Enter a and x. To stop, enter x<0.\n";
cin >> a >> x;
}
cout << "\n P(a,x): " << gammp(a,x);
cout << "\n Q(a,x): " << gammq(a,x);
j=j+1;
}
return 0;
}
/************************************************** *******************
Returns the imcomplete gamma function
P(a,x) = (int_0^x e^{-t} t^{a-1} dt)/Gamma(a) , (a > 0).
C.A. Bertulani May/15/2000
************************************************** *******************/
Number gammp(Number a, Number x)
{
void gcf(Number *gammcf, Number a, Number x, Number *gln);
void gser(Number *gamser, Number a, Number x, Number *gln);
Number gamser,gammcf,gln;

if (x < 0.0 || a <= 0.0) cerr<< "Invalid arguments in routine gammp";
if (x < (a+1.0)) {
gser(&gamser,a,x,&gln);
return gamser;
} else { /* Use the continued fraction representation */
gcf(&gammcf,a,x,&gln); /* and take its complement. */
return 1.0-gammcf;
}
}
/************************************************** *******************
Returns the imcomplete gamma function
Q(a,x) = 1-P(a,x)
= (int_x^infinity e^{-t} t^{a-1} dt)/Gamma(a) , (a > 0).
C.A. Bertulani May/15/2000
************************************************** *******************/
Number gammq(Number a, Number x)
{
void gcf(Number *gammcf, Number a, Number x, Number *gln);
void gser(Number *gamser, Number a, Number x, Number *gln);
Number gamser,gammcf,gln;

if (x < 0.0 || a <= 0.0) cerr << "Invalid arguments in routine gammq";
if (x < (a+1.0)) { /* Use the series representation */
gser(&gamser,a,x,&gln);
return 1.0-gamser; /* and take its complement. */
} else { /* Use the continued fraction representation. */
gcf(&gammcf,a,x,&gln);
return gammcf;
}
}
/************************************************** *******************
Returns the imcomplete gamma function P(a,x) evaluated by its series
representation as gamser.
Also returns ln(Gamma(a)) as gln.
C.A. Bertulani May/15/2000
************************************************** *******************/
#define ITMAX 100
#define EPS 3.0e-7

void gser(Number *gamser, Number a, Number x, Number *gln)
{
Number gamma_ln(Number xx);
int n;
Number sum,del,ap;

*gln=gamma_ln(a);
if (x <= 0.0) {
if (x < 0.0) cerr << "x less than 0 in routine gser";
*gamser=0.0;
return;
} else {
ap=a;
del=sum=1.0/a;
for (n=1;n<=ITMAX;n++) {
++ap;
del *= x/ap;
sum += del;
if (fabs(del) < fabs(sum)*EPS) {
*gamser=sum*exp(-x+a*log(x)-(*gln));
return;
}
}
cerr << "a too large, ITMAX too small in routine gser";
return;
}
}
#undef ITMAX
#undef EPS
/************************************************** *******************
Returns the imcomplete gamma function Q(a,x) evaluated by its
continued fraction representation as gammcf.
Also returns ln(Gamma(a)) as gln.
C.A. Bertulani May/15/2000
************************************************** *******************/
#define ITMAX 100 /* Maximum allowed number of iterations. */
#define EPS 3.0e-7 /* Relative accuracy */
#define FPMIN 1.0e-30 /* Number near the smallest representable */
/* floating point number. */

void gcf(Number *gammcf, Number a, Number x, Number *gln)
{
Number gamma_ln(Number xx);
int i;
Number an,b,c,d,del,h;

*gln=gamma_ln(a);
b=x+1.0-a; /* etup fr evaluating continued fracion by modified Lent'z */
c=1.0/FPMIN; /* method with b_0 = 0. */
d=1.0/b;
h=d;
for (i=1;i<=ITMAX;i++) { /* Iterate to convergence. */
an = -i*(i-a);
b += 2.0;
d=an*d+b;
if (fabs(d) < FPMIN) d=FPMIN;
c=b+an/c;
if (fabs(c) < FPMIN) c=FPMIN;
d=1.0/d;
del=d*c;
h *= del;
if (fabs(del-1.0) < EPS) break;
}
if (i > ITMAX) cerr << "a too large, ITMAX too small in gcf";
*gammcf=exp(-x+a*log(x)-(*gln))*h; /* Put factors in front. */
}
#undef ITMAX
#undef EPS
#undef FPMIN
/************************************************** ******************
Returns the value of ln[Gamma(xx)] for xx > 0
************************************************** ******************/

Number gamma_ln(Number xx)
{
Number x,y,tmp,ser;
static Number cof[6]={76.18009172947146,-86.50532032941677,
24.01409824083091,-1.231739572450155,
0.1208650973866179e-2,-0.5395239384953e-5};
int j;

y=x=xx;
tmp=x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<=5;j++) ser += cof[j]/++y;
return -tmp+log(2.5066282746310005*ser/x);
}
/************************************************** *******************/

[arashkey]

می دونم سوال ریاضی هستش ولی خوب من جایی دیگه رو نمی شناختم توش این سوال رو بپرسم ( این فرم تلار ریاضی که نداره ) حالا اگه یکی فرومی می شناسه که توش در مورد مسائل ریاضی هم بحث بشه و من بتونم سوالم رو اونجا مطرح کنم لطفا معرفی کنید

[mehdi.hh]

سلام دوست عزیز اولین کدی که گذاشتی یه کلاس است به اسم clsGamma اگر بتوانی آن را به صورت فایل dll ذخیره کنی می‌توانی از توابع‌اش استفاده بکنی اگر درست طراحی شده باشد.
منطق توابع ساده‌اس اگر با خود گاما آشنایی داشته باشی می‌توانی از توابع کلاس استفاده بکنی

[arashkey]

ممنون که حداقل یک نفر جواب داد
راستش این کلاس رو خودم نوشتم
نوشتن که نه از روی همون کد ++C ترجمه کردم
مشکل من اینه که این کدی که من ترجمه کردم دوتا ورودی داره a و x
ولی تابع گاما یک تابعی است که یک رورودی می گیرد و یک خروجی می دهد ( این خروجی حاصل یک انتگرال است)
می خواهم وقتی به آن عدد می دهم ( یعنی یک alfa)به آن بدهم تابع یک مقدار به من بدهد که با آن بتوانم یک نمودار شبیه آنچه در لینک بالا می بینید بکشم

[arashkey]

لطفا یکجا را لینک بدهید تا سوال را در آنجا مطرح کنم
اینگلیسی یا فارسی فرق نمی کند