我试图在Python中匹配以下R代码中的正交多项式：
X <- cbind(1, poly(x = x, degree = 9))但是在Python中实现。
为此，我实现了自己的方法来生成正交多项式：
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但似乎并不匹配。有人知道它使用的是哪种正交多项式吗？我查了文档但没有找到相关信息。
为了提供一些背景，我正在尝试将以下R代码转换为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 <- 10x <- seq(from = 0, to = 1, length = N)mu <- sin(2 * pi * x * 4)y <- muplot(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 <- 2for(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)))因为这似乎是唯一一个能够与使用多项式特征的梯度下降和线性回归一起工作的代码。
我觉得任何正交多项式（至少是标准正交）都应该能工作，并给出一个条件数为1的海森矩阵，但在Python中我似乎无法实现。相关问题：如何使用厄米多项式进行随机梯度下降（SGD）？
回答：
poly 使用QR分解，如这个回答中详细描述的那样。
我认为您真正想要的是如何使用Python复制R的poly的输出。
这里我编写了一个基于R实现的函数来实现这一点。我还添加了一些注释，以便您可以看到相应的R语句是什么样的：
import numpy as npdef 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];检查这个是否工作：
x = np.arange(10) + 1degree = 9poly(x, degree)返回矩阵的第一行是
[-0.49543369, 0.52223297, -0.45342519, 0.33658092, -0.21483446, 0.11677484, -0.05269379, 0.01869894, -0.00453516],与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