多元逐步回归分析程序源码

你喜欢那个 2012-11-27

展开全文

多元逐步回归分析程序源码

(2011-04-21 14:30:00)

转载▼

标签：

it

分类：统计与编程

声明：任何人可以修改此代码并个人无偿使用此程序，但不得商用，否则将追究法律责任。此文中包含大量共性程序模块，希望有志于统计学程序代码国产化的同仁与我协作，其它大量代码将来适当时机会公开，最后将以教材方式与大家见面。

try段请置于工作表代码中，其它置于类模块中。
Public Qe As Double, U As Double, SSY As Double, VARY As Double, AVEY As Double
Public fEQUATION As Double, pEQUATION As Double, B0 As Double, SB0 As Double
Public inter As Integer, NAMEY As String
Sub try()

    Dim m As Integer, n As Integer
    m = 142: n = 13
    ReDim XY(m, n) As Double, NAME_(n) As String
    Dim I, J
    Dim TOOL As New LINE20110328ok
    For I = 1 To n
        NAME_(I) = Cells(1, I)
        For J = 1 To m
            XY(J, I) = Val(Cells(1 + J, I))
        Next
    Next
    ReDim ave(n) As Double, var(n) As Double, Uvalue(n) As Double, _
     PVALUE(n) As Double, bi(n) As Double, sbi(n) As Double
    TOOL.ONEY_REGRESS XY, n, m, NAME_, ave, var, Uvalue, PVALUE, bi, sbi
    For I = 1 To TOOL.inter
        Cells(10 + m, I + 1) = ave(I)
        Cells(11 + m, I + 1) = var(I)
        Cells(12 + m, I + 1) = bi(I)
        Cells(13 + m, I + 1) = sbi(I)
        Cells(14 + m, I + 1) = Uvalue(I)
        Cells(15 + m, I + 1) = PVALUE(I)
        Cells(9 + m, I + 1) = NAME_(I)

        Cells(10 + m, 1) = "ave(I)"
        Cells(11 + m, 1) = "var(I)"
        Cells(12 + m, 1) = "bi(I)"
        Cells(13 + m, 1) = "sbi(I)"
        Cells(14 + m, 1) = "Uvalue(I)"
        Cells(15 + m, 1) = "PVALUE(I)"
        Cells(9 + m, 1) = "NAME_(I)"

    Next
    Cells(10 + m, TOOL.inter + 3) = TOOL.AVEY
    Cells(11 + m, TOOL.inter + 3) = TOOL.VARY
    Cells(12 + m, TOOL.inter + 3) = TOOL.B0
    Cells(13 + m, TOOL.inter + 3) = TOOL.SB0
    Cells(14 + m, TOOL.inter + 3) = TOOL.fEQUATION
    Cells(15 + m, TOOL.inter + 3) = TOOL.pEQUATION
    Cells(16 + m, TOOL.inter + 3) = TOOL.Qe
    Cells(17 + m, TOOL.inter + 3) = TOOL.SSY

    Cells(10 + m, TOOL.inter + 2) = "AVEY"
    Cells(11 + m, TOOL.inter + 2) = "VARY"
    Cells(12 + m, TOOL.inter + 2) = "B0"
    Cells(13 + m, TOOL.inter + 2) = "SB0"
    Cells(14 + m, TOOL.inter + 2) = "fEQUATION"
    Cells(15 + m, TOOL.inter + 2) = "pEQUATION"
    Cells(16 + m, TOOL.inter + 2) = "Qe"
    Cells(17 + m, TOOL.inter + 2) = "SSY"

End Sub

Sub ONEY_REGRESS(XY() As Double, ByVal Nxy As Integer, ByVal NSAMPLES As Integer, _
    NAME_() As String, ave() As Double, var() As Double, Uvalue() As Double, _
     PVALUE() As Double, bi() As Double, sbi() As Double)
    ''' 单因变量回归方程
    Dim newInter As Integer '''single_y_x_into single_y_x_out 共用变量

    Dim I As Integer, J As Integer, K As Integer
    newInter = 0 '@@@@@@@@@@@
    ReDim Newrelative(Nxy, Nxy) As Double
    ReDim Newr0(Nxy, Nxy) As Double

    simple_r XY, Nxy, NSAMPLES, Newrelative, ave, var

    For I = 1 To Nxy
        For J = 1 To Nxy
            Newr0(I, J) = Newrelative(I, J)
        Next
    Next

    ReDim BASE(Nxy) As Integer

    Do While single_y_x_into(Newr0, Nxy, 0, BASE, Uvalue, newInter, NSAMPLES)              '''核心段代码
        Do While single_y_x_out(Newr0, Nxy, 0, BASE, Uvalue, newInter, NSAMPLES): Loop '''无变量可出基也无变量可入基时
    Loop                                                '''才终止
    '''以下为结果数据整理程序

    inter = 0

    For I = 1 To Nxy - 1
        If BASE(I) Then inter = inter + 1
    Next
    '''拷回原始数据重新计算，以提高精度----------
    For I = 1 To Nxy
        For J = 1 To Nxy
            Newr0(I, J) = Newrelative(I, J)
        Next
    Next

    For K = 1 To Nxy - 1
        If BASE(K) Then 紧凑变换 Newr0, Nxy, K
    Next

    '''-------------拷回原始数据重新计算，以提高精度
      ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

    ReDim LXX_1(inter, inter) As Double, AVEX(inter) As Double, VARX(inter) As Double, _
        NAMEX(inter) As String, UVALUEX(inter) As Double, PVALUEX(inter) As Double '''生成Lxx-1及相关均值向量
    ReDim order(inter + 1) As Integer '''序数，拷贝数据用
    Dim II As Integer: II = 1:
    For I = 1 To Nxy - 1
        If BASE(I) Then
            order(II) = I

            II = II + 1
        End If
    Next
    order(inter + 1) = Nxy
    For I = 1 To inter '''生成Lxx-1及相关参数向量
        AVEX(I) = ave(order(I)): NAMEX(I) = NAME_(order(I)): VARX(I) = var(order(I))
        UVALUEX(I) = Uvalue(order(I))
        For J = 1 To inter
            LXX_1(I, J) = Newr0(order(I), order(J))
        Next
    Next

    For I = 1 To inter '''Lxx-1还原为协方差分析结果
        For J = 1 To inter
            LXX_1(I, J) = LXX_1(I, J) / VARX(I) / VARX(J)
        Next
    Next
    For I = 1 To Nxy
        Newr0(I, Nxy) = Newr0(I, Nxy) * var(Nxy) / var(I)
        Newr0(Nxy, I) = Newr0(Nxy, I) * var(Nxy) / var(I)
    Next


    SSY = (var(Nxy)) ^ 2 '''相关阵分析时
    Qe = Newr0(Nxy, Nxy) * SSY
    U = SSY - Qe

    For I = 1 To inter
        sbi(I) = Sqr(LXX_1(I, I) * Qe / (NSAMPLES - inter - 1))
        bi(I) = -1 * Newr0(Nxy, order(I))
        PVALUEX(I) = vbaFDIST(UVALUEX(I) / Qe * (NSAMPLES - inter - 2) * SSY, 1, (NSAMPLES - inter - 2))
        UVALUEX(I) = UVALUEX(I) / Qe * (NSAMPLES - inter - 2) * SSY
    Next

    '''计算SB0,B0 AND F OF THE EQUATION
    B0 = ave(Nxy) '''NOT DIM
    For I = 1 To inter
        B0 = B0 - AVEX(I) * bi(I)
    Next


    ReDim TT(1, inter) As Double, xt(1, inter) As Double
    For I = 1 To inter
        xt(1, inter) = AVEX(I)
    Next

    m_mult xt, 1, inter, LXX_1, inter, inter, TT
    ReDim xt(inter, 1) As Double
    Dim T(1, 1) As Double '''TEMP VAR
    For I = 1 To inter
        xt(I, 1) = AVEX(I)
    Next

    m_mult TT, 1, inter, xt, inter, 1, T

    SB0 = T(1, 1) '''返回矩阵运算结果
    SB0 = SB0 + 1 / NSAMPLES
    SB0 = Sqr(SB0 * Qe / (NSAMPLES - inter - 1))
    fEQUATION = U / inter / Qe * (NSAMPLES - inter - 1)
    pEQUATION = vbaFDIST(fEQUATION, inter, NSAMPLES - inter - 1)

    For I = 1 To inter '''原地返回重要参数 X 部分
         ave(I) = AVEX(I)
         NAME_(I) = NAMEX(I)
         var(I) = VARX(I)
         Uvalue(I) = UVALUEX(I)
         PVALUE(I) = PVALUEX(I)
    Next
    AVEY = ave(Nxy)
    VARY = var(Nxy)
End Sub

Function single_y_x_into(R() As Double, ByVal Nxy As Integer, _
    ByVal test As Integer, BASE() As Integer, Uvalue() As Double, _
    newInter As Integer, ByVal NSAMPLES) As Integer
   'MsgBox "single_y_x_into" '''判断x是否引入
    Dim maxI As Integer, into As Integer
    Dim Max As Double, U As Double, Qe As Double, F As Double, P As Double
    Dim I As Integer, df As Integer
    '''统计进入方程自变量与未进入方程的自变量的F值比较结果
    maxI = 0: Max = 0
    For I = 1 To Nxy - 1 '''填F值
        U = R(I, Nxy) ^ 2 / R(I, I)
        Uvalue(I) = U
    Next
    into = 0
    For I = 1 To Nxy - 1 '''统计出基最大值
        If BASE(I) = 1 Then
            into = into + 1
        Else
            If Uvalue(I) > Max Then
                Max = Uvalue(I)
                maxI = I
            End If
        End If
    Next

    If into Then
        Qe = R(Nxy, Nxy)
        U = 1 - Qe
        Uvalue(Nxy) = U / into / Qe * (n - into - 1)
    End If
    If into = Nxy - 1 Then '''满员退出
        single_y_x_into = 0
        Exit Function
    End If

    F = (NSAMPLES - into - 2) * Uvalue(maxI) / (R(Nxy, Nxy) - Uvalue(maxI))

    If vbaFDIST(F, 1, NSAMPLES - into - 2) < 0.05 Then
        BASE(maxI) = 1 '''基变量入基
        newInter = maxI

        If test Then '若不要求执行则退出
            single_y_x_into = 1
            Exit Function
        End If

        紧凑变换 R, Nxy, maxI
        single_y_x_into = 1
    Else
        single_y_x_into = 0
    End If
End Function
Function single_y_x_out(R() As Double, ByVal Nxy As Integer, _
    ByVal test As Integer, BASE() As Integer, Uvalue() As Double, _
    newInter As Integer, ByVal NSAMPLES As Integer)
     'MsgBox "single_y_x_out""判断x是否退出"
    Dim minI As Integer
    Dim Min As Double
    Dim into As Integer
    Dim U As Double, Qe As Double, F As Double, P As Double
    Dim I As Integer, df As Integer

    '''统计进入方程自变量与未进入方程的自变量的F值比较结果
    minI = 0: Min = 100000
    For I = 1 To Nxy - 1 '''填U值
        U = R(I, Nxy) ^ 2 / R(I, I)
        Uvalue(I) = U
    Next
    into = 0
    For I = 1 To Nxy - 1 '''统计入基最小值
        If BASE(I) = 1 Then
            into = into + 1
            If Uvalue(I) < Min Then
                Min = Uvalue(I)
                minI = I
            End If
        End If
    Next

    If newInter = minI And minI <> 0 Then
        single_y_x_out = 0 '''若最小变量为新引入变量，则退出
        Exit Function
    End If

    If into Then
        Qe = R(Nxy, Nxy)
        U = 1 - Qe
        Uvalue(Nxy) = U / into / Qe * (n - into - 1) '''方程显著性
    End If
    If into = 0 Or into = 1 Then
        single_y_x_out = 0 '无变量进入或仅有一个变量进入则退出
        Exit Function
    End If

     F = (NSAMPLES - into - 1) * Uvalue(minI) / R(Nxy, Nxy)

    If vbaFDIST(F, 1, NSAMPLES - into - 1) > 0.1 Then
        BASE(minI) = 0 '''变量出基
        newInter = 0
        If test Then '''若不要求执行则退出
            single_y_x_out = 1
            Exit Function
        End If
        If minI Then 紧凑变换 R, Nxy, minI
        single_y_x_out = 1
        Exit Function
    Else
        single_y_x_out = 0
    End If
End Function

'''**************************************************************************
'''紧凑变换计算程序                          2011.03.28
'''xy() 原始数据, nxy 变量个数,   Nsamples观察值组数
'''**************************************************************************
Sub 紧凑变换(R() As Double, ByVal Num As Integer, ByVal main_n As Integer)
    Dim I As Integer, J As Integer, Main As Double
    ReDim mainrow(Num) As Double, maincolumn(Num) As Double '''
    For I = 1 To Num '''read main row and column value     2011.03.28
        mainrow(I) = R(main_n, I)
        maincolumn(I) = R(I, main_n)
    Next
    Main = R(main_n, main_n)

    For I = 1 To Num '''all
        For J = 1 To Num
            If I = main_n And J = main_n Then
                R(I, J) = 1 / Main
            ElseIf I = main_n And J <> main_n Then
                R(I, J) = R(I, J) / Main
            ElseIf I <> main_n And J = main_n Then
                R(I, J) = R(I, J) / Main * -1
            Else
                R(I, J) = R(I, J) - mainrow(J) * maincolumn(I) / Main
            End If
        Next
    Next

End Sub

'''**************************************************************************
'''简单相关计算程序 2011.03.29
'''xy() 原始数据, nxy 变量个数, Nsamples观察值组数
'''ave() 均值 , var() 方差
'''**************************************************************************

Sub simple_r(XY() As Double, ByVal Nxy As Integer, ByVal NSAMPLES As Integer, _
    r0() As Double, ave() As Double, var() As Double)
    Dim I, J, K As Integer
    Dim sx As Double, sx2 As Double, sy As Double, sy2 As Double, sxy As Double
    For I = 1 To Nxy '''求简单相关系数

        For J = I To Nxy
            If I = J Then
                r0(I, J) = 1
            Else

                sx = 0
                sx2 = 0
                sy = 0
                sy2 = 0
                sxy = 0

                For K = 1 To NSAMPLES
                    sx = sx + XY(K, I)
                    sx2 = sx2 + XY(K, I) * XY(K, I)
                    sy = sy + XY(K, J)
                    sy2 = sy2 + XY(K, J) * XY(K, J)
                    sxy = sxy + XY(K, I) * XY(K, J)
                Next

                r0(I, J) = (sxy - sx * sy / NSAMPLES) / Sqr((sx2 - sx * sx / NSAMPLES) * (sy2 - sy * sy / NSAMPLES))
                r0(J, I) = r0(I, J)
            End If

        Next
    Next
    For I = 1 To Nxy
        sx = 0
        sx2 = 0
        For K = 1 To NSAMPLES
            sx = sx + XY(K, I)
            sx2 = sx2 + XY(K, I) * XY(K, I)
        Next
        ave(I) = sx / NSAMPLES
        var(I) = Sqr((sx2 - sx ^ 2 / NSAMPLES) / (NSAMPLES - 1))
    Next

End Sub
'''**************************************************************************
'''协方差计算程序                          2011.03.28
'''xy() 原始数据, nxy 变量个数,   Nsamples观察值组数
'''**************************************************************************

Sub SP_XXY(XY() As Double, ByVal Nxy As Integer, ByVal NSAMPLES As Integer, sp() As Double)
    Dim I, J, K As Integer
    Dim sx As Double, sx2 As Double, sy As Double, sy2 As Double, sxy As Double
    For I = 1 To Nxy '''求简单相关系数
        For J = I To Nxy
            If I = J Then
                sx = 0
                sx2 = 0
                For K = 1 To NSAMPLES
                    sx = sx + XY(K, I)
                    sx2 = sx2 + XY(K, I) * XY(K, I)
                Next
                sp(I, J) = (sx2 - sx * sx / NSAMPLES)
            Else

                sx = 0
                sx2 = 0
                sy = 0
                sy2 = 0
                sxy = 0

                For K = 1 To NSAMPLES
                    sx = sx + XY(K, I)
                    sx2 = sx2 + XY(K, I) * XY(K, I)
                    sy = sy + XY(K, J)
                    sy2 = sy2 + XY(K, J) * XY(K, J)
                    sxy = sxy + XY(K, I) * XY(K, J)
                Next

                sp(I, J) = (sxy - sx * sy / NSAMPLES)
                sp(J, I) = r0(I, J)
            End If

        Next
    Next
End Sub

'''***********************索引数组生成程序**************************
'''2011.03.07 by Chen Tingmu RIGHT
'''Dim Nindex As LONG 全局变量，最多10个变量，20个水平
'''Dim indexNAME(10) As String 全局变量，最多10个变量，20个水平
'''Dim indexMAX(10) As Long全局变量
'''Dim indexORDER(10, 20) As String全局变量
'''Dim Nsample As Long 样本总数
'''Dim INDEXtable(10,500) As integer 索引化因子表，最多500组值

'''Dim aNAME As Range 由程序中指定作为输入因子及重复名称用
'''Dim aVALUE As Range 由程序中指定作为输入观测值数据用
'''*****************************************************************
Dim Nindex As Long
Dim indexNAME(20) As String
Dim indexMAX(20) As Long
Dim indexORDER(20, 20) As String
Dim Nsample As Long
Dim INDEXtable(10, 500) As Integer

Sub index(ByVal aNAME As Range, Nindex As Long, Nsample As Long, indexMAX() As Long, indexORDER() As String)
'''aName是字符型数据,aIndex是整数型数据
    Dim I, J As Integer
    Dim isadd As Boolean
    isadd = False: Nindex = 1
    Nsample = aNAME.Rows.Count '''samples numbers
    Nindex = aNAME.Columns.Count
    For I = 1 To Nindex
        indexMAX(I) = 1         '''initialize
        indexORDER(I, 1) = aNAME(1, I)
    Next

    For I = 1 To Nindex
        For J = 1 To Nsample '''can be modyfied
            isadd = True
            For K = 1 To indexMAX(I)
                If aNAME(J, I) = indexORDER(I, K) Then
                    isadd = False: Exit For
                End If
            Next

            If isadd Then
                indexMAX(I) = indexMAX(I) + 1
                indexORDER(I, indexMAX(I)) = aNAME(J, I)
            End If
            'MsgBox indexORDER(i, indexMAX(i))
        Next
    Next
End Sub

'''***********************处理水平名称索引化算法程序**************************
'''由索引数组生成
'''2011.03.07 by Chen Tingmu                                       right
'''Dim INDEXtable(10,500) As integer 索引化因子表，最多500组值
'''***************************************************************************
Sub ToIndex(aNAME As Range, INDEXtable() As Integer, Nindex As Long, Nsample As Long, indexMAX() As Long, indexORDER() As String)
    Dim I, J, K As Long
    For I = 1 To Nindex
        For J = 1 To Nsample
            For K = 1 To indexMAX(I)
                If aNAME(J, I) = indexORDER(I, K) Then
                    INDEXtable(I, J) = K
                    Exit For
                End If
            Next
        Next
    Next
End Sub

'''***********************域排序程序:冒泡法**************************
'''多重比较专用,其它用途略改可
'''Data As range,待排序数组，以列分组，不含表头
'''ByVal Nrows As Integer，Data()行数
'''ByVal Ncolumns As Integer，Data()列数
'''ByVal sortField As Integer，Data()第sortField用于排序，从大到小排序
'''2011.03.08 by Chen Tingmu right 14：23

'''***********************************************************
Sub BlockSort(ByVal Data As Range, ByVal sortField As Integer)
    Dim I, J, intT, m As Integer
    Dim maxValue, minValue, doubleT As Double
    Dim maxI, minI As Integer
    Dim index() As Integer, Ndata As Integer
    Dim indexValue() As Double
    Ndata = Data.Rows.Count
    ReDim indexValue(Ndata), index(Ndata)
    For I = 1 To Ndata '''initialize index(i) = i: indexValue(i) = Data(i, sortField)
        index(I) = I: indexValue(I) = Data(I, sortField)
    Next

    I = 1
    Do While I <= Ndata '''?????????????????????????
        maxValue = indexValue(I)
        maxI = I
        J = I
        Do While J <= Ndata '''find the max and min value and order

            If maxValue <= indexValue(J) Then
                maxValue = indexValue(J): maxI = (J)
                'MsgBox maxValue
            End If
            J = J + 1
        Loop

        doubleT = indexValue(I) '''change the larger numble??????????????????
        indexValue(I) = indexValue(maxI)
        indexValue(maxI) = doubleT

        intT = index(I)
        index(I) = index(maxI)
        index(maxI) = intT

        I = I + 1

    Loop


'''index and indexvalue macthing correct

    For I = 1 To Ndata
        Data(I, 4) = index(I)
        Data(I, 5) = indexValue(I)
    Next

End Sub

Private Function GAMA(ByVal alfa As Double) As Double '''gama分布积分函数,correct
    Dim I, nalfa As Long
    Dim done, half As Boolean: done = False: half = False
    nalfa = Int(alfa)
    If VBA.Abs((alfa - nalfa)) < 0.1 Then half = False Else half = True
    Dim F As Double: F = 1
    If half Then
        For I = 1.5 To alfa - 1 Step 1
            F = F * I
        Next
        F = F * VBA.Sqr(3.14159265358979)
    Else
        For I = 1 To alfa - 1 Step 1
            F = F * I
        Next
    End If
    GAMA = F
End Function

'''======================================================
'''函数名：vbaTDIST.c    20110310 BY CHENTINGMU RIGHT
'''功能描述：用t分布检验两个分布的均值是否有显著性差异(方差相同)
'''输入参数：a（简单随机样本一的样本值）
'''      na（样本一的个数）
'''          b（简单随机样本二的样本值）
'''      nb（样本二的个数）
'''      alpha（显著性标准）
'''返回值：1（显著），0（不显著）********实际应用中取补概率********
'''=========================================================

Function vbaTDIST(ByVal T As Double, ByVal F As Integer) As Double

    Dim x, P As Double
    Dim V As Integer
    Dim eps As Double: eps = 0.000000000001

    V = F

    x = V / (V + T * T)

    P = beta2(V / 2, 0.5, x, eps)

    vbaTDIST = P

End Function

'''======================================================
'''函数名：vbaFDIST.c
'''功能描述：用f分布检验两个分布的方差是否有显著性差异
'''输入参数：a（简单随机样本一的样本值）
'''      na（样本一的个数）
'''          b（简单随机样本二的样本值）
'''      nb（样本二的个数）
'''      alpha（显著性标准）
'''返回值：1（显著），0（不显著）********实际应用中取补概率********
'''=======================================================

Function vbaFDIST(ByVal F As Double, ByVal V1 As Integer, ByVal V2 As Integer) As Double

    Dim eps As Double: eps = 0.000000000001
    Dim Q1 As Double
    Q1 = 2# * beta2(0.5 * V2, 0.5 * V1, V2 / (V2 + V1 * F), eps) ''' 调用不完全贝塔函数计算
    If (Q1 > 1#) Then Q1 = 2# - Q1
    vbaFDIST = (Q1) / 2
End Function


'''=============================================================
''' 函数名：vbaCHIDIST
''' 功能描述：求解卡方分布函数的值 2010.06.22 RIGHT
''' 输入参数：x 自变量x的值。
'''           n 整数
''' 返回值：卡方分布函数的值********实际应用中取补概率********
'''==============================================================

Function vbaCHIDIST(ByVal x As Double, ByVal n As Integer) As Double
Dim T As Double

Const eps = 0.000000001 ''' 使用不完全伽马函数需要的数据
Const DMIN = 1E-30

If (x < 0) Then                         ''' 若x<0，则不能计算
    MsgBox ("negative x!")
    chii = (0#)
    Exit Function
End If
T = gamm2(n / 2#, x / 2#, eps, DMIN) ''' 调用不完全伽马函数求值
vbaCHIDIST = 1 - (T)
End Function
'''=============================================================
''' 函数名：gamm2
''' 功能描述：求解不完全伽马函数的值 2010.06.22 RIGHT
''' 输入参数：a 自变量a的值。要求a>0。
'''           x 自变量x的值，要求x>=0。
'''         e1 精度要求，当两次递推的值变化率小于e1时，认为已收敛
'''         e0 极小的数值，接近浮点数能表示的最小数据。为避免除零，将除数设置的值
''' 返回值：不完全伽马函数的值
'''==============================================================

Private Function gamm2(ByVal a, x, e1, e0 As Double) As Double
    Const NMAX = 20000

     Dim n As Integer
     Dim T, del, gln As Double
     Dim an, bn, c, d As Double               ''' 计算连分式级数需要的变量
     If ((x < 0#) Or (a <= 0)) Then
       MsgBox ("x<0.0 or a<=0.0!")
       gamm2 = (0#)
       Exit Function
     End If
     If (x < e0) Then
       gamm2 = (0#)
       Exit Function
     End If
     gln = gammln(a)
     If (x < (a + 1#)) Then            ''' 调用求和级数

       del = 1# / a                '''gamm(a)/gamm(a+1)=1/a
       T = 1# / a
       For n = 1 To NMAX - 1

         del = del * x / (a + n)
         T = T + del
         If (Abs(del) < Abs(T) * e1) Then     ''' 级数部分已经收敛

           T = T * Exp(-x + a * Log(x) - gln)
           gamm2 = (T)
           Exit Function
         End If
       Next
       MsgBox (" iteration too many times!")      ''' 经过NMAX次迭代没有收敛
       gamm2 = (0#)
       Exit Function

     Else                          ''' 使用连分式级数

       bn = x + 1# - a              ''' 已经计算了第一节连分式
       c = 1# / e0
       d = 1# / bn
       T = d
       For n = 1 To NMAX - 1

         an = n * (a - n)            ''' 此节的系数a
         bn = bn + 2#                ''' 此节的系数b
         d = an * d + bn
         c = bn + an / c
         If (Abs(d) < e0) Then d = e0        ''' 若小于e0，则认为是0

         If (Abs(c) < e0) Then c = e0

         d = 1# / d
         del = d * c
         T = T * del
         If (Abs(del - 1#) < e1) Then     ''' 级数部分已经收敛

           T = Exp(-x + a * Log(x) - gln) * T
           T = 1# - T
           gamm2 = (T)
           Exit Function
         End If
       Next
       MsgBox (" iteration too many times!")      ''' 经过NMAX次迭代没有收敛
       gamm2 = (0#)
       Exit Function
     End If
End Function

'''=============================================================
''' 函数名：beta2
''' 功能描述：求解不完全贝塔积分的值2010.06.22 RIGHT
''' 输入参数：a 自变量a的值。要求a>0。
'''           b 自变量b的值。要求b>0。
'''           x 自变量x的值，要求0<=x<=1。
'''         e1 精度要求，当两次递推的值变化率小于e1时，认为已收敛
''' 返回值：不完全贝塔函数的值
'''==============================================================

Private Function beta2(ByVal a As Double, ByVal b As Double, ByVal x As Double, _
    ByVal e1 As Double) As Double
Dim T As Double
          ''' 计算连分式级数需要的变量和函数
If ((x < 0#) Or (x > 1#) Or (a <= 0#) Or (b <= 0#)) Then

    MsgBox ("Bad input parameter!")
    beta2 = (0#)
    Exit Function

ElseIf (x = 0#) Then                               ''' x为0的情况

    T = 0#
    beta2 = (T)
    Exit Function

ElseIf (x = 1#) Then                               ''' x为1的情况

    T = 1#
    beta2 = (T)
    Exit Function

ElseIf (x > (a + 1#) / (a + b + 2#)) Then
'MsgBox "(x > (a + 1#) / (a + b + 2#))"
    T = Exp(gammln(a + b) - gammln(a) - gammln(b) + a * Log(x) + b * Log(1# - x)) ''' 系数
    T = 1# - T * subcf(b, a, 1# - x, e1) / b        ''' 使用连分式级数
    beta2 = (T)
    Exit Function

Else
'MsgBox " Else gammln(a + b)" & gammln(a + b)
    T = Exp(gammln(a + b) - gammln(a) - gammln(b) + a * Log(x) + b * Log(1# - x)) ''' 系数
    T = T * subcf(a, b, x, e1) / a                  ''' 使用连分式级数
    beta2 = (T)
    Exit Function
End If
'MsgBox "beta2 "

End Function

'''''''Function subcf 2010.06.22 RIGHT
Static Function subcf(ByVal a As Double, ByVal b As Double, ByVal x As Double, _
    ByVal e1 As Double) As Double
Const NMAX = 20000 '''2011.03.10 ADD CTM
Const FPMIN = 1E-30
Dim n As Integer
Dim T, del, an, c, d As Double
c = 1#
d = 1# - (a + b) * x / (a + 1#)
If (Abs(d) < FPMIN) Then d = FPMIN
d = 1# / d
T = d
For n = 1 To NMAX - 1

    an = n * (b - n) * x / ((a + 2# * n - 1#) * (a + 2# * n)) ''' 第2n节的系数a,此节的系数b为1
    d = an * d + 1#                             ''' 计算d
    c = 1# + an / c                             ''' 计算c
    If (Abs(d) < FPMIN) Then d = FPMIN                    ''' 检查cd的范围

    If (Abs(c) < FPMIN) Then c = FPMIN
    d = 1# / d
    del = d * c
    T = T * del
    an = -(a + n) * (a + b + n) * x / ((a + 2# * n) * (a + 1# + 2# * n)) ''' 第2n+1节
    d = 1# + an * d
    c = 1# + an / c
    If (Abs(d) < FPMIN) Then d = FPMIN
    If (Abs(c) < FPMIN) Then c = FPMIN
    d = 1# / d
    del = d * c
    T = T * del
    If (Abs(del - 1#) < e1) Then                  ''' 级数部分已经收敛
      subcf = (T)
      Exit Function
    End If
Next
' MsgBox "subcf "

MsgBox ("iteration not converged.")            ''' 没有收敛
subcf = (T)
End Function

'''=============================================================
''' 函数名：gammln
''' 功能描述：求解伽马函数的值的自然对数 2010.06.22 RIGHT
''' 输入参数：x 求值的自变量
''' 返回值：伽马函数的值的自然对数
'''==============================================================

Private Function gammln(ByVal x As Double) As Double
'MsgBox ("gammln")
Dim I As Integer
Dim T, S As Double
Static c(6) As Double
    c(1) = 76.1800917294715

    c(2) = -86.5053203294168
    c(3) = 24.0140982408309
    c(4) = -1.23173957245015
    c(5) = 1.20865097386618E-03
    c(6) = -5.395239384953E-06             ''' 系数
If (x < 0) Then

    MsgBox ("incorrect input parameter!")
    gammln = (0#)
    Exit Function
End If
S = 1.00000000019001
For I = 1 To 6                            ''' 级数和
    S = S + c(I) / (x + I - 1)

Next
T = x + 4.5                               ''' 已取对数
T = T - (x - 0.5) * VBA.Log(T)
T = Log(2.506628274631 * S) - T           ''' 最后结果
gammln = (T)

End Function

'''***********************************************************************************
'''矩阵乘法，同时处理向量与向量，向量与矩阵间，矩阵与矩阵间乘法运算
'''a() As Double,前矩阵
'''ByVal MA As Integer, ByVal NA As Integer 前矩阵规格
''' b() As Double,后矩阵
'''ByVal MB As Integer,    ByVal NB As Integer 后矩阵规格
'''RESULT() As Double 结果矩阵
'''2011.03.17 by chentingmu                              right
'''向量当作特别二维数组，方便几类方法的融合，不可用一维数组表示（会越界）
'''***********************************************************************************
Function m_mult(a() As Double, ByVal MA As Integer, ByVal NA As Integer, b() As Double, ByVal MB As Integer, _
    ByVal NB As Integer, RESULT() As Double) As Boolean   '''right
    Dim I As Long, J As Long, K As Long '''2011.03.17 CTM,重新修正，能处理行向量或列向量.

    If NA <> MB Then '''不合要求
        m_mult = False
        Exit Function
    ElseIf MA = 1 And NA > 1 And MB > 1 And NB = 1 Then '''a为行向量，b为列向量
         Dim T As Double: T = 0
         For I = 1 To NA '''RIGHT
                T = T + a(1, I) * b(I, 1)
         Next
         RESULT(1, 1) = T
         m_mult = True
    ElseIf MA > 1 And NA = 1 And MB = 1 And NB > 1 Then '''a为列向量，b为行向量
         For I = 1 To MA
            For J = 1 To NB
                RESULT(I, J) = a(I, 1) * b(1, J)
            Next
         Next
         m_mult = True
    ElseIf MA = 1 And NA > 1 And MB > 1 And NB > 1 Then '''a为行向量，b为普通矩阵

         For I = 1 To NB '''RIGHT
            RESULT(1, I) = 0
            For J = 1 To MB
                RESULT(1, I) = RESULT(1, I) + a(1, J) * b(J, I)
            Next
         Next
         m_mult = True

    ElseIf MA > 1 And NA > 1 And MB > 1 And NB = 1 Then '''a为普通矩阵，b为列向量

         For I = 1 To MA
            RESULT(I, 1) = 0
            For J = 1 To NA
                RESULT(I, 1) = RESULT(I, 1) + a(I, J) * b(J, 1)
            Next
         Next
         m_mult = True
    ElseIf MA > 1 And NA > 1 And MB > 1 And NB > 1 Then '''a，b均为普通矩阵
        For I = 1 To MA '''RIGHT
            For J = 1 To NB
                RESULT(I, J) = 0
                For K = 1 To NA
                    RESULT(I, J) = RESULT(I, J) + a(I, K) * b(K, J)
                Next
            Next
        Next
        m_mult = True
    End If

End Function

Sub TitleStyle(a As Range)
    a.Select
    With Selection
        .HorizontalAlignment = xlLeft
        .VerticalAlignment = xlCenter
        .WrapText = False
        .Orientation = 0
        .AddIndent = False
        .IndentLevel = 0
        .ShrinkToFit = False
        .ReadingOrder = xlContext
        .MergeCells = True
        .Font.Bold = True
      End With
End Sub

Sub RowHeadStyle()
    Columns("a").Select
    With Selection
        .HorizontalAlignment = xlLeft
        .VerticalAlignment = xlCenter
        .WrapText = False
        .Orientation = 0
        .AddIndent = False
        .IndentLevel = 0
        .ShrinkToFit = False
        .ReadingOrder = xlContext
        .MergeCells = False

    End With
    Range("A1:BA65536").Select
    With Selection.Font
        .NAME = "宋体"
        .FontStyle = "常规"
        .Size = 9
        .Strikethrough = False
        .Superscript = False
        .Subscript = False
        .OutlineFont = False
        .Shadow = False
        .Underline = xlUnderlineStyleNone
        .ColorIndex = xlAutomatic
    End With
End Sub