在书中Matlab源码中, 对最小二乘法的求解是使用t=p\y
语句表示的,在Python中我们可以调用scipy库中的nnls方法来实现同样的功能。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import nnls
n = 50
N = 1000
x = np.linspace(-3, 3, n).T
X = np.linspace(-3, 3, N).T
pix = np.pi * x
y = (np.sin(pix) / pix + 0.1 * x).reshape(n, -1) + np.dot(0.05, np.random.randn(n, 1))
p = np.empty([n, 31])
p[:, 1] = np.ones((n,))
P = np.empty([N, 31])
P[:, 1] = np.ones((N,))
for i in range(1, 16):
p[:, 2 * i - 1] = np.sin(i / 2 * x)
p[:, 2 * i] = np.cos(i / 2 * x)
P[:, 2 * i - 1] = np.sin(i / 2 * X)
P[:, 2 * i] = np.cos(i / 2 * X)
t = nnls(p, y.flatten())[0]
F = P.dot(t)
plt.plot(x, y, 'bo')
plt.plot(X, F, 'g-')
plt.axis([-3.05, 3.05, -0.5, 1.2])
plt.show()