其中a = 0.5, b = 3。我们选取高斯噪声标准差为5,通过如下代码模拟数据并绘制二维散点图:
import numpy as npimport matplotlib.pyplot as plt from scipy.optimize import minimizea =0.5# coefficientb =3# interceptnoiseVar =1# noise的标准差x = np.arange(1, 100, 0.1)# generate datay = a * x + b + noiseVar* np.random.randn(x.size)print(x.size)plt.plot(x, y, 'o', color='r')
x2 = np.vstack((x, np.ones(x.size)))# x2的第二排都是1# 因为numpy都是横向量,我们划成列向量x2 = x2.T # x2 is 99 x 2y_hat = y[:, np.newaxis].copy()# y is 99 x 1res = np.linalg.inv(x2.T @ x2)@ x2.T @ y_hat print('Coefficient is ',res[0][0])print('Intercept is ',res[1][0])
结果显示a = 0.5009766233665898, b = 2.954706934764938。可以看到,这个结果同样也是比较精确的。
方法三:利用算法包直接求解。
有些算法包已经提供了对于特定模型的专用优化函数,我们可以直接调用这些函数来得到最优解。
具体操作步骤如下:
x2 = np.vstack((x, np.ones(x.size)))# x2的第二排都是1# 利用lstsq函数求解res = np.linalg.lstsq(x2.T, y, rcond=None)[0]print('Coefficient is ',res[0])print('Intercept is ',res[1])
结果显示a = 0.5009766233665904, b = 2.9547069347648884。同样得到了较为不错的结果。