Оценка numpy полинома возвращает неправильные значения - PullRequest
Новая
       7

Оценка numpy полинома возвращает неправильные значения

2 голосов
/ 15 марта 2020

У меня есть линейный numpy полином с коэффициентами c = 2.44717404e-03 и m = 1.88697661e + 01, но когда я пытаюсь его оценить, он дает неправильный вывод. 

f=np.polynomial.polynomial.Polynomial([2.44717404e-03, 1.88697661e+01])

например, 

f(0.015)

возвращает 

3.9646796502044523

вместо правильного результата 

0.28549366554

Console output demonstrating the problem

# -*- 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.

Для контекста я делаю некоторый анализ данных и визуализацию некоторых физических данных в полях magneti c. Моя главная проблема заключается в том, что я пытаюсь построить линию наилучшего соответствия для данных, и хотя параметры, которые я получаю, являются правильными для линии наилучшего соответствия, когда я пытаюсь использовать полином numpy для генерации координаты для построения линии наилучшего соответствия, он генерирует неправильные координаты.

1 Ответ

1 голос
/ 15 марта 2020

Похоже, проблема в том, что параметр domain передается в f

Из документов вы передаете значения для домена, но нет значений для окна, поэтому f отображает входные данные в соответствии с картой [min(z), max(z)] -> [-1, 1]

Именно поэтому выходные данные отличаются.

Чтобы увидеть это ясно, посмотрите на следующий пример 

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

Таким образом, любые значения domain, которые не являются [-1,1], дадут многочлен, который возвращает «неправильный» вывод из-за масштабирования ввода

...