Function seektheta(p As Single, sig As Single) As Single
Dim x As Single
Dim z1 As Single
Dim y1 As Single
Dim y2 As Single
Dim y3 As Single
Dim g1 As Single
Dim g2 As Single
y1 = -3 * sig
y2 = 3 * sig
g1 = thetag(y1, sig)
t = 0
Do Until t > 100
g2 = thetag(y2, sig)
y3 = y2 + (p - g2) * (y2 - y1) / (g2 - g1)
y1 = y2
y2 = y3
g1 = g2
If (g2 - p) ^ 2 < 10 ^ (-5) Then t = 1000
t = t + 1
Loop
seektheta = Exp(y2)
End Function
Function thetag(y As Single, sig As Single) As Single
Dim n As Single
Dim h As Single
Dim z As Single
Dim x As Single
h = 0.001
z = Int(y / h)
z1 = 0
For n = -5000 To z
x = h * n
z1 = z1 + h * thetaf(x, sig)
Next
thetag = z1
End Function
Function thetaf(x As Single, sig As Single) As Single
Dim pi As Single
Dim x1 As Single
Dim x2 As Single
pi = 3.1415
x1 = Sqr(2 * pi) * sig
x2 = Exp(-x ^ 2 / (2 * sig ^ 2))
thetaf = x2 / x1
End Function
Private Sub Command1_Click()
Dim n As Single
Dim h As Single
Dim a As Single
Dim beta As Single
Dim k(1 To 100) As Single
Dim cx(1 To 10, 1 To 100) As Single
Dim cp(1 To 10, 1 To 100) As Single
Dim th(1 To 10) As Single
Dim ks As Single
Dim r1 As Single
Dim c1 As Single
Dim n1 As Single
Dim n2 As Single
Dim n3 As Single
Dim e As Single
Dim t As Single
Dim uc As Single
Dim sig As Single
Dim p As Single
sig = 0.1
For n = 1 To 10
p = -0.05 + 0.1 * n
th(n) = seektheta(p, sig)
Debug.Print p, th(n)
Next
beta = 0.95
a = 0.33
ks = ((1 / beta - 1) / a) ^ (1 / (a - 1))
h = 2 * ks / 100
For n = 1 To 100
k(n) = n * h
Next
For n = 1 To 100
For m = 1 To 10
cx(m, n) = th(m) * k(n) ^ a
Next
Next
t = 0
Do Until t > 100
For n = 10 To 90
For m = 1 To 10
uc = 0
For s = 1 To 10
k1 = k(n) + th(m) * k(n) ^ a - cx(m, n)
r1 = th(s) * a * k1 ^ (a - 1)
n1 = k1 / h
n2 = Int(n1)
n3 = n2 + 1
c1 = cx(s, n2) + (n1 - n2) * (cx(s, n3) - cx(s, n2))
uc = uc + (beta * (1 + r1)) / c1
Next
uc = uc / 10
cp(m, n) = 1 / uc
Next
Next
e = 0
For n = 10 To 90
For m = 1 To 10
e = e + (cx(m, n) - cp(m, n)) ^ 2
Next
Next
For n = 10 To 90
For m = 1 To 10
cx(m, n) = cp(m, n)
Next
Next
If e < 10 ^ (-5) Then t = 1000
Debug.Print t, e
t = t + 1
Loop
End Sub
最終更新:2009年04月25日 00:06