よし、楽しかった。
FuncEval というクラスを作成します。
Option Explicit
Dim output_ As Double
Dim input_() As Double
Public Property Get VectArr() As Double()
VectArr = input_
End Property
Public Function Vect(i As Integer)
Vect = input_(i)
End Function
Public Sub SetVect(ByRef newVect() As Double)
Dim i As Integer
ReDim input_(LBound(newVect) To UBound(newVect)) As Double
For i = LBound(newVect) To UBound(newVect)
input_(i) = newVect(i)
Next i
End Sub
Public Property Get Result() As Double
Result = output_
End Property
Public Property Let Result(newRes As Double)
output_ = newRes
End Property
Func というクラス:
Option Explicit
Private cube_ As Double
Public Property Let Cube(newCube As Double)
cube_ = newCube
End Property
Public Function Eval(ByRef val() As Double) As FuncEval
Dim ret As New FuncEval
ret.Result = Abs(cube_ - val(0) * val(0) * val(0))
ret.SetVect val
Set Eval = ret
End Function
このコードを標準モジュールに入れます。
Option Explicit
Function NelderMead(f As Func, _
ByRef guess() As Double) As Double()
'Algorithm follows that outlined here:
'http://www.mathworks.com/help/techdoc/math/bsotu2d.html#bsgpq6p-11
'Used as the perturbation for the initial guess when guess(i) == 0
Dim zeroPert As Double
zeroPert = 0.00025
'The factor each element of guess(i) is multiplied by to obtain the
'initial simplex
Dim pertFact As Double
pertFact = 1.05
'Tolerance
Dim eps As Double
eps = 0.000000000001
Dim shrink As Boolean
Dim i As Integer, j As Integer, n As Integer
Dim simplex() As Variant
Dim origVal As Double, lowest As Double
Dim m() As Double, r() As Double, s() As Double, c() As Double, cc() As Double, diff() As Double
Dim FE As FuncEval, FR As FuncEval, FS As FuncEval, FC As FuncEval, FCC As FuncEval, newFE As FuncEval
n = UBound(guess) - LBound(guess) + 1
ReDim m(LBound(guess) To UBound(guess)) As Double
ReDim r(LBound(guess) To UBound(guess)) As Double
ReDim s(LBound(guess) To UBound(guess)) As Double
ReDim c(LBound(guess) To UBound(guess)) As Double
ReDim cc(LBound(guess) To UBound(guess)) As Double
ReDim diff(LBound(guess) To UBound(guess)) As Double
ReDim simplex(LBound(guess) To UBound(guess) + 1) As Variant
Set simplex(LBound(simplex)) = f.Eval(guess)
'Generate the simplex
For i = LBound(guess) To UBound(guess)
origVal = guess(i)
If origVal = 0 Then
guess(i) = zeroPert
Else
guess(i) = pertFact * origVal
End If
Set simplex(LBound(simplex) + i - LBound(guess) + 1) = f.Eval(guess)
guess(i) = origVal
Next i
'Sort the simplex by f(x)
For i = LBound(simplex) To UBound(simplex) - 1
For j = i + 1 To UBound(simplex)
If simplex(i).Result > simplex(j).Result Then
Set FE = simplex(i)
Set simplex(i) = simplex(j)
Set simplex(j) = FE
End If
Next j
Next i
Do
Set newFE = Nothing
shrink = False
lowest = simplex(LBound(simplex)).Result
'Calculate m
For i = LBound(m) To UBound(m)
m(i) = 0
For j = LBound(simplex) To UBound(simplex) - 1
m(i) = m(i) + simplex(j).Vect(i)
Next j
m(i) = m(i) / n
Next i
'Calculate the reflected point
For i = LBound(r) To UBound(r)
r(i) = 2 * m(i) - simplex(UBound(simplex)).Vect(i)
Next i
Set FR = f.Eval(r)
'Check acceptance conditions
If (simplex(LBound(simplex)).Result <= FR.Result) And (FR.Result < simplex(UBound(simplex) - 1).Result) Then
'Accept r, replace the worst value and iterate
Set newFE = FR
ElseIf FR.Result < simplex(LBound(simplex)).Result Then
'Calculate the expansion point, s
For i = LBound(s) To UBound(s)
s(i) = m(i) + 2 * (m(i) - simplex(UBound(simplex)).Vect(i))
Next i
Set FS = f.Eval(s)
If FS.Result < FR.Result Then
Set newFE = FS
Else
Set newFE = FR
End If
ElseIf FR.Result >= simplex(UBound(simplex) - 1).Result Then
'Perform a contraction between m and the better of x(n+1) and r
If FR.Result < simplex(UBound(simplex)).Result Then
'Contract outside
For i = LBound(c) To UBound(c)
c(i) = m(i) + (r(i) - m(i)) / 2
Next i
Set FC = f.Eval(c)
If FC.Result < FR.Result Then
Set newFE = FC
Else
shrink = True
End If
Else
'Contract inside
For i = LBound(cc) To UBound(cc)
cc(i) = m(i) + (simplex(UBound(simplex)).Vect(i) - m(i)) / 2
Next i
Set FCC = f.Eval(cc)
If FCC.Result < simplex(UBound(simplex)).Result Then
Set newFE = FCC
Else
shrink = True
End If
End If
End If
'Shrink if required
If shrink Then
For i = LBound(simplex) + 1 To UBound(simplex)
For j = LBound(simplex(i).VectArr) To UBound(simplex(i).VectArr)
diff(j) = simplex(LBound(simplex)).Vect(j) + (simplex(i).Vect(j) - simplex(LBound(simplex)).Vect(j)) / 2
Next j
Set simplex(i) = f.Eval(diff)
Next i
End If
'Insert the new element in place
If Not newFE Is Nothing Then
For i = LBound(simplex) To UBound(simplex)
If simplex(i).Result > newFE.Result Then
For j = UBound(simplex) To i + 1 Step -1
Set simplex(j) = simplex(j - 1)
Next j
Set simplex(i) = newFE
Exit For
End If
Next i
End If
Loop Until (simplex(UBound(simplex)).Result - simplex(LBound(simplex)).Result) < eps
NelderMead = simplex(LBound(simplex)).VectArr
End Function
Function test(cube, guess) As Double
Dim f As New Func
Dim guessVec(0 To 0) As Double
Dim Result() As Double
Dim i As Integer
Dim output As String
f.cube = cube
guessVec(0) = guess
Result = NelderMead(f, guessVec)
test = Result(0)
End Function
Func クラスには残差関数が含まれています。NelderMead メソッドは Func クラスの Result メソッドのみを必要とするため、Eval メソッドが最初の推測と同じ長さのベクトルを処理し、FuncEval オブジェクトを返す限り、Func クラスを自由に使用できます。
テスト関数を呼び出して、動作を確認します。注意してください、私は実際に多次元ベクトルでテストしていません。外出する必要があります。問題があればお知らせください!
編集：関数の受け渡しを一般化するための提案
さまざまな問題に対して、多数の異なるクラスを作成する必要があります。つまり、NelderMead 関数を一般的なものにするには、その宣言行を次のように変更する必要があります。
Function NelderMead(f As Object, _
ByRef guess() As Double) As Double()
f が何であれ、double の配列を取る Eval メソッドが常に必要です。
編集：関数の受け渡し、おそらくVBAで行われることを意図した（ばかげた）方法
Function f(x() As Double) As Double
f = x(0) * x(0)
End Function
Sub Test()
Dim x(0 To 0) As Double
x(0) = 5
Debug.Print Application.Run("f", x)
End Sub
このメソッドを使用すると、次の宣言ができます。
Function NelderMead(f As String, _
ByRef guess() As Double) As Double()
次に、上記の Application.Run 構文を使用して f を呼び出します。関数内でもいくつかの変更を加える必要があります。きれいではありませんが、率直に言って、最初はそれほどきれいではありませんでした。