2014-11-05 16 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