linregress and polyfit

Josef Pktd josef.pktd at gmail.com
Wed Sep 18 00:25:17 EDT 2013


On Tuesday, September 17, 2013 10:34:39 PM UTC-4, Krishnan wrote:
> I created an xy pair
> 
> 
> 
> y = slope*x + intercept
> 
> 
> 
> then I added some noise to "y" using
> 
> 
> 
> numpy.random.normal - call it z 
> 
> 
> 
> I could recover the slope, intercept from (x,y) using linregress
> 
> BUT cannot calculate the slope, intercept from (x, z)
> 
> 
> 
> What is puzzling is that for both pairs (x,y) and (x,z) the
> 
> polyfit (order 1) works just fine (gives me the slope, intercept)
> 
> ----------------------------------------------------------------------
> 
> import numpy as np
> 
> import scipy
> 
> # create a straight line, y= slope*x + intercept 
> 
> x = np.linspace(0.0,10.0,21)  
> 
> slope = 1.0  
> 
> intercept = 3.0   
> 
> y = []  
> 
> for i in range(len(x)):  
> 
>     y.append(slope*x[i] + intercept)  
> 
> # now create a data file with noise
> 
> z= []  
> 
> for i in range(len(x)):  
> 
>     z.append(y[i] + 0.1*np.random.normal(0.0,1.0,1))  

When z is converted to a numpy array then it has an extra dimension that linregress cannot handle, because np.random.normal(0.0,1.0, 1) returns an array and not a scalar.

much easier: use vectorized numpy instead of loop

z = y + 0.1*np.random.normal(0.0,1.0, len(y))

which is a one dimensional array and works with linregress

Josef

> 
> # now calculate the slope, intercept using linregress
> 
> from scipy import stats  
> 
> # No error here this is OK, works for x, y
> 
> cslope, cintercept, r_value, p_value, std_err = stats.linregress(x,y)
> 
> print cslope, cintercept
> 
> # I get an error here
> 
> #ValueError: array dimensions must agree except for d_0
> 
> nslope, nintercept, nr_value, np_value, nstd_err = stats.linregress(x,z)  
> 
> print nslope, nintercept
> 
> # BUT polyfit works fine, polynomial or order 1 with both data sets
> 
> ncoeffs = scipy.polyfit(x,z,1)  
> 
> print ncoeffs
> 
> coeffs = scipy.polyfit(x,y,1)  
> 
> print coeffs



More information about the Python-list mailing list