無法適應函數scipy.optimize.curve_fit（）

Question

我想，以適應以下功能：

def invlaplace_stehfest2(time,EL,tau): 
    upsilon=0.25 
    pmax=6.9 
    E0=0.0154 
    M=8 
    results=[] 
    for t in time: 
     func=0 
     for k in range(1,2*M+1): 
      SUM=0 
      for j in range(int(math.floor((k+1)/2)),min(k,M)+1): 
       dummy=j**(M+1)*scipy.special.binom(M,j)*scipy.special.binom(2*j,j)*scipy.special.binom(j,k-j)/scipy.math.factorial(M) 
       SUM+=dummy 
      s=k*math.log(2)/t[enter image description here][1] 
      func+=(-1)**(M+k)*SUM*pmax*EL/(mp.exp(tau*s)*mp.expint(1,tau*s)*E0+EL)/s 

     func=func*math.log(2)/t 
     results.append(func) 
    return [float(i) for i in results] 

這樣做我用以下數據：

data_time=np.array([69.0,99.0,139.0,179.0,219.0,259.0,295.5,299.03]) 
data_relax=np.array([6.2536,6.1652,6.0844,6.0253,5.9782,5.9404,5.9104,5.9066]) 

隨着如下因素的猜測：

guess=np.array([0.1,0.05]) 

而且scipy.optimize.curve_fit（）爲便接踵而來：

Parameter,Covariance=scipy.optimize.curve_fit(invlaplace_stehfest2,data_time,data_relax,guess) 

由於我不明白的原因，我無法正確地擬合數據。我看到下面的圖。

Bad fitting

我的功能，因爲當我使用正確的猜測無疑是正確的：

guess=np.array([0.33226685047281592707364253044085038793404361200072,8.6682623502960394383501102909774397295654841654769]) 

我能夠正確地適合我的數據。

Expected fitting

爲什麼我不能夠正確地安裝任何提示？我應該使用另一種方法嗎？

這裏是空穴程序：

############################################## 
import matplotlib 
import matplotlib.pyplot as plt 
import matplotlib.mlab as mlab 
import matplotlib.pylab as mplab 
import math 
from math import * 
import numpy as np 
import scipy 
from scipy.optimize import curve_fit 
import mpmath as mp 

############################################################################################ 
def invlaplace_stehfest2(time,EL,tau): 
    upsilon=0.25 
    pmax=6.9 
    E0=0.0154 
    M=8 
    results=[] 
    for t in time: 
     func=0 
     for k in range(1,2*M+1): 
      SUM=0 
      for j in range(int(math.floor((k+1)/2)),min(k,M)+1): 
       dummy=j**(M+1)*scipy.special.binom(M,j)*scipy.special.binom(2*j,j)*scipy.special.binom(j,k-j)/scipy.math.factorial(M) 
       SUM+=dummy 
      s=k*math.log(2)/t 
      func+=(-1)**(M+k)*SUM*pmax*EL/(mp.exp(tau*s)*mp.expint(1,tau*s)*E0+EL)/s 

     func=func*math.log(2)/t 
     results.append(func) 
    return [float(i) for i in results] 


############################################################################################  

###Constant### 

####Value#### 
data_time=np.array([69.0,99.0,139.0,179.0,219.0,259.0,295.5,299.03]) 
data_relax=np.array([6.2536,6.1652,6.0844,6.0253,5.9782,5.9404,5.9104,5.9066]) 


###Fitting### 
guess=np.array([0.33226685047281592707364253044085038793404361200072,8.6682623502960394383501102909774397295654841654769]) 
#guess=np.array([0.1,0.05]) 
Parameter,Covariance=scipy.optimize.curve_fit(invlaplace_stehfest2,data_time,data_relax,guess) 
print Parameter 
residu=sum(data_relax-invlaplace_stehfest2(data_time,Parameter[0],Parameter[1])) 


Graph_Curves=plt.figure() 
ax = Graph_Curves.add_subplot(111) 
ax.plot(data_time,invlaplace_stehfest2(data_time,Parameter[0],Parameter[1]),"-") 
ax.plot(data_time,data_relax,"o") 
plt.show() 

Answer 1

非線性鉗工如默認列文伯格 - 馬夸爾特在scipy.optimize.curve_fit（解算器使用），最喜歡的迭代求解器，可以在一個局部最小值停止在錯誤空間。如果錯誤空間是「顛簸」，那麼初始參數估計就變得非常重要，就像在這種情況下一樣。

Scipy已將差分進化遺傳算法添加到優化模塊，該模塊可用於確定curve_fit（）的初始參數估計值。 Scipy的實現使用拉丁超立方體算法來確保對參數空間的徹底搜索，要求參數上下限在其中進行搜索。正如你在下面看到的，我已經使用這個scipy模塊來替換代碼中名爲「guess」的值的硬編碼值。這不會很快運行，但是一個較慢的正確結果要比一個快速錯誤的結果好得多。試試這段代碼：

import matplotlib 
import matplotlib.pyplot as plt 
import matplotlib.mlab as mlab 
import matplotlib.pylab as mplab 
import math 
from math import * 
import numpy as np 
import scipy 
from scipy.optimize import curve_fit 

from scipy.optimize import differential_evolution 

import mpmath as mp 

############################################################################################ 
def invlaplace_stehfest2(time,EL,tau): 
    upsilon=0.25 
    pmax=6.9 
    E0=0.0154 
    M=8 
    results=[] 
    for t in time: 
     func=0 
     for k in range(1,2*M+1): 
      SUM=0 
      for j in range(int(math.floor((k+1)/2)),min(k,M)+1): 
       dummy=j**(M+1)*scipy.special.binom(M,j)*scipy.special.binom(2*j,j)*scipy.special.binom(j,k-j)/scipy.math.factorial(M) 
       SUM+=dummy 
      s=k*math.log(2)/t 
      func+=(-1)**(M+k)*SUM*pmax*EL/(mp.exp(tau*s)*mp.expint(1,tau*s)*E0+EL)/s 

     func=func*math.log(2)/t 
     results.append(func) 
    return [float(i) for i in results] 


############################################################################################  

###Constant### 

####Value#### 
data_time=np.array([69.0,99.0,139.0,179.0,219.0,259.0,295.5,299.03]) 
data_relax=np.array([6.2536,6.1652,6.0844,6.0253,5.9782,5.9404,5.9104,5.9066]) 

# function for genetic algorithm to minimize (sum of squared error) 
def sumOfSquaredError(parameterTuple): 
    return np.sum((data_relax - invlaplace_stehfest2(data_time, *parameterTuple)) ** 2) 

###Fitting### 
#guess=np.array([0.33226685047281592707364253044085038793404361200072,8.6682623502960394383501102909774397295654841654769]) 
#guess=np.array([0.1,0.05]) 

parameterBounds = [[0.0, 1.0], [0.0, 10.0]] 
# "seed" the numpy random number generator for repeatable results 
# note the ".x" here to return only the parameter estimates in this example 
guess = differential_evolution(sumOfSquaredError, parameterBounds, seed=3).x 

Parameter,Covariance=scipy.optimize.curve_fit(invlaplace_stehfest2,data_time,data_relax,guess) 
print Parameter 
residu=sum(data_relax-invlaplace_stehfest2(data_time,Parameter[0],Parameter[1])) 


Graph_Curves=plt.figure() 
ax = Graph_Curves.add_subplot(111) 
ax.plot(data_time,invlaplace_stehfest2(data_time,Parameter[0],Parameter[1]),"-") 
ax.plot(data_time,data_relax,"o") 
plt.show()