We can test this behavior doing the least-squares fit of a y = a*x**0 + b*x**1 + c*x**2 equation that we know the answer should be a=0, b=0, c=1:
np.linalg.lstsq([ [1, 1, 1], [1, 2, 4], [1, 3, 9] ], [1, 4, 9]) #(array([-3.43396424e-15, 3.88578059e-15, 1.00000000e+00]), # array([], dtype = float64), # 3, # array([10.64956309, 1.2507034, 0.15015641]))
Fit a polynomial p(x) = p[0] * x**deg + ... + p[deg] of degree deg to points (x, y). Returns a vector of coefficients p that minimises the squared error in the order deg, deg-1, … 0.,residuals – sum of squared residuals of the least squares fit,The rank of the coefficient matrix in the least-squares fit is deficient. The warning is only raised if full == False.,Switch determining nature of return value. When it is False (the default) just the coefficients are returned, when True diagnostic information from the singular value decomposition is also returned.
>>>
import warnings
>>>
warnings.simplefilter('ignore', np.RankWarning)x[0] ** n * p[0] + ...+x[0] * p[n - 1] + p[n] = y[0] x[1] ** n * p[0] + ...+x[1] * p[n - 1] + p[n] = y[1] ... x[k] ** n * p[0] + ...+x[k] * p[n - 1] + p[n] = y[k]
>>> import warnings >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0]) >>> z = np.polyfit(x, y, 3) >>> z array([0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary
>>> p = np.poly1d(z) >>> p(0.5) 0.6143849206349179 # may vary >>> p(3.5) - 0.34732142857143039 # may vary >>> p(10) 22.579365079365115 # may vary
>>> with warnings.catch_warnings():
...warnings.simplefilter('ignore', np.RankWarning)
...p30 = np.poly1d(np.polyfit(x, y, 30))
...
>>>
p30(4) -
0.80000000000000204 # may vary >>>
p30(5) -
0.99999999999999445 # may vary >>>
p30(4.5) -
0.10547061179440398 # may vary>>>
import matplotlib.pyplot as plt >>>
xp = np.linspace(-2, 6, 100) >>>
_ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--') >>>
plt.ylim(-2, 2)
(-2, 2) >>>
plt.show()2020-02-20
For linear functions, we have this formula:
y = a * x + b
In this equation, usually, a and b are given. E.g:
a = 2 b = 5
So:
y = 2 * x + 5
Nice, we got a line that we can describe with a mathematical equation – this time, with a linear function. The general formula was:
y = a * x + b
And in this specific case, the a and b values of this line are:
a = 2.01 b = -3.9
We can test this behavior doing the least-squares fit of a y = a*x**0 + b*x**1 + c*x**2 equation that we know the answer should be a=0, b=0, c=1:,How to skirt this problem? Well, this version of polyfit accepts weights. You could add another datapoint but weight it at epsilon. This is equivalent to reducing the 2.0 in this formula to a 1.0.,It is possible to do so without iterate along the first axis. However, your second axis is rather short (being just 3), you can really fit no more than 2 coefficients. ,Without loops, this is computed by multiplying the covariance matrix C on both sides by the vector x_powers made of the relevant powers of x: for example, [x**5, x**4, ..., x**0]. The setup of polyfit is included for completeness.
x = [-449., -454., -459., -464., -469.] y = [0.9677024, 0.97341953, 0.97724978, 0.98215678, 0.9876293] x_extra = x + x[-1: ] y_extra = y + y[-1: ] weights = [1.0, 1.0, 1.0, 1.0, 1.0, sys.float_info.epsilon] fit, cov = np.polyfit(x, y, 2, cov = True) fit_extra, cov_extra = np.polyfit(x_extra, y_extra, 2, w = weights, cov = True) print fit == fit_extra print cov_extra
[p, ~, mu] = polyfit(T.year, T.pop, 5);
In[67]:
import numpy as np
import scipy.optimize as so
In[68]:
def MD_ployError(p, x, y):
''
'if x has the shape of (n,m), y must be (n,m), p must be (n*p, ), where p is degree'
''
#d is no.of degree
p_rshp = p.reshape((x.shape[0], -1))
f = y * 1.
for i in range(p_rshp.shape[1]):
f -= p_rshp[: , i][: , np.newaxis] * (x ** i)
return (f ** 2).sum()
In[69]:
X = np.random.random((100, 6))
Y = 4 + 2 * X + 3 * X * X
P = (np.zeros((100, 3)) + [1, 1, 1]).ravel()
In[70]:
MD_ployError(P, X, Y)
Out[70]:
11012.2067606684
In[71]:
R = so.fmin_slsqp(MD_ployError, P, args = (X, Y))
Iteration limit exceeded(Exit mode 9) #you can increase iteration limit, but the result is already good enough.
Current
function value: 0.00243784856039
Iterations: 101
Function evaluations: 30590
Gradient evaluations: 101
In[72]:
R.reshape((100, -1))
Out[72]:
array([
[3.94488512, 2.25402422, 2.74773571],
[4.00474864, 1.97966551, 3.02010015],
[3.99919559, 2.0032741, 2.99753804],
..............................................)np.linalg.lstsq([ [1, 1, 1], [1, 2, 4], [1, 3, 9] ], [1, 4, 9]) #(array([-3.43396424e-15, 3.88578059e-15, 1.00000000e+00]), # array([], dtype = float64), # 3, # array([10.64956309, 1.2507034, 0.15015641]))
xdata = np.arange(10) ydata = np.abs(xdata - 2) # some data to fit degree = 5 p, C = np.polyfit(xdata, ydata, deg = degree, cov = True) x = 5.4 # some value of x to estimate error at x_powers = x ** np.arange(degree, -1, -1) x_error = np.sqrt(x_powers.dot(C).dot(x_powers))
x1 = x - np.mean(x) z1, cov1 = np.polyfit(np.asarray(x1), np.asarray(y), 1, cov = True) std1 = np.sqrt(np.diag(cov1)) print z1 # prints: array([1.56607841e+03, 6.31224162e+06]) print cov1 # prints: array([ [4.56066546e+00, -2.90980285e-07], #[-2.90980285e-07, 3.36480951e+00] ]) print std1 # prints: array([2.13557146, 1.83434171])
However the numpy.polyfit documentation defines the weight as “weights to apply to the y-coordinates”. This definition is not correct. The weights apply to (=multiply) the fit residuals, not only to the y-coordinates.,More importantly, looking at the math in the Numpy (v1.9.1) code, the resulting definition of squared residuals adopted by polyfit is the following, with the optional input weights vector w inside the parenthesis, contrary to standard practice,The confusion in the documentation likely arises from the fact that the Numpy code solves the linear problem below in the last-squares sense, where the w vector does multiply the y-coordinate,I did try to note that the “sigma” way of thinking about the weights was specific to “When these correspond to weighting by Gaussian uncertainties,…” so as not to imply that this was the only use case.
- Replace the doc for
w:
Weights to apply to the unsquared residual of each point.When these correspond to weighting
by Gaussian uncertainties, use w = 1 / sigma.Note that w here is the sqrt of the more typical definition of the weights in weighted least squares, which generally uses W = w ** 2 = 1 / sigma ** 2(or more generally 1 /
var) to weight the squared residual of each point.- Adjust slightly the final sentence of the doc for
cov:
This scaling is omitted
if cov = 'unscaled', as is relevant
for the
case that the weights are w ** 2 = 1 / sigma ** 2,
with sigma known to be a reliable estimate of the uncertainty.