java调用Excel中的GAMMADIST函数返回伽玛分布

java调用Excel中的GAMMADIST函数返回伽玛分布
一、问题描述:想用java调用Excel中的GAMMADIST 函数。
Microsoft office Excel中用法说明如下：https://support.microsoft.com/zh-cn/office/gammadist-%E5%87%BD%E6%95%B0-7327c94d-0f05-4511-83df-1dd7ed23e19e?ui=zh-cn&rs=zh-cn&ad=cn
WPS Excel中用法说明如下：https://www.wps.cn/learning/room/d/328021

二、尝试网上找的vb版GammaDist()和GammaInv()两函数的代码转java
1、从matlab里翻译过来的vb版GammaDist()和GammaInv()两函数的代码如下：

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

http://wenku.baidu.com/link?url=uHtUUETc8q1HAlSV7ChTHRg0SLrBKApRKC-BLFnHHoviqv3G4r4DY8ONIWi_5fKN6DuBURHhAkOfJvGNJg7mFyjjERKRoeWASkFTmK3yvXa从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: 
'%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 
 
GammaInv:(要调用上面的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对应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

https://wenku.baidu.com/link?url=uHtUUETc8q1HAlSV7ChTHRg0SLrBKApRKC-BLFnHHoviqv3G4r4DY8ONIWi_5fKN6DuBURHhAkOfJvGNJg7mFyjjERKRoeWASkFTmK3yvXa&_wkts_=1667896077426
https://bbs.co188.com/thread-1169888-1-1.html
2、转成java后代码如下：

import static java.lang.Double.NEGATIVE_INFINITY; //无穷小
import static java.lang.Double.POSITIVE_INFINITY; //无穷大
import static java.lang.Math.*;


public class GammaDist {

    //对应excel中GammaDist(x,a,b,true)函数,
    public GammaDist(double x,double a,double b, boolean flag){
        if(flag){
            double p = GamCdf(x,a,b);
            if(p>1)
                p=1;
            else
                p=p;
            System.out.println (p);
        }
    }
    //%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 double GamCdf(double x, double a, double b) { //对应excel中GammaDist(x,a,b,true)函数,计算S曲线
        double GamCdf = NEGATIVE_INFINITY; //无穷小
        if(a<=0 || b<=0) GamCdf = NEGATIVE_INFINITY; //无穷小
        GamCdf = GAMMAINC(x / b, a);
        double p = GamCdf;
        if(p>1)
            p=1;
        else
            p=p;
        return GamCdf;
    }
             
    //%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 double GAMMAINC(double x, double a) {
        double amax;
        amax = Math.pow(2,20);
        int ascalar = 1;
        double GAMMAINC = 0.0;

        if(a <=amax) {
            if (a != 0 && x != 0 && x < a + 1) {
                double xk = x;
                double ak = a;
                double ap = ak;
                double Sum = 1 / ap;
                double del = Sum;
                for (int i = 0; i < 10000; i++) {
                    ap = ap + 1;
                    del = xk * del / ap;
                    Sum = Sum + del;
                }
                GAMMAINC = Sum * exp(-xk + ak * log(xk) - GAMMALN(ak));
            }
            else if (a != 0 && x != 0 && x >= a + 1) {
                double xk = x;
                double a0 = 1;
                double a1 = x;
                double b0 = 0;
                double b1 = a0;
                double ak = a;
                double fac = 1;
                double n = 1;
                double g = b1;
                double gold = b0;
                for (int i = 0; i < 10000; i++) {
                    gold = g;
                    double ana = n - ak;
                    a0 = (a1 + a0 * ana) * fac;
                    b0 = (b1 + b0 * ana) * fac;
                    double anf = n * fac;
                    a1 = xk * a0 + anf * a1;
                    b1 = xk * b0 + anf * b1;
                    fac = 1 / a1;
                    g = b1 * fac;
                    n = n + 1;
                }
                GAMMAINC = 1 - exp(-xk + ak * log(xk) - GAMMALN(ak)) * g;
            }
        }
        return GAMMAINC;
    }
             
    //'*****************************************************************************************************
    //**************
    //%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 double GAMMALN( double XX) {
        double[] COF = new double[6];
        double stp, half, one;
        double fpf, x, tmp, ser;
        COF[0] = 76.18009173;
        COF[1] = -86.50532033;
        COF[2] = 24.01409822;
        COF[3] = -1.231739516;
        COF[4] = 0.00120858003;
        COF[5] = -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(int j=0;j<6;j++) {
            x = x + one;
            ser = ser + COF[j] / x;
        }
        double GAMMALN = tmp + log(stp * ser);
        return GAMMALN;
    }
             
    //GammaInv:(要调用上面的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值Φp＝Cs/2*GammaInv(1-p/100,4/Cs^2,1)-2/Cs
    public double gaminv(double p, double a, double b) {
        double gaminv = NEGATIVE_INFINITY;//无穷小
        if(p <0 || p >1 || a <= 0 || b <= 0) {
            gaminv = NEGATIVE_INFINITY;//无穷小
            return gaminv;
        }
        else {
            if(p==0){
                gaminv = 0;
            }
            else if(p==1){
                gaminv = 1;
            }
            else {
                //% Newton' s Method
                //% Permit no more than count_limit interations.
                int count_limit = 100;
                int cunt = 0;
                double pk = p;
                //% Supply a starting guess for the iteration.
                //% Use a method of moments fit to the lognormal distribution.
                double mn = a * b;
                double v = mn * b;
                double temp = log(v + Math.pow(mn,2));
                double mu = 2 * log(mn) - 0.5 * temp;
                double sigma = -2 * log(mn) + temp;
                double xk = exp(MyNormInv(pk, mu, sigma));
                //'' '' '' '' '' '' '' '' '' '' '' '
                double 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.
                double eps = Math.pow(2,-10);
                while((abs(h) > Math.pow(eps,0.5) * abs(xk)) && (abs (h) > Math.pow(eps,0.5)) && (cunt <count_limit)){
                    cunt = tunt + 1;
                    h = (GamCdf(xk, a, b) - pk) / gampdf(xk, a, b);
                    //'' '' '' '' '' '' '
                    double 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 / 10 ths of the way to zero:
                    if(xnew <0) {
                        xnew = xk / 10;
                        h = xk - xnew;
                    }
                    xk = xnew;
                }
                gaminv = xk;
            }
        }
    }
             
             
    //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 double MyNormSInv( double p) {
        //计算频率格纸
        double a1 = -39.6968302866538, a2 = 220.946098424521, a3 = -275.928510446969;
        double a4 = 138.357751867269, a5 = -30.6647980661472, a6 = 2.50662827745924;
        double b1 = -54.4760987982241, b2 = 161.585836858041, b3 = -155.698979859887;
        double b4 = 66.8013118877197, b5 = -13.2806815528857, c1 = -7.78489400243029E-03;
        double c2 = -0.322396458041136, c3 = -2.40075827716184, c4 = -2.54973253934373;
        double c5 = 4.37466414146497, c6 = 2.93816398269878, d1 = 7.78469570904146E-03;
        double d2 = 0.32246712907004, d3 = 2.445134137143, d4 = 3.75440866190742;
        double p_low = 0.02425, p_high = 1 - p_low;
        double q, r, MyNormSInv = NEGATIVE_INFINITY;//无穷小;
        if (p <0 || p >1) {
            System.out.println ("NormSInv: Argument out of range.");
        }
        else if( p <p_low){
            q = sqrt(-2 * log(p));
            MyNormSInv = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
        }
        else if(p <=p_high) {
            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 = sqrt(-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);
        }
        return MyNormSInv;
    }
             
    //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 double gampdf(double x, double a, double b){
        double y = 0.0;
        double gampdf = NEGATIVE_INFINITY;//无穷小
        if(x == 0 && a < 1){
            gampdf = POSITIVE_INFINITY;//无穷大
            return gampdf;
        }
        if((x==0)&&(a==1)){
            gampdf = 1 / b;
            return gampdf;
        }
        if(a <= 0 || b <= 0) {
            y = NEGATIVE_INFINITY;//无穷小
            gampdf = y;
        }
        else if (x > 0) {
            y = (a - 1) * log(x) - x / b - GAMMALN(a) - a * log(b);
            y = exp(y);
        }
        gampdf = y;
        return gampdf;
    }

}

3、结论，转换未成功，函数MyNormInv和变量tunt未找到，有知道的请留言告知。

三、最后使用org.apache.commons.math3.distribution.GammaDistribution类中的方法成功。参考地址：https://vimsky.com/examples/detail/java-class-org.apache.commons.math3.distribution.GammaDistribution.html