评估 numpy 多项式 returns 错误的值
Evaluating numpy polynomial returns the wrong values
我有一个系数为 c=2.44717404e-03 和 m=1.88697661e+01 的线性 numpy 多项式,但是当我尝试对其求值时,它给出了错误的输出。
f=np.polynomial.polynomial.Polynomial([2.44717404e-03, 1.88697661e+01])
例如
f(0.015)
returns
3.9646796502044523
而不是正确的结果
0.28549366554
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 13 17:15:51 2020
@author: tomhe
"""
#Modules
if True:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
import scipy.optimize as scpo
#Functions
def M(params, *args):
c=params[0]
m=params[1]
x=args[0]
y=args[1]
ex=args[2]
ey=args[3]
return(sum((y-(m*x+c))**2/(ey**2+m**2*ex**2)))
#Data
if True:
df=pd.read_excel("Biot Savart Experiment Data 2.xlsx")
B=np.array(df['B'].values)
x=np.arange(0,18,0.5)
a=3.95
x1star=6.35
x2star=10.3
I=20.02
eB=np.array([0.001]*len(B))
ex=0.05
ea=0.1
ex1star=0.05
ex2star=0.05
eI=0.01
#Linearisation
if True:
z=((x-x1star)**2+a**2)**(-3/2)+((x-x2star)**2+a**2)**(-3/2)
ez=3*np.sqrt(((x-x1star)**2+a**2)**(-5)*((x-x1star)**2*(ex**2+ex1star**2)+(a*ea)**2)+((x-x2star)**2+a**2)**(-5)*((x-x2star)**2*(ex**2+ex2star**2)+(a*ea)**2))
#Linear regression
if True:
optrslt=scpo.minimize(M,[0,20],(z,B,ez,eB))
optparams=optrslt.x
f=np.polynomial.polynomial.Polynomial(optparams,domain=(min(z),max(z)))
linfitcoords=f.linspace(1000)
#Resiudals
if True:
residuals=f(z)-B
#Residual level line
if True:
rf=np.polynomial.polynomial.Polynomial([0],domain=(min(z),max(z)))
rfcoords=rf.linspace(1000)
#Plotting
if True:
fig=plt.figure(1)
frame=plt.gca()
axis1=fig.add_axes((0,0,1,1))
axis1.errorbar(z,B,xerr=ez,fmt='o', ms=1, lw=0.5, color='red', capsize=2, capthick=0.5, ecolor='black')
#axis1.plot(z,line(z,optparams[1],optparams[0]),lw=0.5,color='blue')
#axis1.plot(x,y,lw=0.5,color='blue')
axis1.plot(linfitcoords[0],linfitcoords[1],lw=0.5,color='blue')
plt.ylabel('Magnetic Field Strength ($mT$)')
axis1.xaxis.set_visible(False)
axis1.set_xlim(left=0,right=0.025)
axis1.set_ylim(bottom=0,top=0.5)
axis2=fig.add_axes((0,-0.4,1,0.4))
axis2.set_xlim(left=0,right=0.025)
axis2.plot(rfcoords[0],rfcoords[1],lw=0.5, color='blue')
axis2.scatter(z,residuals, s=1, color='red')
plt.ylabel("Residuals ($mT)$")
axis1.yaxis.set_label_coords(-0.1,0.5)
axis2.yaxis.set_label_coords(-0.1,0.5)
plt.xlabel('z $(cm^{-3})$')
plt.xlim([0,0.025])
axis3=fig.add_axes((1,-0.4,0.2,0.4))
axis3.hist(residuals,orientation='horizontal',color='gray')
axis3.yaxis.set_visible(False)
axis3.xaxis.tick_top()
axis3.xaxis.set_label_position('top')
plt.xlabel('Frequency')
plt.xticks([1,3,5,7])
#Note:restarting kernel (close window) fixes x axis label problem where you accidentally set as variable and then it breaks
#Save Figure
if True:
plt.savefig("Linearised_Biot_Savart_Graph.png",dpi=1000,bbox_inches='tight')
#Gradient and mu_0 analysis
if True:
m=optparams[1]
mu0=2*m/(I*a**2)
em=1#need to actually figure this out
emu0=mu0*np.sqrt((em/m)**2+(eI/I)**2+4*(ea/a)**2)
# print(mu0)
#print(emu0)
#Normality of residuals
if True:
print(stats.shapiro(residuals))
print(stats.normaltest(residuals))
#print(stats.kstest(residuals,stats.norm(np.mean(residuals),np.var(residuals)).cdf))
注意:所有的 if 语句只是为了方便 Spyder 中的代码折叠。
对于上下文,我正在对磁场的一些物理数据进行一些数据分析和可视化。我的主要问题是我试图绘制一条最适合数据的线,即使我得到的参数对于最适合的线是正确的,当我尝试使用 numpy 多项式生成坐标时绘制最佳拟合线,它生成错误的坐标。
问题似乎出在传递给 f
的 domain
参数上
从 docs 中,您正在为域传递值,但没有为 window 传递值,因此 f
根据地图 [min(z), max(z)] -> [-1, 1]
[=18 映射输入=]
这就是输出不同的原因。
要清楚地看到这一点,请看下面的示例
import numpy as np
# x**2 with no domain shift
f = np.polynomial.polynomial.Polynomial([0,0,1])
print(f)
>> poly([0. 0. 1.])
print(f(1))
>> 1.0
# x**2 with domain shift [-0.1, 0.1] -> [-1, 1], so an input of 1 maps to 10
g = np.polynomial.polynomial.Polynomial([0,0,1], domain=[-0.1, 0.1])
print(g)
>> poly([0. 0. 1.])
print(g(1))
>> 100.0
因此,任何不是 [-1,1]
的 domain
值都会给出一个多项式,由于输入
的缩放,returns 输出 'wrong'
我有一个系数为 c=2.44717404e-03 和 m=1.88697661e+01 的线性 numpy 多项式,但是当我尝试对其求值时,它给出了错误的输出。
f=np.polynomial.polynomial.Polynomial([2.44717404e-03, 1.88697661e+01])
例如
f(0.015)
returns
3.9646796502044523
而不是正确的结果
0.28549366554
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 13 17:15:51 2020
@author: tomhe
"""
#Modules
if True:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
import scipy.optimize as scpo
#Functions
def M(params, *args):
c=params[0]
m=params[1]
x=args[0]
y=args[1]
ex=args[2]
ey=args[3]
return(sum((y-(m*x+c))**2/(ey**2+m**2*ex**2)))
#Data
if True:
df=pd.read_excel("Biot Savart Experiment Data 2.xlsx")
B=np.array(df['B'].values)
x=np.arange(0,18,0.5)
a=3.95
x1star=6.35
x2star=10.3
I=20.02
eB=np.array([0.001]*len(B))
ex=0.05
ea=0.1
ex1star=0.05
ex2star=0.05
eI=0.01
#Linearisation
if True:
z=((x-x1star)**2+a**2)**(-3/2)+((x-x2star)**2+a**2)**(-3/2)
ez=3*np.sqrt(((x-x1star)**2+a**2)**(-5)*((x-x1star)**2*(ex**2+ex1star**2)+(a*ea)**2)+((x-x2star)**2+a**2)**(-5)*((x-x2star)**2*(ex**2+ex2star**2)+(a*ea)**2))
#Linear regression
if True:
optrslt=scpo.minimize(M,[0,20],(z,B,ez,eB))
optparams=optrslt.x
f=np.polynomial.polynomial.Polynomial(optparams,domain=(min(z),max(z)))
linfitcoords=f.linspace(1000)
#Resiudals
if True:
residuals=f(z)-B
#Residual level line
if True:
rf=np.polynomial.polynomial.Polynomial([0],domain=(min(z),max(z)))
rfcoords=rf.linspace(1000)
#Plotting
if True:
fig=plt.figure(1)
frame=plt.gca()
axis1=fig.add_axes((0,0,1,1))
axis1.errorbar(z,B,xerr=ez,fmt='o', ms=1, lw=0.5, color='red', capsize=2, capthick=0.5, ecolor='black')
#axis1.plot(z,line(z,optparams[1],optparams[0]),lw=0.5,color='blue')
#axis1.plot(x,y,lw=0.5,color='blue')
axis1.plot(linfitcoords[0],linfitcoords[1],lw=0.5,color='blue')
plt.ylabel('Magnetic Field Strength ($mT$)')
axis1.xaxis.set_visible(False)
axis1.set_xlim(left=0,right=0.025)
axis1.set_ylim(bottom=0,top=0.5)
axis2=fig.add_axes((0,-0.4,1,0.4))
axis2.set_xlim(left=0,right=0.025)
axis2.plot(rfcoords[0],rfcoords[1],lw=0.5, color='blue')
axis2.scatter(z,residuals, s=1, color='red')
plt.ylabel("Residuals ($mT)$")
axis1.yaxis.set_label_coords(-0.1,0.5)
axis2.yaxis.set_label_coords(-0.1,0.5)
plt.xlabel('z $(cm^{-3})$')
plt.xlim([0,0.025])
axis3=fig.add_axes((1,-0.4,0.2,0.4))
axis3.hist(residuals,orientation='horizontal',color='gray')
axis3.yaxis.set_visible(False)
axis3.xaxis.tick_top()
axis3.xaxis.set_label_position('top')
plt.xlabel('Frequency')
plt.xticks([1,3,5,7])
#Note:restarting kernel (close window) fixes x axis label problem where you accidentally set as variable and then it breaks
#Save Figure
if True:
plt.savefig("Linearised_Biot_Savart_Graph.png",dpi=1000,bbox_inches='tight')
#Gradient and mu_0 analysis
if True:
m=optparams[1]
mu0=2*m/(I*a**2)
em=1#need to actually figure this out
emu0=mu0*np.sqrt((em/m)**2+(eI/I)**2+4*(ea/a)**2)
# print(mu0)
#print(emu0)
#Normality of residuals
if True:
print(stats.shapiro(residuals))
print(stats.normaltest(residuals))
#print(stats.kstest(residuals,stats.norm(np.mean(residuals),np.var(residuals)).cdf))
注意:所有的 if 语句只是为了方便 Spyder 中的代码折叠。
对于上下文,我正在对磁场的一些物理数据进行一些数据分析和可视化。我的主要问题是我试图绘制一条最适合数据的线,即使我得到的参数对于最适合的线是正确的,当我尝试使用 numpy 多项式生成坐标时绘制最佳拟合线,它生成错误的坐标。
问题似乎出在传递给 f
domain
参数上
从 docs 中,您正在为域传递值,但没有为 window 传递值,因此 f
根据地图 [min(z), max(z)] -> [-1, 1]
[=18 映射输入=]
这就是输出不同的原因。
要清楚地看到这一点,请看下面的示例
import numpy as np
# x**2 with no domain shift
f = np.polynomial.polynomial.Polynomial([0,0,1])
print(f)
>> poly([0. 0. 1.])
print(f(1))
>> 1.0
# x**2 with domain shift [-0.1, 0.1] -> [-1, 1], so an input of 1 maps to 10
g = np.polynomial.polynomial.Polynomial([0,0,1], domain=[-0.1, 0.1])
print(g)
>> poly([0. 0. 1.])
print(g(1))
>> 100.0
因此,任何不是 [-1,1]
的 domain
值都会给出一个多项式,由于输入