vb版GammaDist()和GammaInv()两函数的代码

rt
这两个函数excel里有，也可以在vb里调用，但运算速度有点慢
哪位有vb版的函数代码能不能共享下， no23q@163.com，不盛感激。

自己从matlab里翻译过来，如果有兴趣请大家测试

GammaDist(x,a,b,true):

'%GAMCDF Gamma cumulative distribution function.
'% P = GAMCDF(X,A,B) returns the gamma cumulative distribution
'% function with parameters A and B at the values in X.

Public Function GamCdf(x, a, b) '对应excel中GammaDist(x,a,b,true)函数,计算S曲线
If a <= 0 Or b <= 0 Then GammCdf = "NaN"
GamCdf = GAMMAINC(x / b, a)
p = IIf(p > 1, 1, p)
End Function

'%GAMMAINC Incomplete gamma function.
'% Y = GAMMAINC(X,A) evaluates the incomplete gamma function for
'% corresponding elements of X and A. X and A must be real and the same
'% size (or either can be a scalar). A must also be non-negative.
'% The incomplete gamma function is defined as:
'% http://zanjero.ygblog.com/
'% gammainc(x,a) = 1 ./ gamma(a) .*
'% integral from 0 to x of t^(a-1) exp(-t) dt
'%
'% For any a>=0, as x approaches infinity, gammainc(x,a) approaches 1.
'% For small x and a, gammainc(x,a) ~= x^a, so gammainc(0,0) = 1.
Public Function GAMMAINC(x, a)
Dim amax As Double
amax = 2 ^ 20

ascalar = 1
If a <= amax Then
If a <> 0 And x <> 0 And x < a + 1 Then
xk = x
ak = a
ap = ak
Sum = 1 / ap
del = Sum
Dim i As Double
For i = 1 To 10000
ap = ap + 1
del = xk * del / ap
Sum = Sum + del
Next
GAMMAINC = Sum * Exp(-xk + ak * Log(xk) - GAMMALN(ak))

ElseIf a <> 0 And x <> 0 And x >= a + 1 Then

xk = x
a0 = 1
a1 = x
b0 = 0
b1 = a0
ak = a
fac = 1
n = 1
g = b1
gold = b0

For i = 1 To 10000
gold = g
ana = n - ak
a0 = (a1 + a0 * ana) * fac
b0 = (b1 + b0 * ana) * fac
anf = n * fac
a1 = xk * a0 + anf * a1
b1 = xk * b0 + anf * b1
fac = 1 / a1
g = b1 * fac
n = n + 1

Next
GAMMAINC = 1 - Exp(-xk + ak * Log(xk) - GAMMALN(ak)) * g

End If
Else

End If

End Function

'*******************************************************************************************************************
'%GAMMALN Logarithm of gamma function.
'% Y = GAMMALN(X) computes the natural logarithm of the gamma
'% function for each element of X. GAMMALN is defined as
'%
'% LOG(GAMMA(X))
'%
'% and is obtained without computing GAMMA(X). Since the gamma
'% function can range over very large or very small values, its
'% logarithm is sometimes more useful. http://zanjero.ygblog.com/
Public Function GAMMALN(XX)
Dim COF(6) As Double, stp As Double, half As Double, one As Double
Dim fpf As Double, x As Double, tmp As Double, ser As Double
Dim j As Integer
COF(1) = 76.18009173
COF(2) = -86.50532033
COF(3) = 24.01409822
COF(4) = -1.231739516
COF(5) = 0.00120858003
COF(6) = -0.00000536382
stp = 2.50662827465
half = 0.5
one = 1#
fpf = 5.5
x = XX - one
tmp = x + fpf
tmp = (x + half) * Log(tmp) - tmp
ser = one
For j = 1 To 6
x = x + one
ser = ser + COF(j) / x
Next j
GAMMALN = tmp + Log(stp * ser)
End Function

[ 本帖最后由 no23q 于 2008-10-13 10:52 编辑 ]

no23q

2008年10月13日 10:49:10

2楼

'%GAMINV Inverse of the gamma cumulative distribution function (cdf).
'% X = GAMINV(P,A,B) returns the inverse of the gamma cdf with
'% parameters A and B, at the probabilities in P. http://zanjero.ygblog.com/
'对应Excel中GammaInv(p,a,b)函数，计算PⅢ的Φp值 Φp＝Cs/2*GammaInv(1-p/100,4/Cs^2,1)-2/Cs
Public Function gaminv(p, a, b)
If p < 0 Or p > 1 Or a <= 0 Or b <= 0 Then
gaminv = "NaN"
Exit Function

Else
Select Case p
Case 0
gaminv = 0
Case 1
gaminv = 1
Case Else
'% Newton's Method
'% Permit no more than count_limit interations.

count_limit = 100
cunt = 0
pk = p
'% Supply a starting guess for the iteration.
'% Use a method of moments fit to the lognormal distribution.

mn = a * b
v = mn * b
temp = Log(v + mn ^ 2)
mu = 2 * Log(mn) - 0.5 * temp
sigma = -2 * Log(mn) + temp

xk = Exp(MyNormInv(pk, mu, sigma)) '''''''''''''''''''''''

h = 1
'% Break out of the iteration loop for three reasons:
'% 1) the last update is very small (compared to x)
'% 2) the last update is very small (compared to sqrt(eps))
'% 3) There are more than 100 iterations. This should NEVER happen.
eps = 2 ^ -10
Do While Abs(h) > eps ^ 0.5 * Abs(xk) And Abs(h) > eps ^ 0.5 And cunt < count_limit

cunt = tunt + 1
h = (GamCdf(xk, a, b) - pk) / gampdf(xk, a, b) '''''''''''''
xnew = xk - h
' % Make sure that the current guess stays greater than zero.
' % When Newton's Method suggests steps that lead to negative guesses
' % take a step 9/10ths of the way to zero:

If xnew < 0 Then
xnew = xk / 10
h = xk - xnew
End If
xk = xnew

Loop
gaminv = xk

End Select

End If

End Function

' This function is a replacement for the Microsoft Excel Worksheet function NORMSINV.
' It uses the algorithm of Peter J. Acklam to compute the inverse normal cumulative
' distribution. Refer to http://home.online.no/~pjacklam/notes/invnorm/index.html for
' a description of the algorithm.
' Adapted to VB by Christian d'Heureuse, http://zanjero.ygblog.com/.
Public Function MyNormSInv(ByVal p As Double) '计算频率格纸
Const a1 = -39.6968302866538, a2 = 220.946098424521, a3 = -275.928510446969
Const a4 = 138.357751867269, a5 = -30.6647980661472, a6 = 2.50662827745924
Const b1 = -54.4760987982241, b2 = 161.585836858041, b3 = -155.698979859887
Const b4 = 66.8013118877197, b5 = -13.2806815528857, c1 = -7.78489400243029E-03
Const c2 = -0.322396458041136, c3 = -2.40075827716184, c4 = -2.54973253934373
Const c5 = 4.37466414146497, c6 = 2.93816398269878, d1 = 7.78469570904146E-03
Const d2 = 0.32246712907004, d3 = 2.445134137143, d4 = 3.75440866190742
Const p_low = 0.02425, p_high = 1 - p_low
Dim q As Double, r As Double
If p < 0 Or p > 1 Then
Err.Raise vbObjectError, , "NormSInv: Argument out of range."
ElseIf p < p_low Then
q = Sqr(-2 * Log(p))
MyNormSInv = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _
((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
ElseIf p <= p_high Then
q = p - 0.5: r = q * q
MyNormSInv = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / _
(((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1)
Else
q = Sqr(-2 * Log(1 - p))
MyNormSInv = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / _
((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
End If
End Function

'%GAMPDF Gamma probability density function.
'% Y = GAMPDF(X,A,B) returns the gamma probability density function
'% with parameters A and B, at the values in X. http://zanjero.ygblog.com/
Public Function gampdf(x, a, b)
y = 0
If x = 0 And a < 1 Then
gampdf = "∞"
Exit Function
End If

If x = 0 And a = 1 Then
gampdf = 1 / b
Exit Function
End If

If a <= 0 Or b <= 0 Then
y = "NaN"
gampdf = y
Exit Function
ElseIf x > 0 Then

y = (a - 1) * Log(x) - x / b - GAMMALN(a) - a * Log(b)
y = Exp(y)
End If
gampdf = y

End Function
:victory:

[ 本帖最后由 no23q 于 2008-10-13 10:52 编辑 ]

zhphto

2008年10月14日 08:26:47

3楼

你是想写频率曲线程序还是瞬时单位线推洪水，
我看论坛里这个叫freeskychenjun 估计有