utamaru - プログラム倉庫

Function seekbud(th, c, y) As Single
Dim m As Single
Dim s As Single
Dim b1 As Single
b1 = 0
For s = 1 To 200
m = prefer(s, th, c, y)
b1 = b1 + y(m) - c(m)
Next
seekbud = b1
End Function
Function seekwel(th, c, y) As Single
Dim m As Single
Dim s As Single
Dim w1 As Single
w1 = 0
For s = 1 To 200
m = prefer(s, th, c, y)
w1 = w1 + Log(c(m)) + Log(1 - y(m) / th(s))
Next
seekwel = w1
End Function
Function wel(tl As Single, tr As Single, th, y) As Single
Dim c(0 To 20) As Single
Dim m As Single
Dim s As Single
Dim b1 As Single
For m = 0 To 20
c(m) = (1 - tl) * y(m) + tr
Next
w1 = 0
For s = 1 To 200
m = prefer(s, th, c, y)
w1 = w1 + Log(c(m)) + Log(1 - y(m) / th(s))
Next
wel = w1
End Function

Function trs(tl As Single, th, y) As Single
Dim tr1 As Single
Dim tr2 As Single
Dim tr3 As Single
Dim b1 As Single
Dim b2 As Single
Dim t As Single
tr1 = 0.01
tr2 = 0.5
b1 = bud(tl, tr1, th, y)
b2 = bud(tl, tr2, th, y)
t = 0
Do Until t > 100
tr3 = (tr1 + tr2) / 2
b3 = bud(tl, tr3, th, y)
If b3 < 0 Then tr2 = tr3
If b3 > 0 Then tr1 = tr3
If b3 ^ 2 < 10 ^ (-3) Then t = 1000
t = t + 1
Loop
trs = tr1
End Function
Function bud(tl As Single, tr As Single, th, y) As Single
Dim c(0 To 20) As Single
Dim m As Single
Dim s As Single
Dim b1 As Single
For m = 0 To 20
c(m) = (1 - tl) * y(m) + tr
Next
b1 = 0
For s = 1 To 200
m = prefer(s, th, c, y)
b1 = b1 + y(m) - c(m)
Next
bud = b1
End Function

Function prefer(s As Single, th, c, y) As Single
Dim m As Single
Dim c1 As Single
Dim l1 As Single
Dim pp As Single
Dim u1 As Single
Dim maxu As Single
Dim maxm As Single
maxu = -999
For m = 0 To 20
c1 = c(m)
l1 = y(m) / th(s)
pp = 0
If l1 > 0.99 Then pp = 1
If pp = 1 Then l1 = 0.5
u1 = Log(c1) + Log(1 - l1)
If pp = 1 Then u1 = -999
If u1 > maxu Then maxm = m
If u1 > maxu Then maxu = u1
Next
prefer = maxm
End Function
Private Sub Command1_Click()
Dim th(1 To 200) As Single
Dim y(0 To 20) As Single
Dim curve(0 To 20) As Single
Dim c(0 To 20) As Single
Dim s As Single
Dim m As Single
Dim n As Single
Dim tl As Single
Dim tr As Single
Dim w1 As Single
Dim b1 As Single
Dim cps As Single
Dim cp1 As Single
Dim cp2 As Single
Dim h As Single
Dim mx As Single
Dim maxw As Single
Dim maxtl As Single
Dim opc(0 To 20) As Single
Dim opw(0 To 20) As Single
For s = 1 To 200
th(s) = 0.01 * s
Next
For m = 0 To 20
y(m) = 0.05 * m
Next
maxw = -999
For n = 1 To 50
tl = 0.01 * n
tr = trs(tl, th, y)
w1 = wel(tl, tr, th, y)
If w1 > maxw Then maxtl = tl
If w1 > maxw Then maxw = w1
Next
tl = maxtl
tr = trs(tl, th, y)
For m = 0 To 20
curve(m) = (1 - tl) * y(m) + tr
Next
Debug.Print tl, tr
t2 = 0
Do Until t2 > 10
For m = 0 To 20
t1 = 0
h = 0.01
For mx = 0 To 20
c(mx) = curve(mx)
Next
cps = curve(m)
Do Until t1 > 100
cp1 = cps + h
cp2 = cps - h
c(m) = cps
bs = seekbud(th, c, y)
ws = seekwel(th, c, y)
If bs < 0 Then ws = -999
c(m) = cp1
b1 = seekbud(th, c, y)
w1 = seekwel(th, c, y)
If b1 < 0 Then w1 = -999
c(m) = cp2
b2 = seekbud(th, c, y)
w2 = seekwel(th, c, y)
If b2 < 0 Then w2 = -999
If w1 > ws Then cps = cp1
If w1 > ws Then ws = w1
If w2 > ws Then cps = cp2
If w2 > ws Then ws = w2
If (cps - cpx) ^ 2 < 10 ^ (-5) Then h = h / 2
If h < 10 ^ (-7) Then t1 = 1000
cpx = cps
t1 = t1 + 1
Loop
opc(m) = cps
opw(m) = ws
Next
maxw = -999
For m = 0 To 20
If opw(m) > maxw Then mx = m
If opw(m) > maxw Then maxw = opw(m)
Next
curve(mx) = opc(mx)
t2 = t2 + 1
Debug.Print t2, ws
Loop
For s = 1 To 199
m1 = prefer(s, th, curve, y)
m2 = prefer(s + 1, th, curve, y)
Y1 = y(m2) - y(m1)
c1 = curve(m2) - curve(m1)
If m2 > m1 Then Debug.Print y(m1), 1 - c1 / Y1
Next
End Sub