问题描述
|
在某些迭代优化过程中，以下用于计算双变量正态CDF的VBA代码有时在上层函数的
while
z = hx * hy * c
Function tetrachoric(x As Double,y As Double,rho As Double) As Double
Const FACCURACY As Double = 0.0000000000001
Const MinStopK As Integer = 20
Dim k As Integer
Dim c As Double
Dim z As Double
Dim s As Double
Dim hx As Double
Dim hx1 As Double
Dim hx2 As Double
Dim hy As Double
Dim hy1 As Double
Dim hy2 As Double
Dim CheckPass As Integer
hx = 1
hy = 1
hx1 = 0
hy1 = 0
k = 0
c = rho
z = c
s = z
CheckPass = 0
While CheckPass < MinStopK
k = k + 1
hx2 = hx1
hy2 = hy1
hx1 = hx
hy1 = hy
hx = x * hx1 - (k - 1) * hx2
hy = y * hy1 - (k - 1) * hy2
c = c * rho / (k + 1)
z = hx * hy * c
s = s + z
If Abs(z / s) < FACCURACY Then
CheckPass = CheckPass + 1
Else
CheckPass = 0
End If
Wend
tetrachoric = s
End Function
Public Function bivnor(x As Double,rho As Double) As Double
\'
\' bivnor function
\' Calculates bivariat normal CDF F(x,y,rho) for a pair of standard normal
\' random variables with correlation RHO
\'
If rho = 0 Then
bivnor = Application.WorksheetFunction.normSdist(x) * _
Application.WorksheetFunction.normSdist(y)
Else
bivnor = Application.WorksheetFunction.normSdist(x) * _
Application.WorksheetFunction.normSdist(y) + _
Application.WorksheetFunction.normdist(x,1,False) * _
Application.WorksheetFunction.normdist(y,False) * _
tetrachoric(x,rho)
End If
End Function
解决方法
您正在达到计算机体系结构的极限。由于性能原因和/或溢出时的错误行为，许多复杂的算法无法以其数学表示形式一对一地实现。关于这些问题，有一个非常出色的博客-John D. Cook。
请在此处查看以获得更好的实现。
您还可以尝试绑定外部库，从而为您提供任意精度的数字处理，当然可以使用非常昂贵的（就CPU时间而言）软件算法来实现。在这里可以找到更多。
,使用
On Error Resume Next
On Error Goto
While CheckPass < MinStopK
k = k + 1
hx2 = hx1
hy2 = hy1
hx1 = hx
hy1 = hy
hx = x * hx1 - (k - 1) * hx2
hy = y * hy1 - (k - 1) * hy2
c = c * rho / (k + 1)
On Error Resume Next
z = hx * hy * c
If Err.Number = 0 Then
s = s + z
If Abs(z / s) < FACCURACY Then
CheckPass = CheckPass + 1
Else
CheckPass = 0
End If
Else
Err.Clear
End If
Wend