在python上实现一个永不终止的矩阵公式用于曲线拟合

Question

我正试图编写一个程序来求解一般回归公式：所以我试图实现这个矩阵方程，有没有办法做到这一点，比如让用户决定它可以有多大，而不需要我制造越来越多的if条件(所以只有一段代码可以折叠成用户想要的矩阵(？代码：

#Solving the general matrix for the coefficients
if 3 == n:
a = np.array([[np.sum(np.multiply(FL[1],FL[1])),np.sum(np.multiply(FL[1],FL[2]))],
[np.sum(np.multiply(FL[1],FL[2])),np.sum(np.multiply(FL[2],FL[2]))]])
b = np.array([np.sum(np.multiply(FL[0],FL[1])),np.sum(np.multiply(FL[0],FL[2]))])
x = np.linalg.solve(a, b)
if 4 == n:
a = np.array([[np.sum(np.multiply(FL[1],FL[1])),np.sum(np.multiply(FL[1],FL[2])),np.sum(np.multiply(FL[1],FL[3]))],
[np.sum(np.multiply(FL[1],FL[2])),np.sum(np.multiply(FL[2],FL[2])),np.sum(np.multiply(FL[2],FL[3]))],
[np.sum(np.multiply(FL[1],FL[3])),np.sum(np.multiply(FL[2],FL[3])),np.sum(np.multiply(FL[3],FL[3]))]])
b = np.array([np.sum(np.multiply(FL[0],FL[1])),np.sum(np.multiply(FL[0],FL[2])),np.sum(np.multiply(FL[0],FL[3]))])
x = np.linalg.solve(a, b)

1在该代码中，Phi_0对应于FL[i=1]，FL[0]对应于y。

Answer 1

您可以使算法与多项式的阶数无关。最简单的方法是使用for循环，尽管这些循环会很慢(因为它们不利用NumPy的矢量化(
这里有一个可重复的随机数据示例：

import numpy as np
# Order of polynomial
n = 5
# Random seed for reproducibility
np.random.seed(1)
# Input arrays
phi = np.random.random((100,n))
y = np.random.random(100)
# Output arrays
a = np.zeros((n,n))
b = np.zeros(n)
for i in range(n):
b[i] = np.sum(y * phi[:,i])
for j in range(i,n):
# Exploit that matrix is diagonal
a[i,j] = a[j,i] = np.sum(phi[:,i] * phi[:,j])
# Coefficients array
x = np.linalg.solve(a,b)