2016-08-22 64 views
0

我有以下x,y數據(綠色)。我想獲得適合我曲線的多項式函數。 Python中適合的曲線看起來很好(藍色)。 當我使用多項式的係數,我自己建立函數時,結果不在藍色曲線上。對於X的小值,這可能仍然適合,但對於大值是完全錯誤的。在圖像中,顯示了x = 15和2.5的y(大點)。python多項式曲線擬合 - 係數不正確

enter image description here

數據:

x, y 
0.5883596178 18562.5 
0.6656014904 20850 
0.7407008741 22700 
0.8310800498 24525 
0.9479506185 26370 
1.0768193651 27922 
1.1983161945 29070 
1.3837939534 30410 
1.6650549531 31800 
1.946640319  32740 
2.3811442965 33655 
2.9126326549 34290 
3.6970654824 34800 
4.2868951065 34987.5 
4.8297935972 35102 
5.7876198835 35175 
7.3463468386 35050 
8.9861037519 34725 
10.5490727095 34285 
13.2260016159 33450 
16.5822270413 32795 
20.5352502646 32472 
25.7462680049 32475 

代碼:

data = plb.loadtxt('fig3_1_tiltingRL.dat') 

x = data[:,0] 
y= data[:,1] 
#plt.xscale('log')#plt.set_xscale('log') 
coefs = poly.polyfit(x, y, 10) 
ffit = poly.polyval(x, coefs) 
plt.plot(x, ffit) 
plt.plot(x, y, 'o') 
print(coefs) 
xPoints =15. 
yPt = (-6.98662492e+03 * xPoints**0 + 6.57987934e+04 * xPoints**1 -\ 
     4.65689536e+04 * xPoints**2 + 1.85406629e+04 * xPoints**3 -\ 
     4.49987278e+03 * xPoints**4 + 6.92952944e+02 * xPoints**5 -\ 
     6.87501257e+01 * xPoints**6 + 4.35851202e+00 * xPoints**7 -\ 
     1.69771617e-01 * xPoints**8 + 3.68535224e-03 * xPoints**9 -\ 
     3.39940049e-05 * xPoints**10) 
print(yPt) 
plt.plot(xPoints, yPt , 'or',label="test" ,markersize=18, color='black') 
plt.show() 
+0

請發表一個簡短而自成體系'.py'源文件:http://sscce.org/ – pts

+0

這是一種可疑的是,即使係數均爲負,奇怪的人都是積極的。 'print(coefs)'的確切輸出是什麼? –

回答

0

在我看來,您使用的是poyval不看我的權利的方式。嘗試使用numpy.linspace生成X軸,然後在其上應用polyval。 類似下面的代碼。

import numpy as np 
import matplotlib.pyplot as plt 


data = np.loadtxt('fig3_1_tiltingRL.dat') 

x = data[:,0] 
y= data[:,1] 
#plt.xscale('log')#plt.set_xscale('log') 
coefs = np.polyfit(x, y, 10) 
ffit = np.polyval(coefs, x) 

new_x = np.linspace(0,26) 
new_ffit = np.polyval(coefs, new_x) 

plt.plot(x, y, 'o', label="Raw") 
plt.plot(x, ffit,'x',label="Fit to Raw") 
plt.plot(new_x, new_ffit,label="Fit to LinSpace") 

# This is ugly. I'd use list comprehension here! 
arr = np.linspace(0,26,20) 
new_y = [] 
for xi in arr: 
    total = 0 
    for i,v in enumerate(coefs[::-1]): 
     total += v*xi**i 
    new_y.append(total) 


plt.plot(arr, new_y, '*', label="Polynomial") 

plt.legend(loc=2) 
plt.show() 

enter image description here

正如你所看到的,有沒有出現在你的情節駝峯......

0

你的算法似乎是工作的罰款。你應該剛剛相反:

coefs = poly.polyfit(x, y, 10) 
ffit = poly.polyval(x, coefs) 

此:

coefs = poly.polyfit(x, y, 10) # fit data to the polynomial 
new_x = np.linspace(0, 30, 50) # new x values to evaluate 
ffit = poly.polyval(new_x, coefs) # fitted polynomial evaluated with new data 

因此,函數poly.polyval將評估new_x代替x的所有點座標,你已經知道了。

0

非常感謝您回答我的問題。

由silgon和RicLeal提供的解決方案都可以工作。

最後,由於我有幾條曲線,我已經應用了RicLeal給出的解決方案。

我的數據記錄在X軸上。我只是修改了RicLeal給出的代碼,我對結果感到滿意。

enter image description here

x = data[:,0] 
y= data[:,1] 
plt.xscale('log')#plt.set_xscale('log') 
logx=np.log10(x) 
coefs = np.polyfit(logx, y, 10) 
ffit = np.polyval(coefs, logx) 
print (coefs) 


logxmin=math.log10(0.5883596178) 
logxmax=math.log10(26.) 

new_x = np.logspace(logxmin, logxmax,50) 
lognew_x=np.log10(new_x) 
new_ffit = np.polyval(coefs, lognew_x) 
plt.semilogx(x, y, 'o', label="Raw") 
plt.semilogx(x, ffit,'x',label="Fit to Raw") 
plt.semilogx(new_x, new_ffit,label="Fit to LogSpace") 
print(lognew_x, new_ffit) 
# This is ugly. I'd use list comprehension here! 
arr = np.logspace(logxmin, logxmax,50) 
arrlog= np.log10(arr) 
new_y = [] 
for xi in arrlog: 
    total = 0 
    for i,v in enumerate(coefs[::-1]): 
     #print (v) 
     total += v*xi**i 
    new_y.append(total) 


    plt.semilogx(arr, new_y, '*', label="Polynomial") 

coeffs= [6.85869364,  -92.86678553, 343.39375022,   -555.52532934, 434.18179364, 
    -152.82724751,  9.71300951, 21.68653301, -35.62838377, 28.3985976, 
     27.04762122] 
new_testy = [] 
for xi in arrlog: 
total = 0 
    for i,v in enumerate(coeffs[::-1]): 
     #print (v) 
     total += v*xi**i 
    new_testy.append(total)   
plt.semilogx(arr, new_testy, 'o', label="Polynomial")   
plt.legend(loc=2) 
plt.show()