I am reading some R code about using the poly function but I cannot figure out how it works.
I know this:
ttt = 1:10
tttMod = poly(ttt, degree = 3, raw = FALSE)
tttMod
1 2 3
[1,] -0.49543369 0.52223297 -0.4534252
[2,] -0.38533732 0.17407766 0.1511417
[3,] -0.27524094 -0.08703883 0.3778543
[4,] -0.16514456 -0.26111648 0.3346710
[5,] -0.05504819 -0.34815531 0.1295501
[6,] 0.05504819 -0.34815531 -0.1295501
[7,] 0.16514456 -0.26111648 -0.3346710
[8,] 0.27524094 -0.08703883 -0.3778543
[9,] 0.38533732 0.17407766 -0.1511417
[10,] 0.49543369 0.52223297 0.4534252
attr(,"coefs")
attr(,"coefs")$alpha
[1] 5.5 5.5 5.5
attr(,"coefs")$norm2
[1] 1.0 10.0 82.5 528.0 3088.8
attr(,"degree")
[1] 1 2 3
attr(,"class")
[1] "poly" "matrix"
I cannot understand the following code
ttt2 = 1:5
ttt2Predicted = predict(tttMod, ttt2)
ttt2Predicted
1 2 3
[1,] -0.49543369 0.52223297 -0.4534252
[2,] -0.38533732 0.17407766 0.1511417
[3,] -0.27524094 -0.08703883 0.3778543
[4,] -0.16514456 -0.26111648 0.3346710
[5,] -0.05504819 -0.34815531 0.1295501
I don't know what (and how) coefficients from the tttMod the predict function used to generate the ttt2Predicted. How can I get the same result (i.e. ttt2Predicted) by not using the predict(tttMod, ttt2)? For example, doing an equivalent linear algebra calculation.
I tried using poly but seems it doesn't work in this way
ttt2Predicted_test = poly(ttt2, degree = 3, raw = FALSE)
1 2 3
[1,] -0.6324555 0.5345225 -3.162278e-01
[2,] -0.3162278 -0.2672612 6.324555e-01
[3,] 0.0000000 -0.5345225 -4.095972e-16
[4,] 0.3162278 -0.2672612 -6.324555e-01
[5,] 0.6324555 0.5345225 3.162278e-01
attr(,"coefs")
attr(,"coefs")$alpha
[1] 3 3 3
attr(,"coefs")$norm2
[1] 1.0 5.0 10.0 14.0 14.4
attr(,"degree")
[1] 1 2 3
attr(,"class")
[1] "poly" "matrix"
Update
Though the question has been marked duplicate, I provided the Python implementation of R' poly function and its prediction. Hope you will also find this helpful.
With help from the Stack Overflow community, I have implemented the equivalent Python code for the following R code
ttt = 1:10
tttMod = poly(ttt, degree = 3, raw = FALSE)
ttt2 = 1:5
ttt2Predicted = predict(tttMod, ttt2)
The equivalent Python snippet is translated from here and given below,
import numpy as np
import sys
def poly(x, degree):
'''x: ndarray
degree: scalar'''
## check feasibility
if len(set(x)) < degree:
print('Insufficient unique data points for specified degree!')
sys.exit()
centre = np.mean(x)
x = x - centre
Q, R = np.linalg.qr(np.polynomial.polynomial.polyvander(x, degree))
h, tau = np.linalg.qr(np.polynomial.polynomial.polyvander(x, degree), 'raw')
## rescaling of Q by h.T
X = (Q * np.diag(h.T))[:, 1:]
X2 = X ** 2
norm2 = X2.sum(0)
alpha = np.dot(X2.T, x) / norm2
beta = norm2 / (np.append([len(x)], norm2[:-1]))
scale = np.sqrt(norm2)
X = X * (1 / scale)
return X, centre, scale, alpha[:-1], beta[:-1]
def predict_poly(X, centre, scale, alpha, beta, x):
'''
X, centre, scale, alpha, beta: returns from poly function
x: new data to be predicted
'''
degree = X.shape[1]
x = x - centre
X = np.array(x.tolist() * degree).reshape((len(x), degree), order='F')
if degree > 1:
## second polynomial
X[:, 1] = (x - alpha[0]) * X[:, 0] - beta[0]
## further polynomials obtained from recursion
i = 2
while i < degree:
X[:, i] = (x - alpha[i-1]) * X[:, i-1] - beta[i-1] * X[:, i-2]
i += 1
return X * (1 / scale)
else:
return None
#### test 1
x = np.array([1.160, 1.143, 1.126, 1.109, 1.079, 1.053, 1.040, 1.027, 1.015, 1.004, 0.994, 0.985, 0.977])
X, centre, scale, alpha, beta = poly(x, 5)
x = np.array([1.141096, 1.025588, 1.045099, 1.081832, 1.143202])
predict_poly(X, centre, scale, alpha, beta, x)
#### test 2
x = np.arange(1, 11, 1)
X, centre, scale, alpha, beta = poly(x, 3)
x = np.arange(1, 6, 1)
predict_poly(X, centre, scale, alpha, beta, x)
#### result of test 2
array([[-0.49543369, 0.52223297, -0.45342519],
[-0.38533732, 0.17407766, 0.15114173],
[-0.27524094, -0.08703883, 0.37785433],
[-0.16514456, -0.26111648, 0.33467098],
[-0.05504819, -0.34815531, 0.12955006]])
