Python实现：时间序列趋势外推法应用-龚珀兹曲线拟合

news2024/12/27 2:09:47

龚珀兹曲线

$\hat{y}=ka^{b^{t}}$

下表数据为某跨国公司1989-2021年的年销售量数据，使用适合的模型预测该公司2022年的销售额，并得出理由。

部分数据如下表（具体数据从主页资源下载）：

年份	时序（t）	总额（yt）	时序应该从0开始
1989	1	138.40	0
1990	2	174.00	1
1991	3	190.55	2
1992	4	196.10	3
1993	5	230.50	4
1994	6	237.10	5
1995	7	274.00	6
1996	8	319.00	7
1997	9	348.45	8
1998	10	303.85	9

1、读取数据：

#读取数据
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

data = pd.read_excel('ch4综合分析.xlsx')
data.head()
#%%
# 创建散点图
#显示中文
from matplotlib import rcParams

# 设置中文字体
rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

plt.scatter(data.index, data['总额（yt）'], color='blue', label='原始数据')
plt.xlabel('时序（t）')
plt.ylabel('总额（yt）')
plt.legend()
plt.show()

通过散点图得知数据符合龚珀兹曲线特征，为此采用龚珀兹曲线拟合。

2、原理

（1）取对数：

$lgy=lgk+b^{^{t}}lga$

（2）分组法求解参数

$b^n=\frac{\sum II I\lg y-\sum I I \lg y}{\sum I I \lg y-\sum I \lg y}$

$lg a=\left(\sum II\lg y-\sum I \lg y\right) \cdot \frac{b-1}{\left(b^n-1\right)^2}$

$lg k=\frac{1}{n}\left(\sum I \lg y-\frac{b^n-1}{b-1} \cdot \lg a\right)$

或者

$lg k=\frac{1}{n}\left[\frac{\sum I \lg y \cdot \sum III\lg y-\left(\sum II\lg y\right)^2}{\sum I \lg y+\sum II \lg y-2 \sum III\lg y}\right]$

（3）Python求解

#选择合适的趋势外推法预测销售额
#拟合
# 对 '总额（yt）' 列取对数并写入表格
data['总额（yt）_log'] = np.log10(data['总额（yt）'])
# 显示前几行数据
print(data.head(33))
print('----------------------------------------------------')
#使用分组法估计参数
#分组
n = (2021-1989+1)/3
#读取前n个数据
data1 = data.iloc[:11]
data2 = data.iloc[11:22]
data3 = data.iloc[22:]
#分组求和
data1_sum = data1['总额（yt）_log'].sum()
data2_sum = data2['总额（yt）_log'].sum()
data3_sum = data3['总额（yt）_log'].sum()
#计算参数
b = ((data3_sum-data2_sum)/(data2_sum-data1_sum))**(1/11)
a = (data2_sum-data1_sum)*((b-1)/(b**11-1)**2)
k = (data1_sum-a*(b**11-1)/(b-1))*(1/11)

a1=10**a
k1=10**k
#输出参数
print('k1 = %f, a1 = %f, b = %f' % (k1,a1,b))

#输出拟合函数
print('拟合函数为：y = %f * %f ^ %f^t' % (k1,a1,b))

#计算第2022年的销售量
t1=33
d34 = k1 *a1 ** b**t1
print('2022年的销售量为：%d' % d34)
#计算se
# 计算预测值
data['预测值'] = k1 * a1 ** (b ** data.index)
# 计算残差
data['残差'] = data['总额（yt）'] - data['预测值']
# 计算残差的标准误差
se = np.sqrt(np.sum(data['残差'] ** 2) / (len(data) - 2))
# 输出标准误差
print('标准误差 (SE) = %f' % se)

（4）龚珀兹曲线方程

$y=110.036 \times 1.831^{1.047^{t}}$

3、曲线拟合

import matplotlib.pyplot as plt
# 计算拟合值
data['拟合值'] = k1 * a1 ** (b ** data.index)
# 创建散点图和拟合曲线
plt.scatter(data.index, data['总额（yt）'], color='blue', label='原始数据')
plt.plot(data.index, data['拟合值'], color='red', label='拟合曲线')
#添加预测点
plt.scatter(t1, d34, color='green', label='2022年预测点')
plt.xlabel('时序（t）')#时序从第0开始计时
plt.ylabel('总额（yt）')
plt.legend()
plt.show()