python - evaluar - ¿Qué tipo de polinomios ortogonales usa R?

evaluar polinomios en python

Question

Estaba tratando de hacer coincidir los polinomios ortogonales en el siguiente código en R:

X <- cbind(1, poly(x = x, degree = 9))

pero en Python.

Para hacer esto, implementé mi propio método para dar polinomios ortogonales:

def get_hermite_poly(x,degree):
    #scipy.special.hermite()
    N, = x.shape
    ##
    X = np.zeros( (N,degree+1) )
    for n in range(N):
        for deg in range(degree+1):
            X[n,deg] = hermite( n=deg, z=float(x[deg]) )
    return X

aunque no parece coincidir. ¿Alguien sabe el tipo de polinomio ortogonal que usa? Intenté buscar en la documentación pero no dije nada.

Para dar un contexto, intento implementar el siguiente código R en python ( https://stats.stackexchange.com/questions/313265/issue-with-convergence-with-sgd-with-function-approximation-using-polynomial-lin/315185#comment602020_315185 ):

set.seed(1234)

N <- 10
x <- seq(from = 0, to = 1, length = N)
mu <- sin(2 * pi * x * 4)
y <- mu
plot(x,y)

X <- cbind(1, poly(x = x, degree = 9))
# X <- sapply(0:9, function(i) x^i)
w <- rnorm(10)

learning_rate <- function(t) .1 / t^(.6)

n_samp <- 2
for(t in 1:100000) {
  mu_hat <- X %*% w
  idx <- sample(1:N, n_samp)
  X_batch <- X[idx,]
  y_batch <- y[idx]
  score_vec <- t(X_batch) %*% (y_batch - X_batch %*% w)

  change <- score_vec * learning_rate(t)
  w <- w + change
}

plot(mu_hat, ylim = c(-1, 1))
lines(mu)
fit_exact <- predict(lm(y ~ X - 1))
lines(fit_exact, col = ''red'')
abs(w - coef(lm(y ~ X - 1)))

porque parece ser el único que trabaja con pendiente de gradiente con regresión lineal con características polinomiales.

Siento que cualquier polinomio ortogonal (o al menos ortonormal) debería funcionar y dar una arpillera con la condición número 1, pero parece que no puedo hacer que funcione en python. Pregunta relacionada: ¿Cómo se usan polinomios de Hermite con Gradiente de Gradiente Estocástico (SGD)?

Answer 1

poly usa la factorización QR, como se describe en detalle en esta respuesta .

Creo que lo que realmente parece estar buscando es cómo replicar el resultado de R''s poly usando python.

Aquí he escrito una función para hacer eso en función de la implementación de R. También agregué algunos comentarios para que pueda ver a qué se parecen las declaraciones equivalentes en R:

import numpy as np

def poly(x, degree):
    xbar = np.mean(x)
    x = x - xbar

    # R: outer(x, 0L:degree, "^")
    X = x[:, None] ** np.arange(0, degree+1)

    #R: qr(X)$qr
    q, r = np.linalg.qr(X)

    #R: r * (row(r) == col(r))
    z = np.diag((np.diagonal(r)))  

    # R: Z = qr.qy(QR, z)
    Zq, Zr = np.linalg.qr(q)
    Z = np.matmul(Zq, z)

    # R: colSums(Z^2)
    norm1 = (Z**2).sum(0)

    #R: (colSums(x * Z^2)/norm2 + xbar)[1L:degree]
    alpha = ((x[:, None] * (Z**2)).sum(0) / norm1 +xbar)[0:degree]

    # R: c(1, norm2)
    norm2 = np.append(1, norm1)

    # R: Z/rep(sqrt(norm1), each = length(x))
    Z = Z / np.reshape(np.repeat(norm1**(1/2.0), repeats = x.size), (-1, x.size), order=''F'')

    #R: Z[, -1]
    Z = np.delete(Z, 0, axis=1)
    return [Z, alpha, norm2];

Comprobando que esto funciona:

x = np.arange(10) + 1
degree = 9
poly(x, degree)

La primera fila de la matriz devuelta es

[-0.49543369,  0.52223297, -0.45342519,  0.33658092, -0.21483446,
  0.11677484, -0.05269379,  0.01869894, -0.00453516],

en comparación con la misma operación en R

poly(1:10, 9)
# [1] -0.495433694  0.522232968 -0.453425193  0.336580916 -0.214834462
# [6]  0.116774842 -0.052693786  0.018698940 -0.004535159