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

(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)
多元逐步回归分析程序源码 spss多元逐步线性回归

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

  

爱华网本文地址 » http://www.413yy.cn/a/25101013/157061.html

更多阅读

股票公式全解析:16 股票软件源码引入

股票公式全解析:[16]股票软件源码引入——简介我们上一篇文章主要说明了大智慧的源码引入的基本方法,大家在使用这个源码的时候一定要注意具体的设置,按照我的要求一步一步实现,源码的编写我会专门有一个介绍,现在我继续说明其他股票软件

spss教程:线性回归分析

spss教程:线性回归分析——简介回归分析是一种应用很广的数量分析方法,用于分析事物间的统计关系,侧重数量关系变化。回归分析在数据分析中占有比较重要的位置。一元线性回归模型:指只有一个解释变量的线性回归模型,用来揭示被解释变量

易语言进度条源码 精易论坛

易语言进度条源码——简介今天,我给大家带来如何弄进度条!易语言进度条源码——工具/原料电脑易语言易语言进度条源码——方法/步骤易语言进度条源码 1、打开易语言!拉

声明:《多元逐步回归分析程序源码 spss多元逐步线性回归》为网友雲煙雨巷分享!如侵犯到您的合法权益请联系我们删除