|
|
|
|
|
|
|
|
|
|
|
import numpy as np
|
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
from X1 import random_points, compute_curveData
|
|
|
|
|
|
from X2 import draw_axis
|
|
|
|
|
|
|
|
|
|
|
|
# # (1) 定义x和y坐标数据
|
|
|
|
|
|
# x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
|
|
|
|
|
|
# y = np.array([1, 2, 3, 5, 8, 13, 21, 34, 55, 89])
|
|
|
|
|
|
#
|
|
|
|
|
|
# # (2) 定义多项式阶数m，这里以3阶多项式为例
|
|
|
|
|
|
# m = 2
|
|
|
|
|
|
#
|
|
|
|
|
|
# # 使用x的幂次来构造A矩阵
|
|
|
|
|
|
# A = np.vander(x, m + 1)
|
|
|
|
|
|
#
|
|
|
|
|
|
# # (3) 计算拟合系数theta，使用最小二乘法的正规方程
|
|
|
|
|
|
# theta = np.linalg.inv(A.T @ A) @ A.T @ y
|
|
|
|
|
|
# print(theta)
|
|
|
|
|
|
# # (4) 使用theta和x的幂次计算拟合曲线
|
|
|
|
|
|
# x_fit = np.linspace(min(x), max(x), 100)
|
|
|
|
|
|
# A_fit = np.vander(x_fit, m + 1)
|
|
|
|
|
|
# y_fit = A_fit @ theta
|
|
|
|
|
|
#
|
|
|
|
|
|
# # 创建图形对象并绘制拟合曲线和样本点
|
|
|
|
|
|
# fig, ax = plt.subplots()
|
|
|
|
|
|
# ax.plot(x_fit, y_fit, label=f'Polynomial Fit of degree {m}')
|
|
|
|
|
|
# ax.scatter(x, y, color='red', label='Sample Points')
|
|
|
|
|
|
# ax.set_xlabel('X')
|
|
|
|
|
|
# ax.set_ylabel('Y')
|
|
|
|
|
|
# ax.legend()
|
|
|
|
|
|
# plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# 编程4.5-----------------------------------------
|
|
|
|
|
|
|
|
|
|
|
|
# m次多项式的最小二乘法
|
|
|
|
|
|
def least_square_method(m, sampleData):
|
|
|
|
|
|
x = sampleData[:,0]
|
|
|
|
|
|
y = sampleData[:,1]
|
|
|
|
|
|
|
|
|
|
|
|
# (2) 多项式阶数m，这里以3阶多项式为例
|
|
|
|
|
|
# 使用x的幂次来构造A矩阵
|
|
|
|
|
|
A = np.vander(x, m + 1)
|
|
|
|
|
|
|
|
|
|
|
|
# (3) 计算拟合系数theta，使用最小二乘法的正规方程
|
|
|
|
|
|
theta = np.linalg.inv(A.T @ A) @ A.T @ y
|
|
|
|
|
|
|
|
|
|
|
|
# (4) 使用theta和x的幂次计算拟合曲线
|
|
|
|
|
|
x_fit = np.linspace(min(x), max(x), 100)
|
|
|
|
|
|
A_fit = np.vander(x_fit, m + 1)
|
|
|
|
|
|
y_fit = A_fit @ theta
|
|
|
|
|
|
|
|
|
|
|
|
plt.plot(x_fit, y_fit)
|
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
# 计算残差平方和
|
|
|
|
|
|
residuals = y - A @ theta
|
|
|
|
|
|
residuals_squared_sum = residuals.T @ residuals
|
|
|
|
|
|
# 计算均方误差
|
|
|
|
|
|
degrees_of_freedom = len(y) - (m + 1)
|
|
|
|
|
|
mean_squared_error = residuals_squared_sum / degrees_of_freedom
|
|
|
|
|
|
# 计算协方差矩阵
|
|
|
|
|
|
covariance_matrix = np.linalg.inv(A.T @ A) * mean_squared_error
|
|
|
|
|
|
|
|
|
|
|
|
return theta, covariance_matrix
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def curve_draw_line(curveData,sampleData):
|
|
|
|
|
|
x_fit = curveData[:,0]
|
|
|
|
|
|
y_fit = curveData[:,1]
|
|
|
|
|
|
x = sampleData[:,0]
|
|
|
|
|
|
y = sampleData[:,1]
|
|
|
|
|
|
# 创建图形对象并绘制拟合曲线和样本点
|
|
|
|
|
|
ax = draw_axis(0, 1000,step=250)
|
|
|
|
|
|
ax.plot(x_fit, y_fit, label=f'FitCurve')
|
|
|
|
|
|
ax.scatter(x, y, color='red', label='SampleData')
|
|
|
|
|
|
# ax.set_xlabel('X')
|
|
|
|
|
|
# ax.set_ylabel('Y')
|
|
|
|
|
|
plt.legend()
|
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
|
|
|
|
LimitNum = 1000
|
|
|
|
|
|
sampleData = random_points(20, 0, 1000)
|
|
|
|
|
|
# 这里为了便于观察使用设计好的数据
|
|
|
|
|
|
# x = np.array([0, 150, 250, 350, 450, 550, 650, 700, 750, 850])
|
|
|
|
|
|
# y = np.array([6, 30, 90, 160, 220, 340, 490, 620, 730, 1000])
|
|
|
|
|
|
# x = np.array([10, 100, 200, 300, 400, 500, 600, 700, 800, 900])
|
|
|
|
|
|
# y = np.array([10, 20, 30, 50, 80, 130, 210, 340, 550, 890])
|
|
|
|
|
|
# x = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90])
|
|
|
|
|
|
# y = np.array([1, 2, 1, 5, 8, 13, 21, 34, 55, 89])
|
|
|
|
|
|
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
|
|
|
|
|
|
y = np.array([1, 2, 1, 5, 8, 13, 21, 34, 55, 89])
|
|
|
|
|
|
sampleData = np.array(list(zip(x, y))) # 将两个一维数组拼接成二维数组
|
|
|
|
|
|
m = 3
|
|
|
|
|
|
theta, covariance_matrix = least_square_method(m, sampleData)
|
|
|
|
|
|
# print(theta)
|
|
|
|
|
|
curveData = compute_curveData(LimitNum, 1, theta, m, )
|
|
|
|
|
|
# print("curveData",curveData)
|
|
|
|
|
|
curve_draw_line(curveData,sampleData)
|
|
|
|
|
|
|
|
|
|
|
|
# 编程4.5 END-----------------------------------------
|
|
|
|
|
|
|
