最小二乘法
使用最小二乘法拟合大气二氧化碳浓度数据
数据保存在monthly_co2.xls
文件中(只截取部分)
python需要安装的库
- xlrd
- numpy
- pandas
- matplotlib
绘制图像代码(绘制整体数据趋势图)
# -*- coding: utf-8 -*-
"""
@File : 绘制趋势图.py
@Time : 2024/11/27 23:52:44
@Version :
@Desc :
""""""
@Python version : 3.8.7
@matplotlib version : 3.5.1
@pandas version : 1.4.2
@numpy version : 1.22.3
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator #从ticker中导入MultipleLocator类df = pd.read_excel(io="./monthly_co2.xls", header=2) # 数据在第一个sheet中# 删除索引为0的行
df.drop(index=0, inplace=True) # 删除第三行数据 直接在原始dataframe上修改
df.replace(-99.99, np.nan, inplace=True) # 处理缺失值
df.dropna(inplace=True) # 删除缺失值的行# df['Year_Month'] = df['Yr'] + (df['Mn'] - 1) / 12 # 将年份和月份转换为小数形式
# [:, 0] 表示选择所有行
df['Year_Month'] = df.iloc[:, 0] + (df.iloc[:, 1] - 1) / 12 # 使用列index索引而不是名称索引x_vec = df['Year_Month'].values # x
y_vec = df.iloc[:, 2].values # y# matplotlib字体设置
plt.rcParams['font.family'] = "Times New Roman" # 设置全局字体# marker='o' 散点图绘制为圆形
# edgecolor='r' 设置圆圈的边缘颜色为红色
custom_color = '#f0a1a8'
plt.scatter(x_vec, y_vec, marker='o', edgecolors=custom_color, facecolors='none', label='Actual Data') # 原始数据
plt.xlabel("Year_Month", fontsize=16)
plt.ylabel("CO2 Concentration [ppm]", fontsize=16)
plt.title("Trend Chart of CO2 Concentration Over Time", fontsize=20) # 二氧化碳浓度随时间变化趋势图
plt.legend(loc='upper left')# 美化图表plt.rcParams['axes.facecolor'] = 'lightgray'
plt.rcParams['legend.fontsize'] = 10 # 设置图例字体大小# 刻度值设置
ax = plt.gca() # 获取当前坐标轴
ax.xaxis.set_major_locator(MultipleLocator(10)) # 设置x轴的主刻度间隔
ax.xaxis.set_minor_locator(MultipleLocator(5)) # 设置x轴的次刻度间隔ax.yaxis.set_major_locator(MultipleLocator(20)) # 设置y轴的主刻度间隔
ax.yaxis.set_minor_locator(MultipleLocator(10)) # 设置y轴的次刻度间隔plt.grid(True, color='gray', linestyle='--', linewidth=0.5) # 添加自定义样式网格线
plt.show()
结果(请自己修改线形、标题等内容)
二次函数拟合
# -*- coding: utf-8 -*-
"""
@File : 二次函数拟合结果.py
@Time : 2024/11/28 11:19:38
@Version :
@Desc :
""""""
@Python version : 3.8.7
@matplotlib version : 3.5.1
@pandas version : 1.4.2
@numpy version : 1.22.3
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator #从ticker中导入MultipleLocator类df = pd.read_excel(io="./monthly_co2.xls", header=2) # 数据在第一个sheet中# 删除索引为0的行
df.drop(index=0, inplace=True) # 删除第三行数据 直接在原始dataframe上修改
df.replace(-99.99, np.nan, inplace=True) # 处理缺失值
df.dropna(inplace=True) # 删除缺失值的行# df['Year_Month'] = df['Yr'] + (df['Mn'] - 1) / 12 # 将年份和月份转换为小数形式
# [:, 0] 表示选择所有行
df['Year_Month'] = df.iloc[:, 0] + (df.iloc[:, 1] - 1) / 12 # 使用列index索引而不是名称索引x_vec = df['Year_Month'].values # x
y_vec = df.iloc[:, 2].values # y# 这里添加代码M = df.shape[0] # 数据点个数
# 方程组的系数矩阵Phi = np.zeros((M, 3)) # 创建空矩阵
for i in range(M):Phi[i][0] = 1Phi[i][1] = x_vec[i]Phi[i][2] = x_vec[i]**2# 正规方程组A = Phi.T @ Phi
b = Phi.T @ y_vecx = np.linalg.solve(A, b)print(x)a_0 = x[0]
a_1 = x[1]
a_2 = x[2]# f_x = a_2x^2 + a_1x + a_0
varphi_x = lambda x: a_2 * x**2 + a_1 * x + a_0 # 拟合函数# 可视化图像y_vec_model = varphi_x(x_vec)# 均方根误差
RMSE = np.sqrt(np.sum((y_vec - y_vec_model)**2) / M)
print("均方根误差: ", RMSE)# matplotlib字体设置
plt.rcParams['font.family'] = "Times New Roman" # 设置全局字体plt.plot(x_vec, y_vec_model-279, label='Fitted Line', color='red')# marker='o' 散点图绘制为圆形
# edgecolor='r' 设置圆圈的边缘颜色为红色
custom_color = '#f0a1a8'
plt.scatter(x_vec, y_vec-279, marker='o', edgecolors=custom_color, facecolors='none', label='Actual Data') # 原始数据
plt.xlabel("Year_Month", fontsize=16)
plt.ylabel("CO2 Concentration [ppm]", fontsize=16)
plt.title("Trend Chart of CO2 Concentration Over Time", fontsize=20) # 二氧化碳浓度随时间变化趋势图
plt.legend(loc='upper left')# 美化图表plt.rcParams['axes.facecolor'] = 'lightgray'
plt.rcParams['legend.fontsize'] = 10 # 设置图例字体大小# 刻度值设置
ax = plt.gca() # 获取当前坐标轴
ax.xaxis.set_major_locator(MultipleLocator(10)) # 设置x轴的主刻度间隔
ax.xaxis.set_minor_locator(MultipleLocator(5)) # 设置x轴的次刻度间隔ax.yaxis.set_major_locator(MultipleLocator(20)) # 设置y轴的主刻度间隔
ax.yaxis.set_minor_locator(MultipleLocator(10)) # 设置y轴的次刻度间隔plt.grid(True, color='gray', linestyle='--', linewidth=0.5) # 添加自定义样式网格线# 添加RMSE文本信息
# (0.02, 0.85) 以轴的百分比来表示 分别对应x轴和y轴的为转移
plt.text(0.09, 0.95, f"RMSE: {RMSE:.2f}", transform=plt.gca().transAxes, fontsize=18)plt.show()
程序运行结果