2014-11-05 96 views
1

我试图最小化功能mymodel与内尔德 - 米德算法,以适应我的数据。这是在myfit功能与SciPy的的optimize.fmin完成。我觉得我很接近,但我一定是失去了一些东西,因为我不断收到一个错误:Python的最小化功能与内尔德 - 米德算法

“操作数无法与形状(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() 
+2

包含在问题的完整回溯。这将让我们知道哪条线路导致错误。 – 2014-11-06 02:29:04

+0

你从来没有定义'y'。 – dbliss 2016-03-10 02:46:28

回答

2

的问题是,fmin,称为与full_output=true参数,不仅返回优化的参数,但包含的参数和其他的东西的元组。

的参数是在元组中的第一个值。你可以只用[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