使用 Python 和 Statsmodels 进行 Cobb-Douglas 函数回归分析：模型拟合、R 平方计算和可视化

from cProfile import label
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from mpl_toolkits import mplot3d
# 以上导入绘图 matplotlib，数据处理 numpy，pandas，拟合 statsmodels 各个模块

def calculate_r_squared(actual, predicted):  # 计算 R**2
    ss_total = np.sum((actual - np.mean(actual))**2)
    ss_residual = np.sum((actual - predicted)**2)
    r_squared = 1 - (ss_residual / ss_total)
    return r_squared

data = pd.read_csv('cobb_douglas.csv')  # 下载数据

# model = smf.ols('lny ~ lnk + lnl', data=data)  # 以 lny 为响应变量，lnk 与 lnl 为特征变量进行 (1) 的回归
# results = model.fit()
# print(results.params)  # 输出回归结果
# print(results.params[0] + data.lnk * results.params[1] + data.lnl * results.params[2])  # 输出拟合值
# 以下为 (1) 绘图
# xx = np.linspace(data.lnk.min(), data.lnk.max(), 100)
# yy = np.linspace(data.lnl.min(), data.lnl.max(), 100)
# XX, YY = np.meshgrid(xx, yy)
# ZZ = results.params[0] + XX * results.params[1] + YY * results.params[2]
# fig = plt.figure()
# ax = plt.axes(projection='3d')
# ax.plot_surface(XX, YY, ZZ, rstride=8, cstride=8, alpha=0.4, cmap='pink_r')
# ax.scatter(data.lnk, data.lnl, data.lny, c='slateblue')
# ax.set_xlabel('lnk')
# ax.set_ylabel('lnl')
# ax.set_zlabel('lny')
# plt.show()
# 以上为 (1) 绘图
# 以下为 (2) 绘图
# ZZ = results.params[0] + data.lnk * results.params[1] + data.lnl * results.params[2]
# fig = plt.figure()
# plt.xlabel('Year')
# plt.ylabel('Output logarithm')
# plt.plot(data.year, ZZ, label='Fitted value')
# plt.plot(data.year, data.lny, label='Actual value')
# plt.legend()
# plt.show()
# 以上为 (2) 绘图
# (3)
# model = smf.ols('lny ~ lnk + lnl + lnk*lnl', data=data)  # 以 lny 为响应变量，lnk 与 lnl 为特征变量进行 (2) 的回归
# results = model.fit()
# print(results.params)  # 输出回归结果
# 以下为 (3) 绘图, 方法与 (1) 形同
# xx = np.linspace(data.lnk.min(), data.lnk.max(), 100)
# yy = np.linspace(data.lnl.min(), data.lnl.max(), 100)
# XX, YY = np.meshgrid(xx, yy)
# ZZ = results.params[0] + XX * results.params[1] + YY * results.params[2] + XX * YY * results.params[3]
# fig = plt.figure()
# ax = plt.axes(projection='3d')
# ax.plot_surface(XX, YY, ZZ, rstride=8, cstride=8, alpha=0.4, cmap='pink_r')
# ax.scatter(data.lnk, data.lnl, data.lny, c='slateblue')
# ax.set_xlabel('lnk')
# ax.set_ylabel('lnl')
# ax.set_zlabel('lny')
# plt.show()
# 以上为 (3) 绘图
# (4)
model = smf.ols('lny ~ I(lnk**2) + lnk + lnk*lnl + lnl + I(lnl**2) ', data=data)  # 以 lny 为响应变量，lnk 与 lnl 为特征变量进行 (2) 的回归
results = model.fit()
print(results.params)  # 输出回归结果
# 以下为 (4) 绘图, 方法与 (3) 形同
xx = np.linspace(data.lnk.min(), data.lnk.max(), 100)
yy = np.linspace(data.lnl.min(), data.lnl.max(), 100)
XX, YY = np.meshgrid(xx, yy)
# ZZ = results.params[0] + XX * XX * results.params[1] + XX * results.params[2] + YY * results.params[3] + XX * YY * results.params[4] + YY * YY * results.params[5]
# fig = plt.figure()
# ax = plt.axes(projection='3d')
# ax.plot_surface(XX, YY, ZZ, rstride=8, cstride=8, alpha=0.4, cmap='pink_r')
# ax.scatter(data.lnk, data.lnl, data.lny, c='slateblue')
# ax.set_xlabel('lnk')
# ax.set_ylabel('lnl')
# ax.set_zlabel('lny')
# plt.show()
ZZ = results.params[0] + data.lnk * data.lnk * results.params[1] + data.lnk * results.params[2] + data.lnl * results.params[3] + data.lnk * data.lnl * results.params[4] + data.lnl * data.lnl * results.params[5]
R2 = calculate_r_squared(data.lny, ZZ)
fig = plt.figure()
plt.xlabel('Year')
plt.ylabel('Output logarithm')
plt.plot(data.year, ZZ, label='Fitted value')
plt.plot(data.year, data.lny, label='Actual value')
plt.legend()
# plt.annotate(r'$R^2 = $' + str(R2), (1, 1), fontsize=12, color='red')
plt.text(4, 5, r'$R^2 = {:.2f}$'.format(R2), fontsize=12, color='red')
plt.show()
# 以上为 (4) 绘图
# plt.text(4, 5, r'$R^2 = $' + str(R2), fontsize=12, color='red') 注释在生成的图中显示 4 位小数，为什么，如何解决？内容：在使用 plt.text 函数时，可以通过使用字符串的格式化方法，将 R2 的值格式化为两位小数。修改代码如下：

plt.text(4, 5, r'$R^2 = {:.2f}$'.format(R2), fontsize=12, color='red')

这样就可以将 R2 的值保留两位小数并显示在图中。