自学内容网 自学内容网

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

龚珀兹曲线

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

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

部分数据如下表(具体数据从主页资源下载):

年份时序(t)总额(yt)时序应该从0开始
19891138.400
19902174.001
19913190.552
19924196.103
19935230.504
19946237.105
19957274.006
19968319.007
19979348.458
199810303.859

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()


原文地址:https://blog.csdn.net/2301_76574743/article/details/142579999

免责声明:本站文章内容转载自网络资源,如本站内容侵犯了原著者的合法权益,可联系本站删除。更多内容请关注自学内容网(zxcms.com)!