2017-02-11 13 views
2

這是我第一次發佈一個問題,我會盡量讓它儘可能清晰,但可以隨時提問。當其中一個擬合參數更改xdata輸入數組值時,我可以使用適合python的scipy.curve嗎?

我嘗試使用下面的方法scipy.curve_fit以適應模型曲線:

import numpy as np 
import matplotlib.pyplot as pyplot 
import scipy 
from scipy.optimize import curve_fit 




def func2(x,EM): 
    return (((4.0*EM*(np.sqrt(8*10**-9)))/(3.0*(1.0-(0.5**2))*8*10**-9))*(((((x))*1*10**-9)**((3.0/2.0))))) 


ydata=[-0.003428768, -0.009050058, -0.0037997673999999996, -0.0003833233, -0.007557649, -0.0034860994, -0.0009856887, -0.0017508664, -0.00036931394999999996, 
     -0.0040713947, -0.005737315000000001, 0.0005120568, -0.007336486, -0.00719302, -0.0039941817, -0.0029785274, -0.0013044578, -0.008190335, -0.00833507, 
     -0.0074282060000000006, -0.009629990000000001, -0.009425125, -0.008662485999999999, -0.0019445216, -0.008331748, -0.009513038, -0.0047609017, -0.004364422, 
     -0.010325097, -0.0036570733, -0.0060091914, -0.005655772, -0.0045517069999999995, -0.00066998035, 0.006374902, 0.006445733, 0.0019101816, 
     0.010262737999999999, 0.011139007, 0.018161469, 0.016963122, 0.022915895, 0.027177791, 0.028707139, 0.040105638, 0.044088004, 0.041657403, 
     0.052325636999999994, 0.062399405, 0.07020844, 0.076979915, 0.08888523, 0.099634745, 0.10961602, 0.12188646, 0.13677225, 0.15639512, 0.16833586, 
     0.18849944000000002, 0.21515548, 0.23989769000000002, 0.26319308, 0.29388397, 0.321042, 0.35637776, 0.38564656999999997, 0.4185209, 0.44986692, 
     0.48931552999999994, 0.52583893, 0.5626885, 0.6051665, 0.6461075, 0.69644346, 0.7447817, 0.7931281, 0.8381386000000001, 0.8883482, 0.9395609999999999, 
     0.9853629, 1.0377034, 1.0889026, 1.1334094] 


xdata=[34.51388, 33.963736999999995, 
     33.510695, 33.04127, 32.477253, 32.013624, 31.536019999999997, 31.02925, 30.541649999999997, 
     30.008646, 29.493828, 29.049707, 28.479668, 27.980956, 27.509590000000003, 27.018721, 26.533737, 25.972296, 
     25.471065, 24.979228000000003, 24.459624, 23.961517, 23.46839, 23.028454, 22.471411, 21.960924, 21.503428000000003, 
     21.007033, 20.453855, 20.013475, 19.492528, 18.995746999999998, 18.505670000000002, 18.040403, 17.603387, 17.104082, 
     16.563634, 16.138298000000002, 15.646187, 15.20897, 14.69833, 14.25156, 13.789688, 13.303409, 12.905278, 12.440909, 11.919262, 
     11.514609, 11.104646, 10.674512, 10.235055, 9.84145, 9.437523, 9.026733, 8.63639, 8.2694065, 7.944733, 7.551445, 7.231599999999999, 
     6.9697434, 6.690793299999999, 6.3989780000000005, 6.173159, 5.9157856, 5.731453, 5.4929328, 5.2866156, 5.066648000000001, 4.9190496, 
     4.745381399999999, 4.574569599999999, 4.4540283, 4.3197597000000005, 4.2694026, 4.2012034, 4.133134, 4.035212, 3.9837262, 3.9412007, 3.8503475999999996, 
     3.8178950000000005, 3.7753053999999997, 3.6728842] 


dstart=20.0 

xdata=np.array(xdata[::-1]) 
xdata=xdata-dstart 
xdata=list(xdata) 

xdata1=[] 
ydata1=[] 
for i in range(len(xdata)): 
    if xdata[i]>0: 
     xdata1.append(xdata[i]) 
     ydata1.append(ydata[i]) 

xdata=np.array(xdata1) 
ydata=np.array(ydata1) 

popt, pcov = curve_fit(func2, xdata, ydata) 
a=popt[0] 

print "E=", popt[0]/10**6 


t=func2(xdata,a) 

ax=pyplot.figure().add_subplot(1,1,1) 
ax.plot(xdata,t, color="blue",mew=2.0,label="Hertz Fit") 
ax.plot(xdata,ydata,ls="",marker="x",color="red",mew=2.0,label="Data") 
ax.legend(loc=2) 
pyplot.show() 

的「DSTART」值基本上切斷了代碼,我不想下部以適應,因爲它是負面的,模型不喜歡負數。目前,我必須在運行代碼之前手動設置「dstart」,然後才能看到最終結果。

我開始在Excel中使用Solver進行擬合,通過嵌套用「dstart」調整xdata的代碼來同時改變「EM」變量和「dstart」變量,並將負值切入函數適合。

基本上我要的是:

import numpy as np 
import matplotlib.pyplot as pyplot 
import scipy 
from scipy.optimize import curve_fit 




def func2(x,EM,dstart): 

    xdata=np.array(x[::-1]) 
    xdata=dstart-xdata 
    xdata=list(xdata) 

    xdata1=[] 
    for i in range(len(xdata)): 
     if xdata[i]>0: 
      xdata1.append(xdata[i]) 

    global xdata2 
    xdata2=np.array(xdata1) 






    return (((4.0*EM*(np.sqrt(8*10**-9)))/(3.0*(1.0-(0.5**2))*8*10**-9))*(((((xdata2))*1*10**-9)**((3.0/2.0))))) 


ydata=[-0.003428768, -0.009050058, -0.0037997673999999996, -0.0003833233, -0.007557649, -0.0034860994, -0.0009856887, -0.0017508664, -0.00036931394999999996, 
     -0.0040713947, -0.005737315000000001, 0.0005120568, -0.007336486, -0.00719302, -0.0039941817, -0.0029785274, -0.0013044578, -0.008190335, -0.00833507, 
     -0.0074282060000000006, -0.009629990000000001, -0.009425125, -0.008662485999999999, -0.0019445216, -0.008331748, -0.009513038, -0.0047609017, -0.004364422, 
     -0.010325097, -0.0036570733, -0.0060091914, -0.005655772, -0.0045517069999999995, -0.00066998035, 0.006374902, 0.006445733, 0.0019101816, 
     0.010262737999999999, 0.011139007, 0.018161469, 0.016963122, 0.022915895, 0.027177791, 0.028707139, 0.040105638, 0.044088004, 0.041657403, 
     0.052325636999999994, 0.062399405, 0.07020844, 0.076979915, 0.08888523, 0.099634745, 0.10961602, 0.12188646, 0.13677225, 0.15639512, 0.16833586, 
     0.18849944000000002, 0.21515548, 0.23989769000000002, 0.26319308, 0.29388397, 0.321042, 0.35637776, 0.38564656999999997, 0.4185209, 0.44986692, 
     0.48931552999999994, 0.52583893, 0.5626885, 0.6051665, 0.6461075, 0.69644346, 0.7447817, 0.7931281, 0.8381386000000001, 0.8883482, 0.9395609999999999, 
     0.9853629, 1.0377034, 1.0889026, 1.1334094] 


xdata=[34.51388, 33.963736999999995, 
     33.510695, 33.04127, 32.477253, 32.013624, 31.536019999999997, 31.02925, 30.541649999999997, 
     30.008646, 29.493828, 29.049707, 28.479668, 27.980956, 27.509590000000003, 27.018721, 26.533737, 25.972296, 
     25.471065, 24.979228000000003, 24.459624, 23.961517, 23.46839, 23.028454, 22.471411, 21.960924, 21.503428000000003, 
     21.007033, 20.453855, 20.013475, 19.492528, 18.995746999999998, 18.505670000000002, 18.040403, 17.603387, 17.104082, 
     16.563634, 16.138298000000002, 15.646187, 15.20897, 14.69833, 14.25156, 13.789688, 13.303409, 12.905278, 12.440909, 11.919262, 
     11.514609, 11.104646, 10.674512, 10.235055, 9.84145, 9.437523, 9.026733, 8.63639, 8.2694065, 7.944733, 7.551445, 7.231599999999999, 
     6.9697434, 6.690793299999999, 6.3989780000000005, 6.173159, 5.9157856, 5.731453, 5.4929328, 5.2866156, 5.066648000000001, 4.9190496, 
     4.745381399999999, 4.574569599999999, 4.4540283, 4.3197597000000005, 4.2694026, 4.2012034, 4.133134, 4.035212, 3.9837262, 3.9412007, 3.8503475999999996, 
     3.8178950000000005, 3.7753053999999997, 3.6728842] 

xdata2=list(xdata2) 
ydata1=[] 
for i in range(len(xdata2)): 
    if xdata2[i]>0: 
     ydata1.append(ydata[i]) 




popt, pcov = curve_fit(func2, xdata, ydata) 

但是,當我得到一個值錯誤,這並不工作「ValueError異常:操作數不能連同形狀(28)廣播(30)」。我認爲我需要的是curve_fit引入xdata,通過第一個猜測的「dstart」進行調整,猜測EM並檢查是否合適並且最小化了錯誤,嘗試新的「dstart」來調整xdata,猜測EM並檢查是否合適並最小化錯誤,等等。由於我對Python還是比較新的,所以我肯定沒有使用曲線擬合的元素,如果我沒有可能運行數千條曲線,我只會使用Excel。

任何幫助,將不勝感激!

回答

0

我會一分爲二,這樣的:概念和編碼相關

概念:

讓我們改寫你的問題開始。現在的答案是:顯然,是的。只需吸收目標函數中x的參數相關變化。但那不會解決你的問題。你真的似乎感興趣的是如何處理一些x不能被你的函數處理的參數。沒有一個通用的解決方案。

您可以選擇將這些參數視爲不可接受,在這種情況下,您不得不求助於約束優化。 scipy中有幾個解決方案可以做到這一點。

您可以選擇在擬合前從數據集中移除難點。

你可以引入軟約束和壞值,而不是完全排除它們。

編程風格:

for數字程序中的循環。還有在這個網站上的職位gazillions,所以我只舉一個例子:

xdata2=list(xdata2) 
ydata1=[] 
for i in range(len(xdata2)): 
    if xdata2[i]>0: 
     ydata1.append(ydata[i]) 

可以在一行中寫入,將執行速度更快,並返回的array而不是list

ydata1 = ydata[xdata2 > 0] 

看看numpy教程/文檔或搜索此網站的「矢量化」,如果你想學習這種技術。

除此之外,沒有投訴。

爲什麼你的第二個程序不起作用。

您正在篩選您的x和您的y,因此它們應具有相同的形狀。但是,你繼續使用舊的副本,而不是新的y,而你使用新的x。這就是爲什麼形狀不匹配

Btw。你設置它的方式(修改xfunc2之內)或多或少地實現了前面提到的吸收策略。只有,由於您無法訪問y,因此無法更改x的形狀。

相關問題