在書(shū)中Matlab源碼中, 對(duì)最小二乘法的求解是使用t=p\y
語(yǔ)句表示的青团,在Python中我們可以調(diào)用scipy庫(kù)中的nnls方法來(lái)實(shí)現(xiàn)同樣的功能。
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()