vb版GammaDist()和GammaInv()两函数的代码.doc
vb版GammaDist()和GammaInv()两函数的代码从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 nonnegative.% The incomplete gamma function is defined as:'http:/zanjero。% 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） = xa, 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 IfEnd 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 DoubleDim fpf As Double， x As Double， tmp As Double， ser As DoubleDim 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 FunctionGammaInv：(要调用上面的GamCdf和GAMMALN函数)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值 pCs/2GammaInv(1-p/100，4/Cs2，1)2/CsPublic 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 Newtons 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 IfEnd 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 dHeureuse, http：/zanjero。ygblog。com/.Public Function MyNormSInv(ByVal p As Double） '计算频率格纸Const a1 = 39.6968302866538, a2 = 220.946098424521， a3 = 275.928510446969Const a4 = 138.357751867269， a5 = 30.6647980661472， a6 = 2。50662827745924Const b1 = -54.4760987982241, b2 = 161.585836858041， b3 = 155。698979859887Const b4 = 66。8013118877197， b5 = -13.2806815528857, c1 = -7。78489400243029E-03Const c2 = 0.322396458041136, c3 = 2。40075827716184， c4 = 2.54973253934373Const c5 = 4。37466414146497， c6 = 2.93816398269878, d1 = 7.78469570904146E03Const d2 = 0.32246712907004, d3 = 2.445134137143, d4 = 3。75440866190742Const p_low = 0.02425， p_high = 1 p_lowDim 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 IfEnd FunctionGAMPDF 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 = yEnd Function
