爱华网标签： 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
