python - scipy 最適化 fmin 構文

Question

numseq = ['0012000', '0112000', '0212000', '0312000', '1012000', '1112000',                                                                                   '1212000', '1312000', '2012000', '2112000', '2212000', '2312000', '3012000', '3112000',          '3212000', '3312000', '0002000', '0022000', '0032000', '1002000', '1022000', '1032000',     '2002000', '2022000', '2032000', '3002000', '3022000', '3032000', '0010000', '0011000', '0013000', '1010000', '1011000', '1013000', '2010000', '2011000', '2013000', '3010000', '3011000', '3013000', '0012100', '0012200', '0012300', '1012100', '1012200', '1012300', '2012100', '2012200', '2012300', '3012100']
prob = [-0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.78361598908750163, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212]

numseqとprobは、それぞれ長さが 50 のリストです。それらは収集された実験データです。numseqX 軸の値にprob対応し、Y 軸の値に対応します。

最小化したい機能は次のとおりです。

def residue(allparams, xdata, ydata):
    chi2 = 0.0
    for i in range(0,len(xdata)):
        x = xdata[i]
        y = 0
        for j in range(len(x)):
            y = y-allparams[int(x[j])][j]
            chi2 = chi2 + (ydata[i]-y)**2
return chi2

そう：

• allparams最適化するすべてのパラメーターを含む 4 × 7 行列です。
• xdataは X 軸の値、つまりnumseq
• ydata単なる数字のリストです。prob

chi2は、実験値とモデル値の差の 2 乗です。これは最小限に抑える必要があります。

パラメーターの初期推定値は次のように与えられます。

x0 = [[-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6]]

fminでは、この関数を呼び出すにはどうすればよいでしょうか。私は試した

fmin(residue, x0, args=(numseq, prob))

しかし、私はエラーが発生し続けます:

Traceback (most recent call last):
  File "<pyshell#362>", line 1, in <module>
    fmin(residue, x0, args=(numseq, prob))
  File "C:\Python31\lib\site-packages\scipy\optimize\optimize.py", line 258, in fmin
    fsim[0] = func(x0)
  File "C:\Python31\lib\site-packages\scipy\optimize\optimize.py", line 177, in function_wrapper
    return function(x, *args)
  File "<pyshell#361>", line 7, in residue
    y = y-allparams[int(x[j])][j]
IndexError: invalid index to scalar variable.

これはなぜですか？fmin2D 配列を最初の推測として受け入れることができないためですか? 次に、1D 配列のパラメーターで動作するようにコード全体を変更する必要がありますか?

この問題を説明できなくても、少なくともfminモジュールがどのように機能するかを正確に教えていただけますか? fminつまり、 N 次元配列の最適化を実装する方法の構文は? とは何か説明していただけますargs()か？私は最適化が初めてで、それを実装する方法がわかりません:(

Answer 1

「fmin」ルーチンは、2 次元配列を最初の推測として受け入れることができます。しかし、最初に行うことは、この配列を平坦化することです [(4,7) --> (28)]。そのため、剰余関数が (4,7) 配列を入力として取り、"fmin" ルーチンが長さ 28 のフラット化された "x0" を与えることになります。そのため、エラーが表示されます:
y = y-allparams[int(x[j])][j]
IndexError: invalid index to scalar variable.

See the source code here.

そのため、配列ではなくベクトルを受け入れるように剰余関数を変更する必要があるようです。ただし、これはそれほど悪くはないようです。私はうまくいくように見えた次のことを試しました（注：再確認してください！）

def residue_alternative(allparams, inshape, xdata, ydata):
    m, n = inshape
    chi2 = 0.0
    for i in range(0,len(xdata)):
        x = xdata[i]
        y = 0
        for j in range(len(x)):
            idx = int(x[j]) * n +  j #Double check this to 
            y = y-allparams[idx]     #make sure it does what you want
            chi2 = chi2 + (ydata[i]-y)**2
    return chi2

私はそれを使用して呼び出しました：

x0 = -0.6 * np.ones((4,7), dtype=np.double)
[xopt, fopt, iter, funcalls, warnflag] = \
    fmin(residue_alternative, x0, args=(x0.shape, numseq, prob),
       maxiter = 100000,
       maxfun  = 100000,
       full_output=True,
       disp=True)

そして、次の結果を受け取りました。

Optimization terminated successfully.
         Current function value: 7.750523
         Iterations: 21570
         Function evaluations: 26076

>>>xopt
array([ 0.57669042, -0.21965861,  0.2635061 , -0.08284016, -0.0779489 ,
   -0.10358114,  0.14041582,  0.72469391, -0.43190214,  0.31269757,
   -0.0338726 , -0.14919739, -2.58314651,  2.74251214,  0.57695759,
   -0.49574628,  0.1490926 ,  0.04912353,  0.02420988,  1.17924051,
   -7.2147027 ,  0.57860843, -0.28386938,  0.2431877 , -0.22674694,
   -0.58308225, -6.05706775, -2.06350063])    

4x7 配列に変形できます。これを試してみて、うまくいくかどうか教えてください。