Evaluating numpy polynomial returns the wrong values

Question

I have a linear numpy polynomial with coefficients c=2.44717404e-03 and m=1.88697661e+01, but when I try to evaluate it, it gives the wrong output.

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

e.g.

f(0.015)

returns

3.9646796502044523

instead of the correct result

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))

Note: all the if statements are just there to facilitate code folding in Spyder.

For context, I am doing some data analysis and visualisation of some physical data on magnetic fields. My main problem is that I'm trying to plot a line of best fit for the data, and even though the parameters I'm getting are correct for the line of best fit, when I try to use a numpy polynomial to generate coordinates to plot the best fit line, it generates the wrong coordinates.

Perhaps the order of coefficients the function expects is turned around? — phipsgabler, Mar 15 '20 at 09:52
Trying the code you put gives me the correct answer 0.28549..., are you sure you are not overwriting the function `f`? — KyleL, Mar 15 '20 at 09:53
@phipsgabler I've tried that and it doesn't help unfortunately. — Rational Function, Mar 15 '20 at 10:01
@KyleL I'm pretty certain I'm not overwriting it, I'll add a screenshot to demonstrate some of the weird behaviour. — Rational Function, Mar 15 '20 at 10:02
try printing out the values of `optparam` and `min(z)`, `max(z)` before you define `f`. It may be related to the `domain` parameter which is also being added — KyleL, Mar 15 '20 at 11:21
@KyleL Just tried that, and min(z) and max(z) are returning the correct values. — Rational Function, Mar 15 '20 at 11:43

KyleL · Accepted Answer · 2020-03-15T11:51:06.247

3

It looks like the issue is the domain parameter being passed to f

From the docs you are passing in values for the domain but no values for the window so f maps inputs according to the map [min(z), max(z)] -> [-1, 1]

This is why the output is different.

To see this clearly look at the following example

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

So any domain values which aren't [-1,1] will give a polynomial which returns the 'wrong' output due to scaling of the input

edited Mar 15 '20 at 11:51

answered Mar 15 '20 at 11:45

KyleL

855
6
24

Thank you very much for your reply! I vaguely understand the problem now, but I'm still not sure how to fix it. Do you know what I should set my window to in order to get the desired behaviour? – Rational Function Mar 15 '20 at 11:52
That depends what you want to do, if you set the window to the same as the domain it will do no scaling of the input, but this would change the current behaviour of the code and would be equivalent to not passing in any domain or window parameters – KyleL Mar 15 '20 at 12:00
By passing in `domain=[min(z), max(z)]` you are scaling the values in `z` to be -1 and 1 when applying the polynomial `f` – KyleL Mar 15 '20 at 12:02
1

using window=[min(z),max(z)] worked for me. I basically needed to set a domain in order to use the linspace method of the polynomial to generate my points, but I didn't actually want to do any rescaling, so setting the domain and window equal to each other so that their effects cancelled out resolved the issue. Thank you very much for your help!! – Rational Function Mar 15 '20 at 12:08

Evaluating numpy polynomial returns the wrong values

1 Answers1