使用 Nelder-Mead 算法的 Python 最小化函数

Question

我正在尝试使用 Nelder-Mead 算法最小化函数 mymodel 以适应我的数据。这是通过 scipy 的 optimize.fmin 在 myfit 函数中完成的。我想我已经很接近了，但我一定错过了一些东西，因为我不断收到错误:

“操作数无法与形状 (80,) (5,) 一起广播”。

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy import special

def mymodel(c,t,y):

    """
    This function specifies the form to be minimized by fmins in myfit.
    c is a 1 x 5 array containing the fit parameters.    
    """

    m = (np.sin(np.exp(-c[1]*t)*c[0]/2.0))**2

    # compute complete elliptic integral of the first kind with ellipk
    w = np.pi*c[2]/2.0/special.ellipk(m)

    dt = t[1] - t[0]    
    phase = np.cumsum(w)*dt

    z = np.sum((y - c[0] * np.exp(-c[1]*t) * np.cos(phase+c[3])-c[4])**2)

    return z


def myfit(c, pos):

    """
    Fitting procedure for the amplitude decay of the undriven pendulum

    initial fit parameters:
    c[0]=theta_m, c[1]=alpha, c[2]=omega_0, c[3]=phi, c[4]=const.
    pos = the position data
    """    

    # convert data to seconds
    t = 0.001*np.arange(0,len(pos))

    dt = t[1] - t[0]

    # Minimise the function mymodel using Nelder-Mead algorithm
    c = optimize.fmin(mymodel, c, args=(t,y), maxiter=5000, full_output=True)   

    m = (np.sin(np.exp(-c[1]*t)*c[0]/2.0))**2

    # change of frequency with amplitude
    w = np.pi*c[2]/2.0/special.ellipk(m)
    phase = np.cumsum(w)*dt

    # use values from fmin
    fit = c[0]*np.exp(-c[1]*t)*np.cos(phase+c[3])+c[4]

    return t, c, fit

t = np.array([    0.,    15.,    30.,    45.,    60.,    75.,    90.,   105.,
         120.,   135.,   150.,   165.,   180.,   195.,   210.,   225.,
         240.,   255.,   270.,   285.,   300.,   315.,   330.,   345.,
         360.,   375.,   390.,   405.,   420.,   435.,   450.,   465.,
         480.,   495.,   510.,   525.,   540.,   555.,   570.,   585.,
         600.,   615.,   630.,   645.,   660.,   675.,   690.,   705.,
         720.,   735.,   750.,   765.,   780.,   795.,   810.,   825.,
         840.,   855.,   870.,   885.,   900.,   915.,   930.,   945.,
         960.,   975.,  1005.,  1020.,  1035.,  1050.,  1065.,  1080.,
        1095.,  1110.,  1125.,  1140.,  1155.,  1170.,  1185.,  1200.,
        ])

pos = np.array([ 28.95,  28.6 ,  28.1 ,  27.5 ,  26.75,  25.92,  24.78,  23.68,
        22.5 ,  21.35,  20.25,  19.05,  17.97,  16.95,  15.95,  15.1 ,
        14.45,  13.77,  13.3 ,  13.  ,  12.85,  12.82,  12.94,  13.2 ,
        13.6 ,  14.05,  14.65,  15.45,  16.1 ,  16.9 ,  17.75,  18.7 ,
        19.45,  20.3 ,  21.1 ,  21.9 ,  22.6 ,  23.25,  23.75,  24.2 ,
        24.5 ,  24.75,  24.88,  24.9 ,  24.8 ,  24.65,  24.35,  23.9 ,
        23.55,  22.95,  22.5 ,  21.98,  21.3 ,  20.65,  20.05,  19.4 ,
        18.85,  18.3 ,  17.8 ,  17.35,  16.95,  16.6 ,  16.35,  16.2 ,
        16.1 ,  16.1 ,  16.35,  16.5 ,  16.75,  17.02,  17.4 ,  17.75,
        18.3 ,  18.65,  19.1 ,  19.55,  20.  ,  20.45,  20.85,  21.25,
        ])

# fitting with myfit function
c = np.array([1,1,1,1,1]) # initial guess
t, c, fit = myfit(c, pos)

plt.plot(t,fit)
plt.show()

Answer 1

问题在于，使用 full_output=true 参数调用的 fmin 不仅返回优化后的参数，还返回包含参数和其他内容的元组。

参数是元组中的第一个值。您可以使用 [0] 选择此选项:

c = optimize.fmin(mymodel, c, args=(t,pos), maxiter=5000, full_output=True)[0]

或者只需删除 full_output=true 参数即可:

c = optimize.fmin(mymodel, c, args=(t,pos), maxiter=5000)

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.fmin.html