【数学建模】 时间序列与投资模型

文章目录

【数学建模】 时间序列与投资模型

1. 时间序列的基本概念

时间序列(Time Series)是指按照时间顺序排列的一系列数据点。这些数据点通常是以等时间间隔记录的。时间序列分析是统计学和数据科学中的一个重要领域,用于分析、建模和预测随时间变化的数据。

1.1 时间序列的典型应用

时间序列分析在许多领域都有广泛的应用,包括但不限于以下几个方面:

  1. 经济学:如股票价格、利率、失业率等的分析和预测。
  2. 金融学:如金融市场的风险管理、投资组合管理等。
  3. 气象学:如天气预报、气候变化研究等。
  4. 工程学:如信号处理、质量控制等。
  5. 生物学:如心电图、脑电图等生物信号分析。
  6. 商业:如销售额预测、需求预测等。

1.2 时间序列的描述与分解

时间序列可以通过分解模型来描述,分解模型将时间序列分解为不同的成分,以便更好地理解和分析数据。这些成分通常包括:

  1. 趋势成分(Trend Component, T [ t ] T[t] T[t]:表示时间序列中长期的增长或下降趋势。
  2. 季节成分(Seasonal Component, S [ t ] S[t] S[t]:表示时间序列中固定周期内的波动,如一年中的四季变化。
  3. 循环成分(Cyclical Component, C [ t ] C[t] C[t]:表示时间序列中超过一年以上的周期性波动。
  4. 不规则成分(Irregular Component, I [ t ] I[t] I[t]:表示时间序列中的随机波动或噪声。

根据不同的组合方式,分解模型可以分为加法模型和乘法模型。

加法模型

在加法模型中,时间序列的各个成分是相互独立的,四个成分都有相同的量纲。加法模型的形式如下:

Y [ t ] = T [ t ] + S [ t ] + C [ t ] + I [ t ] Y[t] = T[t] + S[t] + C[t] + I[t] Y[t]=T[t]+S[t]+C[t]+I[t]

乘法模型

在乘法模型中,时间序列的输出部分和趋势项有相同的量纲,季节项和循环项是比例数,不规则变动项为独立随机变量序列,服从正态分布。乘法模型的形式如下:

Y [ t ] = T [ t ] × S [ t ] × C [ t ] × I [ t ] Y[t] = T[t] \times S[t] \times C[t] \times I[t] Y[t]=T[t]×S[t]×C[t]×I[t]

组合模型

当然,也可以将加法和乘法的分解模式进行组合,以便更好地描述复杂的时间序列。例如,某些情况下,可以使用以下形式的组合模型:

Y [ t ] = ( T [ t ] + S [ t ] ) × C [ t ] × I [ t ] Y[t] = (T[t] + S[t]) \times C[t] \times I[t] Y[t]=(T[t]+S[t])×C[t]×I[t]

Y [ t ] = T [ t ] × ( S [ t ] + C [ t ] ) + I [ t ] Y[t] = T[t] \times (S[t] + C[t]) + I[t] Y[t]=T[t]×(S[t]+C[t])+I[t]

通过对时间序列进行分解,我们可以更好地理解数据的内部结构,并使用适当的模型进行预测和分析。

加法模型代码示例
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# 生成示例时间序列数据
# 创建一个包含 120 个数据点的时间序列数据,包含趋势、季节性和不规则成分
np.random.seed(0)
n = 120
t = np.arange(n)
trend = 0.05 * t
seasonal = 10 * np.sin(2 * np.pi * t / 12)
irregular = np.random.normal(0, 1, n)
time_series = trend + seasonal + irregular

# 创建时间序列对象
# 使用 pandas 创建时间序列对象,并设定日期范围
date_range = pd.date_range(start='2020-01-01', periods=n, freq='M')
ts = pd.Series(time_series, index=date_range)

# 分解时间序列
# 使用 seasonal_decompose 函数对时间序列进行分解,选择加法模型(additive)
result = seasonal_decompose(ts, model='additive', period=12)

# 绘制分解结果
# 绘制分解后的趋势、季节性和残差成分
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 8))
result.observed.plot(ax=ax1)
ax1.set_ylabel('Observed')
result.trend.plot(ax=ax2)
ax2.set_ylabel('Trend')
result.seasonal.plot(ax=ax3)
ax3.set_ylabel('Seasonal')
result.resid.plot(ax=ax4)
ax4.set_ylabel('Residual')
plt.tight_layout()
plt.show()

# 组合模型示例
# 创建一个组合模型,将趋势、季节性和残差成分相加
combined_model = result.trend + result.seasonal + result.resid

# 绘制原始数据和组合模型的比较
plt.figure(figsize=(10, 4))
plt.plot(ts, label='Original Time Series')
plt.plot(combined_model, label='Combined Model', linestyle='--')
plt.legend()
plt.show()

在这里插入图片描述
在这里插入图片描述

乘法模型代码示例
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# 生成示例时间序列数据
np.random.seed(0)
n = 120
t = np.arange(n)
trend = 1 + 0.05 * t  # 确保趋势成分不为零或负值
seasonal = 1 + 0.1 * np.sin(2 * np.pi * t / 12)  # 确保季节性成分不为零或负值
irregular = 1 + 0.1 * np.random.normal(0, 1, n)  # 确保不规则成分不为零或负值
time_series = trend * seasonal * irregular

# 创建时间序列对象
date_range = pd.date_range(start='2020-01-01', periods=n, freq='M')
ts = pd.Series(time_series, index=date_range)

# 分解时间序列
result = seasonal_decompose(ts, model='multiplicative', period=12)

# 绘制分解结果
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 8))
result.observed.plot(ax=ax1)
ax1.set_ylabel('Observed')
result.trend.plot(ax=ax2)
ax2.set_ylabel('Trend')
result.seasonal.plot(ax=ax3)
ax3.set_ylabel('Seasonal')
result.resid.plot(ax=ax4)
ax4.set_ylabel('Residual')
plt.tight_layout()
plt.show()

# 组合模型示例
combined_model = result.trend * result.seasonal * result.resid

# 绘制原始数据和组合模型的比较
plt.figure(figsize=(10, 4))
plt.plot(ts, label='Original Time Series')
plt.plot(combined_model, label='Combined Model', linestyle='--')
plt.legend()
plt.show()

在这里插入图片描述
在这里插入图片描述

组合模型案例:零售销售数据

我们来看一个零售销售数据的案例,通过对销售额进行分解和组合建模,以便更好地理解和预测销售趋势。

假设我们有一个零售店的月度销售数据,数据包括了季节性变化(如假期和促销活动)、长期趋势(如整体销售增长)和一些不规则波动(如市场环境变化)。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose

# 生成假设的月度销售数据
np.random.seed(42)
n = 120
t = np.arange(n)
trend = 500 + 5 * t  # 假设一个线性增长的趋势
seasonal = 1 + 0.2 * np.sin(2 * np.pi * t / 12)  # 年周期性波动
cyclical = 1 + 0.1 * np.sin(2 * np.pi * t / 60)  # 长期周期性波动
irregular = 1 + 0.1 * np.random.normal(0, 1, n)  # 不规则波动

# 组合数据,确保没有零值或负值
sales_data = (trend + 100 * seasonal) * cyclical * irregular

# 创建时间序列对象
date_range = pd.date_range(start='2010-01-01', periods=n, freq='M')
ts = pd.Series(sales_data, index=date_range)

# 绘制原始数据
plt.figure(figsize=(10, 4))
plt.plot(ts, label='Monthly Sales Data')
plt.title('Monthly Sales Data')
plt.legend()
plt.show()

# 分解时间序列,使用加法模型来分解趋势和季节性
result = seasonal_decompose(ts, model='additive', period=12)

# 提取趋势和季节性成分
trend = result.trend
seasonal = result.seasonal

# 构建组合模型 Y[t] = (T[t] + S[t]) * C[t] * I[t]
combined_model = (trend + seasonal) * cyclical * irregular

# 绘制分解结果
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 8))
result.observed.plot(ax=ax1)
ax1.set_ylabel('Observed')
result.trend.plot(ax=ax2)
ax2.set_ylabel('Trend')
result.seasonal.plot(ax=ax3)
ax3.set_ylabel('Seasonal')
result.resid.plot(ax=ax4)
ax4.set_ylabel('Residual')
plt.tight_layout()
plt.show()

# 绘制原始数据和组合模型的比较
plt.figure(figsize=(10, 4))
plt.plot(ts, label='Original Sales Data')
plt.plot(combined_model, label='Combined Model', linestyle='--')
plt.legend()
plt.show()

在组合模型中,销售数据 Y [ t ] Y[t] Y[t] 被分解并重新组合为:

Y [ t ] = ( T [ t ] + S [ t ] ) × C [ t ] × I [ t ] Y[t] = (T[t] + S[t]) \times C[t] \times I[t] Y[t]=(T[t]+S[t])×C[t]×I[t]

  • 趋势成分(Trend Component):显示销售数据中的长期增长趋势。
  • 季节成分(Seasonal Component):显示销售数据中的年周期性波动,如假期和促销活动。
  • 循环成分(Cyclical Component):显示销售数据中的长期周期性波动,例如经济周期。
  • 不规则成分(Irregular Component):显示销售数据中的随机波动或噪声。

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

2. 移动平均法与指数平均法

2.1 移动平均法

移动平均法是一种时间序列分析的方法,通过取固定时间窗口内数据的平均值来平滑时间序列。它适用于平滑短期波动,突出长期趋势或周期。

一次移动平均

一次移动平均的计算公式是:
M ( 1 ) t = y t + y t − 1 + … + y t − N + 1 N M_{(1)}^t = \frac{y_t + y_{t-1} + \ldots + y_{t-N+1}}{N} M(1)t=Nyt+yt1++ytN+1

其中, y t y_t yt 是第 t t t 期的实际值, N N N 是移动平均的周期。

二次移动平均

二次移动平均是在一次移动平均的基础上再次进行移动平均,计算公式是:
M ( 2 ) t = M ( 1 ) t + M ( 1 ) t − 1 + … + M ( 1 ) t − N + 1 N M_{(2)}^t = \frac{M_{(1)}^t + M_{(1)}^{t-1} + \ldots + M_{(1)}^{t-N+1}}{N} M(2)t=NM(1)t+M(1)t1++M(1)tN+1

其中, M ( 1 ) t M_{(1)}^t M(1)t 是第 t t t 期的一次移动平均值。

代码示例
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 生成假设的月度销售数据
np.random.seed(42)
n = 120
t = np.arange(n)
sales_data = 100 + 0.5 * t + 5 * np.random.normal(size=n)  # 生成具有线性趋势的销售数据

# 创建时间序列对象
date_range = pd.date_range(start='2010-01-01', periods=n, freq='M')
ts = pd.Series(sales_data, index=date_range)

# 定义移动平均的周期
N = 12

# 计算一次移动平均
M1 = ts.rolling(window=N).mean()

# 计算二次移动平均
M2 = M1.rolling(window=N).mean()

# 绘制原始数据,一次移动平均和二次移动平均
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Sales Data')
plt.plot(M1, label='Once Moving Average')
plt.plot(M2, label='Twice Moving Average')
plt.legend()
plt.title('Original Sales Data and Moving Averages')
plt.show()

2.2 指数平均法

指数平均法(Exponential Smoothing)是一种加权移动平均方法,越近期的数据权重越大。它能够较好地应对时间序列中的趋势和季节性变化。

1)简单指数平滑

简单指数平滑的公式为:
S t = α y t + ( 1 − α ) S t − 1 S_t = \alpha y_t + (1 - \alpha) S_{t-1} St=αyt+(1α)St1

其中, S t S_t St 是第 t t t 期的平滑值, y t y_t yt 是第 t t t 期的实际值, α \alpha α 是平滑常数(0 < \alpha < 1)。

2)双指数平滑

双指数平滑考虑了数据的趋势,公式为:
S t = α y t + ( 1 − α ) ( S t − 1 + T t − 1 ) S_t = \alpha y_t + (1 - \alpha) (S_{t-1} + T_{t-1}) St=αyt+(1α)(St1+Tt1)
T t = β ( S t − S t − 1 ) + ( 1 − β ) T t − 1 T_t = \beta (S_t - S_{t-1}) + (1 - \beta) T_{t-1} Tt=β(StSt1)+(1β)Tt1

其中, T t T_t Tt 是第 t t t 期的趋势估计, β \beta β 是趋势平滑常数(0 < β \beta β < 1)。

代码示例
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# 生成假设的月度销售数据
np.random.seed(42)
n = 120
t = np.arange(n)
sales_data = 100 + 0.5 * t + 20 * np.sin(2 * np.pi * t / 12) + 5 * np.random.normal(size=n)  # 生成具有线性趋势和季节性的销售数据

# 创建时间序列对象
date_range = pd.date_range(start='2010-01-01', periods=n, freq='M')
ts = pd.Series(sales_data, index=date_range)

# 简单指数平滑
alpha = 0.2
simple_exp_smoothing = ts.ewm(alpha=alpha).mean()

# 双指数平滑
double_exp_smoothing = ExponentialSmoothing(ts, trend='add', seasonal=None).fit()
double_exp_smoothing_fitted = double_exp_smoothing.fittedvalues

# 绘制原始数据,简单指数平滑和双指数平滑
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Sales Data')
plt.plot(simple_exp_smoothing, label='Simple Exponential Smoothing')
plt.plot(double_exp_smoothing_fitted, label='Double Exponential Smoothing')
plt.legend()
plt.title('Original Sales Data and Exponential Smoothing')
plt.show()

在这里插入图片描述

代码解释
  1. 一次移动平均二次移动平均

    • 使用 rolling 方法计算一次和二次移动平均。
    • 一次移动平均平滑短期波动,二次移动平均进一步平滑一次移动平均的结果。
  2. 简单指数平滑

    • 使用 Pandas 中的 ewm 方法进行简单指数平滑。
    • 参数 alpha 控制权重分配,越接近1,越强调近期数据。
  3. 双指数平滑

    • 使用 statsmodels 中的 ExponentialSmoothing 方法进行双指数平滑。
    • trend='add' 参数表示添加趋势成分,用于捕捉时间序列中的线性趋势。

通过这些方法,可以有效地平滑时间序列数据,捕捉其趋势和季节性变化,从而更准确地进行预测。

3)三次指数平滑

三次指数平滑(Triple Exponential Smoothing),也称为霍尔特-温特斯(Holt-Winters)方法,是一种考虑了趋势和季节性成分的指数平滑方法。它适用于具有明显趋势和季节性变化的时间序列数据。

公式

三次指数平滑模型包含三个主要成分:水平(Level)、趋势(Trend)和季节性(Seasonal)。其递推公式如下:

  1. 水平成分(Level)
    L t = α y t S t − m + ( 1 − α ) ( L t − 1 + T t − 1 ) L_t = \alpha \frac{y_t}{S_{t-m}} + (1 - \alpha)(L_{t-1} + T_{t-1}) Lt=αStmyt+(1α)(Lt1+Tt1)
    其中, L t L_t Lt 是第 t t t 期的水平成分, y t y_t yt 是第 t t t 期的实际值, S t − m S_{t-m} Stm 是第 t − m t-m tm 期的季节性成分, α \alpha α 是平滑常数(0 < \alpha < 1)。

  2. 趋势成分(Trend)
    T t = β ( L t − L t − 1 ) + ( 1 − β ) T t − 1 T_t = \beta (L_t - L_{t-1}) + (1 - \beta) T_{t-1} Tt=β(LtLt1)+(1β)Tt1
    其中, T t T_t Tt 是第 t t t 期的趋势成分, β \beta β 是趋势平滑常数(0 < \beta < 1)。

  3. 季节性成分(Seasonal)
    S t = γ y t L t + ( 1 − γ ) S t − m S_t = \gamma \frac{y_t}{L_t} + (1 - \gamma) S_{t-m} St=γLtyt+(1γ)Stm
    其中, S t S_t St 是第 t t t 期的季节性成分, γ \gamma γ 是季节性平滑常数(0 < \gamma < 1), m m m 是季节周期长度。

  4. 预测值(Forecast)
    y ^ t + h = ( L t + h T t ) S t + h − m ( k + 1 ) \hat{y}_{t+h} = (L_t + h T_t) S_{t+h-m(k+1)} y^t+h=(Lt+hTt)St+hm(k+1)
    其中, y ^ t + h \hat{y}_{t+h} y^t+h 是第 t + h t+h t+h 期的预测值, h h h 是预测步长, k k k 是整数部分的商数 h m \frac{h}{m} mh

使用场景

三次指数平滑适用于以下情况:

  1. 时间序列具有明显的趋势和季节性:例如,零售销售数据、旅游行业的客流量等。
  2. 需要短期预测:三次指数平滑模型在短期预测中的表现通常优于长期预测。
  3. 需要较高的平滑效果:通过考虑趋势和季节性成分,三次指数平滑可以更好地捕捉数据的复杂变化。
代码示例

下面是使用 statsmodels 库实现三次指数平滑的示例代码:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# 生成假设的月度销售数据
np.random.seed(42)
n = 120
t = np.arange(n)
seasonal_pattern = np.sin(2 * np.pi * t / 12)
sales_data = 100 + 0.5 * t + 10 * seasonal_pattern + 5 * np.random.normal(size=n)  # 生成具有线性趋势和季节性的销售数据

# 创建时间序列对象
date_range = pd.date_range(start='2010-01-01', periods=n, freq='M')
ts = pd.Series(sales_data, index=date_range)

# 三次指数平滑
triple_exp_smoothing = ExponentialSmoothing(ts, trend='add', seasonal='add', seasonal_periods=12).fit()
triple_exp_smoothing_fitted = triple_exp_smoothing.fittedvalues

# 预测未来12个月的数据
forecast_period = 12
forecast = triple_exp_smoothing.forecast(forecast_period)

# 绘制原始数据和三次指数平滑结果
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Sales Data')
plt.plot(triple_exp_smoothing_fitted, label='Triple Exponential Smoothing (Fitted)')
plt.plot(forecast, label='Forecast', linestyle='--')
plt.legend()
plt.title('Original Sales Data and Triple Exponential Smoothing')
plt.show()

在这里插入图片描述

代码解释
  1. 生成假设的月度销售数据:创建一个具有线性趋势和季节性的时间序列数据。
  2. 三次指数平滑
    • 使用 ExponentialSmoothing 方法进行三次指数平滑,指定 trend='add'seasonal='add' 来捕捉数据中的线性趋势和加性季节性成分。
    • 调用 fit 方法进行模型拟合,得到平滑后的时间序列。
  3. 预测未来数据
    • 使用 forecast 方法预测未来 12 个月的数据。
  4. 绘图
    • 绘制原始数据、平滑后的数据和预测数据,以可视化三次指数平滑的效果。

3. ARIMA系列模型

3.1 AR模型 (Autoregressive Model)

自回归模型 (AR) 是时间序列分析中最简单和最常用的模型之一。AR模型假设时间序列的当前值是其过去若干值的线性组合,可以表示为:

Y t = ϕ 1 Y t − 1 + ϕ 2 Y t − 2 + ⋯ + ϕ p Y t − p + ϵ t Y_t = \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \cdots + \phi_p Y_{t-p} + \epsilon_t Yt=ϕ1Yt1+ϕ2Yt2++ϕpYtp+ϵt

其中:

  • Y t Y_t Yt 是时间序列在时间 t t t 的值
  • ϕ 1 , ϕ 2 , … , ϕ p \phi_1, \phi_2, \ldots, \phi_p ϕ1,ϕ2,,ϕp 是模型参数
  • p p p 是模型的阶数
  • ϵ t \epsilon_t ϵt 是白噪声误差项,通常假设其均值为0,方差为常数

Python中的代码实现:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.ar_model import AutoReg

# 生成假设的时间序列数据
np.random.seed(42)
n = 200
t = np.arange(n)
Y = 0.5 * np.random.normal(size=n)
for i in range(2, n):
    Y[i] = 1.5 * Y[i-1] - 0.7 * Y[i-2] + np.random.normal()

# 创建时间序列对象
ts = pd.Series(Y, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))

# 拟合AR模型
ar_model = AutoReg(ts, lags=2).fit()

# 预测
ar_pred = ar_model.predict(start=2, end=n-1)

# 绘制原始数据和预测值
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Time Series')
plt.plot(ts.index[2:], ar_pred, label='AR Model Predictions')
plt.legend()
plt.title('AR Model for Time Series Data')
plt.show()

在这里插入图片描述

3.2 MA模型 (Moving Average Model)

移动平均模型 (MA) 假设时间序列的当前值是其过去误差项的线性组合,可以表示为:

Y t = ϵ t + θ 1 ϵ t − 1 + θ 2 ϵ t − 2 + ⋯ + θ q ϵ t − q Y_t = \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q} Yt=ϵt+θ1ϵt1+θ2ϵt2++θqϵtq

其中:

  • Y t Y_t Yt 是时间序列在时间 t t t 的值
  • θ 1 , θ 2 , … , θ q \theta_1, \theta_2, \ldots, \theta_q θ1,θ2,,θq 是模型参数
  • q q q 是模型的阶数
  • ϵ t , ϵ t − 1 , … , ϵ t − q \epsilon_t, \epsilon_{t-1}, \ldots, \epsilon_{t-q} ϵt,ϵt1,,ϵtq 是白噪声误差项

Python中的代码实现:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA

# 生成假设的时间序列数据(具有MA过程的特性)
np.random.seed(42)
n = 200
theta1, theta2 = 0.5, 0.2
epsilon = np.random.normal(size=n)
Y = np.zeros(n)
for t in range(2, n):
    Y[t] = epsilon[t] + theta1 * epsilon[t-1] + theta2 * epsilon[t-2]

# 创建时间序列对象
ts = pd.Series(Y, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))

# 拟合MA模型
ma_model = ARIMA(ts, order=(0, 0, 2)).fit()

# 预测
ma_pred = ma_model.predict(start=0, end=n-1)

# 绘制原始数据和预测值
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Time Series')
plt.plot(ts.index, ma_pred, label='MA Model Predictions')
plt.legend()
plt.title('MA Model for Time Series Data')
plt.show()

在这里插入图片描述

3.3 ARMA模型和ARIMA模型

1) ARMA模型 (Autoregressive Moving Average Model)

ARMA模型结合了AR模型和MA模型的特性,用于描述平稳时间序列的数据,可以表示为:

Y t = ϕ 1 Y t − 1 + ϕ 2 Y t − 2 + ⋯ + ϕ p Y t − p + ϵ t + θ 1 ϵ t − 1 + θ 2 ϵ t − 2 + ⋯ + θ q ϵ t − q Y_t = \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \cdots + \phi_p Y_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q} Yt=ϕ1Yt1+ϕ2Yt2++ϕpYtp+ϵt+θ1ϵt1+θ2ϵt2++θqϵtq

Python中的代码实现:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA

# 生成假设的时间序列数据
np.random.seed(42)
n = 200
Y = np.random.normal(size=n)
for i in range(2, n):
    Y[i] = 1.5 * Y[i-1] - 0.7 * Y[i-2] + np.random.normal() + 0.5 * np.random.normal()

# 创建时间序列对象
ts = pd.Series(Y, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))

# 拟合ARMA模型
arma_model = ARIMA(ts, order=(2, 0, 1)).fit()

# 预测
arma_pred = arma_model.predict(start=2, end=n-1)

# 绘制原始数据和预测值
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Time Series')
plt.plot(ts.index[2:], arma_pred, label='ARMA Model Predictions')
plt.legend()
plt.title('ARMA Model for Time Series Data')
plt.show()

在这里插入图片描述

代码详解

在上面代码中,生成的数据模拟了一个 ARMA(2,1) 模型。这种模型结合了自回归 (AR) 和移动平均 (MA) 两部分。具体来说:

  • 自回归部分 (AR) 使用前两个时间点的值 Y [ i − 1 ] Y[i-1] Y[i1] Y [ i − 2 ] Y[i-2] Y[i2]
  • 移动平均部分 (MA) 使用当前和前一个时间点的随机噪声项。

更正式的 ARMA(2,1) 模型公式可以写成:

Y [ t ] = c + ϕ 1 Y [ t − 1 ] + ϕ 2 Y [ t − 2 ] + ϵ t + θ 1 ϵ t − 1 Y[t] = c + \phi_1 Y[t-1] + \phi_2 Y[t-2] + \epsilon_t + \theta_1 \epsilon_{t-1} Y[t]=c+ϕ1Y[t1]+ϕ2Y[t2]+ϵt+θ1ϵt1

其中:

  • Y [ t ] Y[t] Y[t] 是时间序列值。
  • c c c 是常数项(在你的例子中,忽略了常数项)。
  • ϕ 1 \phi_1 ϕ1 ϕ 2 \phi_2 ϕ2 是 AR 部分的系数,分别为 1.5 和 -0.7。
  • ϵ t \epsilon_t ϵt 是当前的误差项(噪声)。
  • θ 1 \theta_1 θ1 是 MA 部分的系数。

为了清楚地体现移动平均部分的公式,下面是一个更明确的解释和示例代码:

数据生成代码如下:

for i in range(2, n):
    Y[i] = 1.5 * Y[i-1] - 0.7 * Y[i-2] + np.random.normal() + 0.5 * np.random.normal()

这段代码生成的是包含AR(2)和MA(1)部分的时间序列数据。为了准确体现移动平均部分,我们可以对噪声项进行如下处理:

  • ϵ t \epsilon_t ϵt:当前时刻的误差项(噪声)。
  • ϵ t − 1 \epsilon_{t-1} ϵt1:前一个时刻的误差项(噪声)。

假设我们生成的数据是根据以下公式:

Y [ t ] = 1.5 Y [ t − 1 ] − 0.7 Y [ t − 2 ] + ϵ t + 0.5 ϵ t − 1 Y[t] = 1.5 Y[t-1] - 0.7 Y[t-2] + \epsilon_t + 0.5 \epsilon_{t-1} Y[t]=1.5Y[t1]0.7Y[t2]+ϵt+0.5ϵt1

为了更清晰地模拟这个过程,我们可以显式地存储误差项并在生成数据时使用:

np.random.seed(42)
n = 200
Y = np.zeros(n)
epsilon = np.random.normal(size=n)
for i in range(2, n):
    Y[i] = 1.5 * Y[i-1] - 0.7 * Y[i-2] + epsilon[i] + 0.5 * epsilon[i-1]

生成数据接下来便可以使用 statsmodels 库拟合 ARMA 模型,并进行预测

2) ARIMA模型 (Autoregressive Integrated Moving Average Model)

ARIMA模型扩展了ARMA模型,适用于非平稳时间序列。ARIMA模型通过差分操作将非平稳时间序列转化为平稳序列,然后再应用ARMA模型。其模型表示为:

Y t = ϕ 1 Y t − 1 + ϕ 2 Y t − 2 + ⋯ + ϕ p Y t − p + ϵ t + θ 1 ϵ t − 1 + θ 2 ϵ t − 2 + ⋯ + θ q ϵ t − q Y_t = \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \cdots + \phi_p Y_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q} Yt=ϕ1Yt1+ϕ2Yt2++ϕpYtp+ϵt+θ1ϵt1+θ2ϵt2++θqϵtq

其中,差分操作表示为 Δ Y t = Y t − Y t − 1 \Delta Y_t = Y_t - Y_{t-1} ΔYt=YtYt1

Python中的代码实现:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA

# 生成假设的时间序列数据
np.random.seed(42)
n = 200
t = np.arange(n)
Y = 0.5 * t + 5 * np.random.normal(size=n)

# 创建时间序列对象
ts = pd.Series(Y, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))

# 拟合ARIMA模型
arima_model = ARIMA(ts, order=(2, 1, 2)).fit()

# 预测
arima_pred = arima_model.predict(start=2, end=n-1)

# 绘制原始数据和预测值
plt.figure(figsize=(10, 6))
plt.plot(ts, label='Original Time Series')
plt.plot(ts.index[2:], arima_pred, label='ARIMA Model Predictions')
plt.legend()
plt.title('ARIMA Model for Time Series Data')
plt.show()
代码详解
  1. 数据生成与时间序列创建

    np.random.seed(42)
    n = 200
    t = np.arange(n)
    Y = 0.5 * t + 5 * np.random.normal(size=n)  # 生成具有线性趋势的时间序列数据
    
    • 这段代码生成了一个长度为200的时间序列数据,包含线性趋势和随机噪声。
    • np.random.seed(42) 用于设置随机种子,以保证结果可重复。
    • t 是一个从0到199的数组,用于模拟时间。
    • Y 是时间序列数据,具有线性趋势 0.5 * t 和随机噪声 5 * np.random.normal(size=n)
    ts = pd.Series(Y, index=pd.date_range(start='2020-01-01', periods=n, freq='D'))
    
    • pd.Series 创建一个时间序列对象,索引为从2020年1月1日开始的日期,频率为天 ('D')。
  2. 拟合ARIMA模型

    arima_model = ARIMA(ts, order=(2, 1, 2)).fit()
    
    • ARIMA 函数用于创建并拟合ARIMA模型。
    • order=(2, 1, 2) 参数说明了ARIMA模型的阶数:
      • p=2: 自回归项的阶数。
      • d=1: 差分次数,用于将非平稳序列转换为平稳序列。
      • q=2: 移动平均项的阶数。
  3. 预测与绘图

    arima_pred = arima_model.predict(start=2, end=n-1)
    
    • predict 函数用于生成从第2个数据点到第n-1个数据点的预测值。
    • start=2end=n-1 指定了预测范围。
    plt.figure(figsize=(10, 6))
    plt.plot(ts, label='Original Time Series')
    plt.plot(ts.index[2:], arima_pred, label='ARIMA Model Predictions')
    plt.legend()
    plt.title('ARIMA Model for Time Series Data')
    plt.show()
    
    • 使用 matplotlib 绘制原始时间序列数据和ARIMA模型预测值。
    • label 参数用于在图例中标识不同的数据。
    • plt.legend() 显示图例,plt.title() 设置图表标题。
  4. ARIMA模型的order参数设置的作用

  • p: 自回归 (AR) 项的阶数,表示模型中使用了多少个过去的值来预测当前值。
  • d: 差分次数,表示对原始数据进行几次差分操作以使其平稳。如果数据本身已经平稳,则 d=0。如果数据呈现一次趋势 (如线性趋势),则 d=1。如果数据呈现二次趋势,则 d=2,以此类推。
  • q: 移动平均 (MA) 项的阶数,表示模型中使用了多少个过去的预测误差项来预测当前值。

3.4 ARIMA模型order参数

ARIMA模型(AutoRegressive Integrated Moving Average Model,自回归积分滑动平均模型)是一种广泛用于时间序列分析和预测的统计模型。它结合了自回归 (AR) 和移动平均 (MA) 组件,并引入了差分以使时间序列平稳。ARIMA模型的order参数由三个部分组成:(p, d, q),分别表示自回归阶数、差分阶数和移动平均阶数。

1) 自回归 (AR) 部分 - p

自回归部分表示时间序列自身滞后值的线性组合。p表示模型中使用的滞后项的数量。

公式:
Y t = c + ∑ i = 1 p ϕ i Y t − i + ϵ t Y_t = c + \sum_{i=1}^{p} \phi_i Y_{t-i} + \epsilon_t Yt=c+i=1pϕiYti+ϵt

其中:

  • Y t Y_t Yt 是时间序列值
  • c c c 是常数项
  • ϕ i \phi_i ϕi 是自回归系数
  • ϵ t \epsilon_t ϵt 是白噪声误差项
2) 差分 (I) 部分 - d

差分部分用于使非平稳时间序列变得平稳。d表示进行几次差分以使时间序列平稳。

公式:
Y t ′ = Y t − Y t − 1 Y_t' = Y_t - Y_{t-1} Yt=YtYt1

如果d=1,表示进行一次差分。如果仍然不平稳,可能需要进行更多次差分。

3) 移动平均 (MA) 部分 - q

移动平均部分表示过去白噪声误差项的线性组合。q表示模型中使用的误差项的数量。

公式:
Y t = c + ϵ t + ∑ i = 1 q θ i ϵ t − i Y_t = c + \epsilon_t + \sum_{i=1}^{q} \theta_i \epsilon_{t-i} Yt=c+ϵt+i=1qθiϵti

其中:

  • Y t Y_t Yt 是时间序列值
  • c c c 是常数项
  • θ i \theta_i θi 是移动平均系数
  • ϵ t \epsilon_t ϵt 是白噪声误差项
ARIMA模型公式

将上述三个部分结合起来,即可得到ARIMA模型:

Y t = c + ∑ i = 1 p ϕ i Y t − i + ϵ t + ∑ j = 1 q θ j ϵ t − j Y_t = c + \sum_{i=1}^{p} \phi_i Y_{t-i} + \epsilon_t + \sum_{j=1}^{q} \theta_j \epsilon_{t-j} Yt=c+i=1pϕiYti+ϵt+j=1qθjϵtj

设置order参数的作用
  1. 自回归阶数 ( p ):确定多少个过去的时间序列值用于预测当前值。如果时间序列中存在显著的自相关,较高的p值可能更合适。
  2. 差分阶数 (d):确定需要进行几次差分以使时间序列平稳。如果时间序列存在趋势或不平稳,可能需要进行一次或多次差分。
  3. 移动平均阶数 (q):确定多少个过去的误差项用于预测当前值。如果时间序列中的噪声具有显著的自相关,较高的q值可能更合适。

4. 灰色系统模型

灰色系统理论是由我国学者邓聚龙教授于1982年提出的。灰色系统理论主要解决系统中由于信息不完全、不确定性较大而导致的研究难题,特别适用于数据较少、信息不完全的系统。灰色系统模型主要包括灰色预测模型和灰色关联模型。

4.1 灰色预测模型

灰色预测模型(GM, Grey Model)主要用于时间序列数据的预测。最常用的是GM(1,1)模型。GM(1,1)模型是一阶单变量模型,常用于短期预测。

GM(1,1)模型公式推导
  1. 数据生成:对原始数据序列 X ( 0 ) = { x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , … , x ( 0 ) ( n ) } X^{(0)} = \{x^{(0)}(1), x^{(0)}(2), \ldots, x^{(0)}(n)\} X(0)={x(0)(1),x(0)(2),,x(0)(n)} 进行一次累加生成,得到累加生成序列 X ( 1 ) X^{(1)} X(1)
    x ( 1 ) ( k ) = ∑ i = 1 k x ( 0 ) ( i ) , k = 1 , 2 , … , n x^{(1)}(k) = \sum_{i=1}^{k} x^{(0)}(i), \quad k = 1, 2, \ldots, n x(1)(k)=i=1kx(0)(i),k=1,2,,n

  2. 建立灰色微分方程:累加生成序列 X ( 1 ) X^{(1)} X(1) 满足如下的灰色微分方程:
    d x ( 1 ) ( t ) d t + a x ( 1 ) ( t ) = b \frac{dx^{(1)}(t)}{dt} + a x^{(1)}(t) = b dtdx(1)(t)+ax(1)(t)=b
    其中, a a a b b b 为待估参数。

  3. 离散化:将灰色微分方程离散化,得到:
    x ( 0 ) ( k ) + a z ( 1 ) ( k ) = b , k = 2 , 3 , … , n x^{(0)}(k) + a z^{(1)}(k) = b, \quad k = 2, 3, \ldots, n x(0)(k)+az(1)(k)=b,k=2,3,,n
    其中, z ( 1 ) ( k ) z^{(1)}(k) z(1)(k) 是累加生成序列 x ( 1 ) ( k ) x^{(1)}(k) x(1)(k) 的均值生成序列:
    z ( 1 ) ( k ) = 1 2 ( x ( 1 ) ( k ) + x ( 1 ) ( k − 1 ) ) z^{(1)}(k) = \frac{1}{2} (x^{(1)}(k) + x^{(1)}(k-1)) z(1)(k)=21(x(1)(k)+x(1)(k1))

  4. 参数估计:根据最小二乘法估计参数 a a a b b b
    θ ^ = [ a b ] = ( B T B ) − 1 B T Y \hat{\theta} = \begin{bmatrix} a \\ b \end{bmatrix} = (B^T B)^{-1} B^T Y θ^=[ab]=(BTB)1BTY
    其中,
    B = [ − z ( 1 ) ( 2 ) 1 − z ( 1 ) ( 3 ) 1 ⋮ ⋮ − z ( 1 ) ( n ) 1 ] , Y = [ x ( 0 ) ( 2 ) x ( 0 ) ( 3 ) ⋮ x ( 0 ) ( n ) ] B = \begin{bmatrix} -z^{(1)}(2) & 1 \\ -z^{(1)}(3) & 1 \\ \vdots & \vdots \\ -z^{(1)}(n) & 1 \end{bmatrix}, \quad Y = \begin{bmatrix} x^{(0)}(2) \\ x^{(0)}(3) \\ \vdots \\ x^{(0)}(n) \end{bmatrix} B= z(1)(2)z(1)(3)z(1)(n)111 ,Y= x(0)(2)x(0)(3)x(0)(n)

  5. 还原:通过累加生成序列 X ( 1 ) X^{(1)} X(1) 预测原始序列 X ( 0 ) X^{(0)} X(0),得到预测模型:
    x ( 1 ) ( k + 1 ) = ( x ( 1 ) ( 1 ) − b a ) e − a k + b a x^{(1)}(k+1) = \left( x^{(1)}(1) - \frac{b}{a} \right) e^{-a k} + \frac{b}{a} x(1)(k+1)=(x(1)(1)ab)eak+ab
    通过累减还原,得到对原始数据的预测值:
    x ( 0 ) ( k + 1 ) = x ( 1 ) ( k + 1 ) − x ( 1 ) ( k ) x^{(0)}(k+1) = x^{(1)}(k+1) - x^{(1)}(k) x(0)(k+1)=x(1)(k+1)x(1)(k)

GM(1,1)模型的Python代码实现

以下是一个使用Python实现GM(1,1)模型进行时间序列预测的示例:

import numpy as np
import pandas as pd

def GM11(x0):
    x1 = np.cumsum(x0)
    z1 = (x1[:-1] + x1[1:]) / 2.0
    B = np.vstack([-z1, np.ones_like(z1)]).T
    Y = x0[1:]
    [[a], [b]] = np.linalg.inv(B.T @ B) @ B.T @ Y
    def predict(k):
        return (x0[0] - b / a) * np.exp(-a * k) + b / a
    return predict

# 示例数据
data = np.array([7.94, 8.15, 8.06, 8.05, 7.60, 7.39, 7.74, 7.38, 8.32, 7.69, 8.15, 6.88, 8.00, 7.94, 8.01, 7.91, 8.04])

# 建立GM(1,1)模型
gm = GM11(data)

# 预测未来5个时刻的值
predictions = [gm(k) for k in range(len(data), len(data) + 5)]
print("未来5个时刻的预测值:", predictions)

4.2 灰色关联模型

灰色关联模型用于测度不同系统或同一系统内部不同因素之间的关联程度。其基本思想是通过比较各因素的变化趋势来确定它们之间的关联程度。灰色关联模型常用于多指标评价和决策分析。

灰色关联模型公式推导
  1. 确定参考序列和比较序列:设参考序列为 X 0 = { x 0 ( 1 ) , x 0 ( 2 ) , … , x 0 ( n ) } X_0 = \{x_0(1), x_0(2), \ldots, x_0(n)\} X0={x0(1),x0(2),,x0(n)},比较序列为 X i = { x i ( 1 ) , x i ( 2 ) , … , x i ( n ) } X_i = \{x_i(1), x_i(2), \ldots, x_i(n)\} Xi={xi(1),xi(2),,xi(n)},其中 i = 1 , 2 , … , m i = 1, 2, \ldots, m i=1,2,,m

  2. 求差序列:计算参考序列与比较序列的差序列:
    Δ i ( k ) = ∣ x 0 ( k ) − x i ( k ) ∣ \Delta_i(k) = |x_0(k) - x_i(k)| Δi(k)=x0(k)xi(k)

  3. 求最大差和最小差
    Δ max = max ⁡ i , k Δ i ( k ) , Δ min = min ⁡ i , k Δ i ( k ) \Delta_{\text{max}} = \max_{i,k} \Delta_i(k), \quad \Delta_{\text{min}} = \min_{i,k} \Delta_i(k) Δmax=i,kmaxΔi(k),Δmin=i,kminΔi(k)

  4. 计算关联系数:计算每个时刻的关联系数:
    ξ i ( k ) = Δ min + ρ Δ max Δ i ( k ) + ρ Δ max \xi_i(k) = \frac{\Delta_{\text{min}} + \rho \Delta_{\text{max}}}{\Delta_i(k) + \rho \Delta_{\text{max}}} ξi(k)=Δi(k)+ρΔmaxΔmin+ρΔmax
    其中, ρ \rho ρ 为分辨系数,通常取值为 0.5。

  5. 计算关联度:计算比较序列与参考序列的关联度:
    r i = 1 n ∑ k = 1 n ξ i ( k ) r_i = \frac{1}{n} \sum_{k=1}^{n} \xi_i(k) ri=n1k=1nξi(k)

灰色关联模型的Python代码实现

以下是一个使用Python实现灰色关联分析的示例:

import numpy as np

def grey_relation(reference, comparables, rho=0.5):
    reference = np.array(reference)
    comparables = np.array(comparables)
    n = len(reference)
    
    delta = np.abs(comparables - reference[:, None])
    delta_min = np.min(delta)
    delta_max = np.max(delta)
    
    xi = (delta_min + rho * delta_max) / (delta + rho * delta_max)
    r = xi.mean(axis=0)
    return r

# 示例数据
reference = [7.94, 8.15, 8.06, 8.05, 7.60, 7.39, 7.74, 7.38, 8.32, 7.69, 8.15, 6.88, 8.00, 7.94, 8.01, 7.91, 8.04]
comparables = [
    [9.47, 9.00, 8.45, 9.16, 7.93, 7.12, 7.29, 6.51, 8.47, 8.50, 9.88, 7.59, 8.15, 7.48, 7.76, 7.93, 8.34],
    [1.63, 1.40, 2.83, 3.33, 2.07, 2.23, 1.77, 3.63, 1.60, 2.73, 2.00, 1.77, 4.87, 3.30, 2.67, 5.47, 3.87],
    [0.08, 0.42, 0.20, 0.29, 0.13, 0.20, 0.

06, 0.41, 0.14, 0.28, 0.08, 0.92, 0.33, 0.13, 6.36, 0.21, 0.20],
    [6.78, 5.27, 2.50, 5.65, 5.26, 6.21, 6.44, 3.17, 3.32, 7.25, 7.22, 2.94, 4.68, 5.81, 4.52, 7.48, 3.73],
    [4.94, 4.47, 7.03, 5.26, 6.39, 7.50, 3.93, 5.75, 6.29, 8.21, 3.82, 8.03, 5.01, 4.87, 3.22, 2.40, 4.05]
]

# 计算关联度
relation_degrees = grey_relation(reference, comparables)
print("各指标与参考序列的关联度:", relation_degrees)

上述代码实现了灰色关联分析,计算了各比较序列与参考序列之间的关联度。关联度越大,表示该比较序列与参考序列的变化趋势越一致。

灰色关联模型的应用

灰色关联模型(Grey Relational Analysis, GRA)在测度关联程度方面有广泛的应用,特别是在多指标评价和决策分析中。以下是一些常见的应用领域和具体案例:

应用领域
  1. 经济管理

    • 企业绩效评价:通过灰色关联分析对多个企业的绩效指标进行综合评价,找出最优企业。
    • 投资风险分析:评估不同投资项目的风险与收益,选择最佳投资方案。
  2. 工程与技术

    • 产品质量评价:对比多个产品的质量指标,确定质量最佳的产品。
    • 工艺参数优化:在制造过程中,优化工艺参数以提高产品质量和生产效率。
  3. 环境科学

    • 环境质量评价:分析不同地区的环境监测数据,评价环境质量。
    • 污染源识别:识别和评价污染源对环境质量的影响。
  4. 社会科学

    • 教育质量评价:对比不同学校或教育机构的教学质量指标,找出最优的教育方案。
    • 社会经济发展评价:评估不同地区的经济社会发展水平,制定相应的发展政策。
  5. 医学与健康

    • 疾病诊断与治疗:分析患者的多项健康指标,辅助医生进行疾病诊断和治疗方案选择。
    • 药物效果评价:比较不同药物的治疗效果,选择最优药物。
案例一:企业绩效评价

假设有三家企业A、B、C,它们的绩效评价指标包括销售额、利润率、市场占有率和员工满意度。通过灰色关联分析,可以评估每个企业在这些指标上的综合表现,从而得出哪家企业的整体绩效最佳。

import numpy as np

# 参考序列:理想的绩效指标
reference = [100, 20, 30, 90]

# 比较序列:企业A、B、C的绩效指标
comparables = [
    [95, 18, 28, 85],   # 企业A
    [90, 22, 26, 88],   # 企业B
    [88, 19, 29, 87]    # 企业C
]

def grey_relation(reference, comparables, rho=0.5):
    reference = np.array(reference)
    comparables = np.array(comparables)
    n = len(reference)
    
    delta = np.abs(comparables - reference[:, None])
    delta_min = np.min(delta)
    delta_max = np.max(delta)
    
    xi = (delta_min + rho * delta_max) / (delta + rho * delta_max)
    r = xi.mean(axis=0)
    return r

# 计算关联度
relation_degrees = grey_relation(reference, comparables)
print("各企业与理想绩效指标的关联度:", relation_degrees)

通过计算得到每个企业的关联度,从而可以评估哪个企业的绩效与理想状态最接近。

案例二:环境质量评价

假设有三个城市X、Y、Z,它们的环境质量指标包括空气质量指数、水质指数、噪声指数和绿化覆盖率。通过灰色关联分析,可以评估每个城市在这些指标上的环境质量,从而得出哪个城市的环境质量最佳。

# 参考序列:理想的环境质量指标
reference = [50, 90, 30, 70]

# 比较序列:城市X、Y、Z的环境质量指标
comparables = [
    [55, 85, 35, 75],   # 城市X
    [52, 88, 32, 72],   # 城市Y
    [48, 92, 28, 68]    # 城市Z
]

# 计算关联度
relation_degrees = grey_relation(reference, comparables)
print("各城市与理想环境质量指标的关联度:", relation_degrees)

通过计算得到每个城市的关联度,从而可以评估哪个城市的环境质量与理想状态最接近。

案例三:药物效果评价

假设有三种药物A、B、C,它们在治疗某种疾病时的效果指标包括治愈率、副作用发生率、患者满意度和恢复时间。通过灰色关联分析,可以评估每种药物的综合效果,从而选择最优的药物。

# 参考序列:理想的药物效果指标
reference = [95, 5, 90, 10]

# 比较序列:药物A、B、C的效果指标
comparables = [
    [92, 6, 85, 12],   # 药物A
    [90, 7, 88, 11],   # 药物B
    [93, 5, 89, 9]     # 药物C
]

# 计算关联度
relation_degrees = grey_relation(reference, comparables)
print("各药物与理想效果指标的关联度:", relation_degrees)

通过计算得到每种药物的关联度,从而可以评估哪个药物的效果与理想状态最接近。

灰色关联模型通过定量分析多指标数据之间的关联程度,能够为决策提供科学依据,广泛应用于各个领域的评价与决策分析中。

5. 组合投资问题的一些策略

组合投资策略是通过分散投资来降低风险并优化回报的金融方法。以下是三种常见的组合投资策略:马科维兹均值方差模型、最大夏普比率模型和风险平价模型。

5.1 马科维兹均值方差模型

马科维兹均值方差模型(Markowitz Mean-Variance Model)是现代投资组合理论的基础,旨在通过优化投资组合的权重来平衡预期收益和风险。

公式推导
  1. 预期收益率:假设有 n n n 个资产,每个资产的预期收益率为 μ i \mu_i μi,资产权重为 w i w_i wi,投资组合的预期收益率为:
    μ p = ∑ i = 1 n w i μ i \mu_p = \sum_{i=1}^{n} w_i \mu_i μp=i=1nwiμi

  2. 风险计算:投资组合的风险用方差表示,资产收益的协方差矩阵为 Σ \Sigma Σ,投资组合的方差为:
    σ p 2 = ∑ i = 1 n ∑ j = 1 n w i w j σ i j = w T Σ w \sigma_p^2 = \sum_{i=1}^{n} \sum_{j=1}^{n} w_i w_j \sigma_{ij} = w^T \Sigma w σp2=i=1nj=1nwiwjσij=wTΣw

  3. 目标函数:通过优化投资组合的权重 w i w_i wi 来最小化风险(方差),目标函数为:
    min ⁡ w w T Σ w \min_w w^T \Sigma w wminwTΣw
    约束条件为:
    ∑ i = 1 n w i μ i = μ p , ∑ i = 1 n w i = 1 , w i ≥ 0 \sum_{i=1}^{n} w_i \mu_i = \mu_p, \quad \sum_{i=1}^{n} w_i = 1, \quad w_i \geq 0 i=1nwiμi=μp,i=1nwi=1,wi0

代码示例
import numpy as np
import scipy.optimize as sco

# 示例数据:资产的预期收益率和协方差矩阵
returns = np.array([0.12, 0.18, 0.15])
cov_matrix = np.array([
    [0.1, 0.02, 0.04],
    [0.02, 0.08, 0.02],
    [0.04, 0.02, 0.09]
])

# 目标收益率
target_return = 0.15

def portfolio_variance(weights, cov_matrix):
    return weights.T @ cov_matrix @ weights

def constraint_return(weights, returns, target_return):
    return weights.T @ returns - target_return

def constraint_sum(weights):
    return np.sum(weights) - 1

# 约束条件
constraints = (
    {'type': 'eq', 'fun': constraint_return, 'args': (returns, target_return)},
    {'type': 'eq', 'fun': constraint_sum}
)

# 初始猜测
init_guess = np.array([1/3, 1/3, 1/3])

# 边界条件
bounds = tuple((0, 1) for _ in range(returns.shape[0]))

# 优化
opt_result = sco.minimize(portfolio_variance, init_guess, args=(cov_matrix,), method='SLSQP', constraints=constraints, bounds=bounds)

# 最优权重
optimal_weights = opt_result.x
print("最优权重:", optimal_weights)

5.2 最大夏普比率模型

夏普比率(Sharpe Ratio)是衡量投资组合单位风险所获得收益的指标,定义为投资组合的超额收益与其标准差的比值。最大夏普比率模型通过最大化夏普比率来优化投资组合。

公式推导
  1. 夏普比率:投资组合的夏普比率定义为:
    S = μ p − r f σ p S = \frac{\mu_p - r_f}{\sigma_p} S=σpμprf
    其中, μ p \mu_p μp 是投资组合的预期收益率, r f r_f rf 是无风险收益率, σ p \sigma_p σp 是投资组合的标准差。

  2. 目标函数:通过优化投资组合的权重 w i w_i wi 来最大化夏普比率,目标函数为:
    max ⁡ w w T μ − r f w T Σ w \max_w \frac{w^T \mu - r_f}{\sqrt{w^T \Sigma w}} wmaxwTΣw wTμrf
    约束条件为:
    ∑ i = 1 n w i = 1 , w i ≥ 0 \sum_{i=1}^{n} w_i = 1, \quad w_i \geq 0 i=1nwi=1,wi0

代码示例
import numpy as np
import scipy.optimize as sco

# 示例数据:资产的预期收益率、协方差矩阵和无风险收益率
returns = np.array([0.12, 0.18, 0.15])
cov_matrix = np.array([
    [0.1, 0.02, 0.04],
    [0.02, 0.08, 0.02],
    [0.04, 0.02, 0.09]
])
risk_free_rate = 0.03

def sharpe_ratio(weights, returns, cov_matrix, risk_free_rate):
    portfolio_return = weights.T @ returns
    portfolio_std = np.sqrt(weights.T @ cov_matrix @ weights)
    return (portfolio_return - risk_free_rate) / portfolio_std

# 负夏普比率(因为我们需要最大化,所以取负数进行最小化)
def neg_sharpe_ratio(weights, returns, cov_matrix, risk_free_rate):
    return -sharpe_ratio(weights, returns, cov_matrix, risk_free_rate)

# 约束条件
constraints = {'type': 'eq', 'fun': lambda weights: np.sum(weights) - 1}

# 初始猜测
init_guess = np.array([1/3, 1/3, 1/3])

# 边界条件
bounds = tuple((0, 1) for _ in range(returns.shape[0]))

# 优化
opt_result = sco.minimize(neg_sharpe_ratio, init_guess, args=(returns, cov_matrix, risk_free_rate), method='SLSQP', constraints=constraints, bounds=bounds)

# 最优权重
optimal_weights = opt_result.x
print("最优权重:", optimal_weights)

5.3 风险平价模型

风险平价模型(Risk Parity Model)是一种投资组合优化策略,其目标是使投资组合中每种资产的风险贡献相等。风险平价模型在降低组合风险的同时追求稳健的收益。

公式推导
  1. 风险贡献:资产 i i i 的风险贡献定义为:
    R C i = w i ∂ σ p ∂ w i RC_i = w_i \frac{\partial \sigma_p}{\partial w_i} RCi=wiwiσp
    其中, ∂ σ p ∂ w i \frac{\partial \sigma_p}{\partial w_i} wiσp 是资产 i i i 对投资组合标准差的边际贡献。

  2. 边际风险贡献:资产 i i i 的边际风险贡献为:
    M R C i = ∂ σ p ∂ w i = ( Σ w ) i σ p MRC_i = \frac{\partial \sigma_p}{\partial w_i} = \frac{(\Sigma w)_i}{\sigma_p} MRCi=wiσp=σp(Σw)i

  3. 风险平价条件:风险平价模型的目标是使每种资产的风险贡献相等,即:
    R C i = R C j , ∀ i , j RC_i = RC_j, \quad \forall i, j RCi=RCj,i,j

代码示例
import numpy as np
import scipy.optimize as sco

# 示例数据:资产的协方差矩阵
cov_matrix = np.array([
    [0.1, 0.02, 0.04],
    [0.02, 0.08, 0.02],
    [0.04, 0.02, 0.09]
])

def portfolio_std(weights, cov_matrix):
    return np.sqrt(weights.T @ cov_matrix @ weights)

def risk_contribution(weights, cov_matrix):
    portfolio_std_val = portfolio_std(weights, cov_matrix)
    mrc = cov_matrix @ weights / portfolio_std_val
    rc = weights * mrc
    return rc

def risk_parity_objective(weights, cov_matrix):
    rc = risk_contribution(weights, cov_matrix)
    target = np.mean(rc)
    return np.sum((rc - target) ** 2)

# 约束条件
constraints = {'type': 'eq', 'fun': lambda weights: np.sum(weights) - 1}

# 初始猜测
init_guess = np.array([1/3, 1/3, 1/3])

# 边界条件
bounds = tuple((0, 1) for _ in range(cov_matrix.shape[0]))

# 优化
opt_result = sco.minimize(risk_parity_objective, init_guess, args=(cov_matrix,), method='SLSQP', constraints=constraints, bounds=bounds)

# 最优权重
optimal_weights = opt_result.x
print("最优权重:", optimal_weights)

6. 马尔可夫模型

马尔可夫模型(Markov Model)是一类用于描述随机过程的数学模型,在许多领域(如语音识别、自然语言处理、金融建模等)有广泛应用。马尔可夫模型基于马尔可夫性质,即未来状态只依赖于当前状态,而与过去状态无关。

6.1 一些马尔可夫模型的概念

马尔可夫链

马尔可夫链(Markov Chain)是马尔可夫模型的基本形式,用于描述系统在不同状态之间的转移。马尔可夫链由以下要素组成:

  1. 状态空间 S = { s 1 , s 2 , . . . , s n } S = \{s_1, s_2, ..., s_n\} S={s1,s2,...,sn}:系统可能的状态集合。
  2. 转移概率矩阵 P = [ p i j ] P = [p_{ij}] P=[pij]:从状态 s i s_i si 转移到状态 s j s_j sj 的概率,满足 ∑ j p i j = 1 \sum_{j} p_{ij} = 1 jpij=1
隐马尔可夫模型

隐马尔可夫模型(Hidden Markov Model, HMM)是马尔可夫链的扩展,适用于无法直接观察状态的情况。HMM由以下要素组成:

  1. 状态空间 S = { s 1 , s 2 , . . . , s n } S = \{s_1, s_2, ..., s_n\} S={s1,s2,...,sn}
  2. 观察空间 O = { o 1 , o 2 , . . . , o m } O = \{o_1, o_2, ..., o_m\} O={o1,o2,...,om}:可能的观察结果集合。
  3. 初始概率分布 π = { π i } \pi = \{\pi_i\} π={πi}:初始时刻在各状态的概率。
  4. 状态转移概率矩阵 A = [ a i j ] A = [a_{ij}] A=[aij]:从状态 s i s_i si 转移到状态 s j s_j sj 的概率。
  5. 观测概率矩阵 B = [ b j k ] B = [b_{jk}] B=[bjk]:在状态 s j s_j sj 下观察结果为 o k o_k ok 的概率。

6.2 马尔可夫模型的实现

马尔可夫模型的实现可以分为三类主要问题:评估问题、解码问题和学习问题。

1)评估问题

评估问题是指给定模型参数( π , A , B \pi, A, B π,A,B)和观察序列 O = { o 1 , o 2 , . . . , o T } O = \{o_1, o_2, ..., o_T\} O={o1,o2,...,oT},计算观察序列出现的概率。常用前向算法(Forward Algorithm)解决该问题。

前向算法
α t ( i ) = P ( o 1 , o 2 , . . . , o t , s t = s i ∣ λ ) \alpha_t(i) = P(o_1, o_2, ..., o_t, s_t = s_i \mid \lambda) αt(i)=P(o1,o2,...,ot,st=siλ)
其中, α t ( i ) \alpha_t(i) αt(i) 表示在时刻 t t t 观察到 o 1 , o 2 , . . . , o t o_1, o_2, ..., o_t o1,o2,...,ot 且处于状态 s i s_i si 的概率。

递推公式:
α t + 1 ( j ) = [ ∑ i = 1 N α t ( i ) a i j ] b j ( o t + 1 ) \alpha_{t+1}(j) = \left[ \sum_{i=1}^{N} \alpha_t(i) a_{ij} \right] b_j(o_{t+1}) αt+1(j)=[i=1Nαt(i)aij]bj(ot+1)

初始化:
α 1 ( i ) = π i b i ( o 1 ) \alpha_1(i) = \pi_i b_i(o_1) α1(i)=πibi(o1)

终止:
P ( O ∣ λ ) = ∑ i = 1 N α T ( i ) P(O \mid \lambda) = \sum_{i=1}^{N} \alpha_T(i) P(Oλ)=i=1NαT(i)

代码示例
import numpy as np

# 定义模型参数
states = {'Rainy': 0, 'Sunny': 1}
observations = ['walk', 'shop', 'clean']
start_probability = np.array([0.6, 0.4])
transition_probability = np.array([[0.7, 0.3], [0.4, 0.6]]) # 转移概率
emission_probability = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]]) #

# 观测序列,用索引表示
obs_seq = [0, 1, 2]  # 对应于 'walk', 'shop', 'clean'

# 初始化前向概率矩阵
alpha = np.zeros((len(obs_seq), len(states)))

# 初始化
alpha[0, :] = start_probability * emission_probability[:, obs_seq[0]]

# 递推计算
for t in range(1, len(obs_seq)):
    for j in range(len(states)):
        alpha[t, j] = np.dot(alpha[t-1, :], transition_probability[:, j]) * emission_probability[j, obs_seq[t]]

# 序列的总概率为最后一步的概率之和
total_prob = np.sum(alpha[-1, :])

print("Forward Probability Matrix:")
print(alpha)
print("\nTotal Probability of Observations:", total_prob)

定义模型参数
states = {'Rainy': 0, 'Sunny': 1}
observations = ['walk', 'shop', 'clean']
start_probability = np.array([0.6, 0.4])
transition_probability = np.array([[0.7, 0.3], [0.4, 0.6]])  # 转移概率
emission_probability = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
  • states:定义了状态空间,包括“Rainy”和“Sunny”,并给它们分配了索引 0 和 1。
  • observations:定义了观测空间,包括“walk”(散步)、“shop”(购物)和“clean”(打扫)。
  • start_probability:初始状态概率向量,表示初始时刻处于每个状态的概率。
  • transition_probability:状态转移概率矩阵,表示从一个状态转移到另一个状态的概率。
  • emission_probability:观测概率矩阵,表示在某个状态下观察到某个观测的概率。
观测序列
obs_seq = [0, 1, 2]  # 对应于 'walk', 'shop', 'clean'
  • obs_seq:给定的观测序列,这里使用索引表示,分别对应于“walk”、“shop”和“clean”。
初始化前向概率矩阵
alpha = np.zeros((len(obs_seq), len(states)))
  • alpha:前向概率矩阵,大小为(观测序列长度,状态数)。alpha[t, j]表示在时刻 t 处于状态 j 并观测到序列前 t 个观测的概率。
初始化
alpha[0, :] = start_probability * emission_probability[:, obs_seq[0]]
  • 初始化前向概率矩阵的第一行。alpha[0, j]表示在时刻 0 处于状态 j 并观测到第一个观测的概率,计算公式为:
    α [ 0 , j ] = start_probability [ j ] × emission_probability [ j , obs_seq [ 0 ] ] \alpha[0, j] = \text{start\_probability}[j] \times \text{emission\_probability}[j, \text{obs\_seq}[0]] α[0,j]=start_probability[j]×emission_probability[j,obs_seq[0]]

隐马尔可夫模型的定义
在隐马尔可夫模型中,有两个主要的概率分布:

  1. 初始状态概率(Initial State Probability):表示在时刻 0 处于每个状态的概率。记为 π i \pi_i πi,其中 i i i 表示状态。
  2. 观测概率(Emission Probability):表示在某个状态下观察到某个观测值的概率。记为 B j ( o t ) B_j(o_t) Bj(ot),其中 j j j 表示状态, o t o_t ot 表示在时刻 t t t 的观测值。

初始化前向概率

在前向算法中,初始化步骤是计算在时刻 0 处于每个状态并观测到第一个观测值的概率。这一步可以分解为两个部分:

  1. 初始状态概率:在时刻 0 处于状态 j j j 的概率,记为 π j \pi_j πj
  2. 观测概率:在状态 j j j 下观测到第一个观测值的概率,记为 B j ( o 0 ) B_j(o_0) Bj(o0),其中 o 0 o_0 o0 是第一个观测值。

具体公式

因此, α 0 ( j ) \alpha_0(j) α0(j) 表示在时刻 0 处于状态 j j j 并观测到第一个观测值的概率,其计算公式为:

α 0 ( j ) = π j × B j ( o 0 ) \alpha_0(j) = \pi_j \times B_j(o_0) α0(j)=πj×Bj(o0)

对应代码

在代码中,这一计算步骤对应为:

alpha[0, :] = start_probability * emission_probability[:, obs_seq[0]]

具体解释如下:

  • start_probability 对应于 π \pi π,即初始状态概率。
  • emission_probability[:, obs_seq[0]] 对应于 B ( o 0 ) B(o_0) B(o0),即在不同状态下观测到第一个观测值的概率。
  • obs_seq[0] 是观测序列的第一个观测值的索引。

通过逐元素相乘,得到在时刻 0 处于每个状态并观测到第一个观测值的概率,即 α 0 ( j ) \alpha_0(j) α0(j)

举例说明

假设有以下参数:

  • start_probability = np.array([0.6, 0.4])
  • emission_probability = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
  • obs_seq = [0, 1, 2](对应于 ‘walk’, ‘shop’, ‘clean’)

根据公式:

  • 对于状态 ‘Rainy’ ( j = 0 j = 0 j=0):
    α 0 ( 0 ) = π 0 × B 0 ( o 0 ) = 0.6 × 0.1 = 0.06 \alpha_0(0) = \pi_0 \times B_0(o_0) = 0.6 \times 0.1 = 0.06 α0(0)=π0×B0(o0)=0.6×0.1=0.06

  • 对于状态 ‘Sunny’ ( j = 1 j = 1 j=1):
    α 0 ( 1 ) = π 1 × B 1 ( o 0 ) = 0.4 × 0.6 = 0.24 \alpha_0(1) = \pi_1 \times B_1(o_0) = 0.4 \times 0.6 = 0.24 α0(1)=π1×B1(o0)=0.4×0.6=0.24

最终得到:

alpha[0, :] = [0.06, 0.24]

这表示在初始时刻观测到 ‘walk’ 时处于 ‘Rainy’ 状态的概率是 0.06,处于 ‘Sunny’ 状态的概率是 0.24。

递推计算
for t in range(1, len(obs_seq)):
    for j in range(len(states)):
        alpha[t, j] = np.dot(alpha[t-1, :], transition_probability[:, j]) * emission_probability[j, obs_seq[t]]
  • 外循环遍历观测序列的每一个时刻 t t t(从 1 开始,因为第 0 个时刻已经初始化)。

  • len(obs_seq) 是观测序列的长度。

  • 外循环遍历所有可能的状态 j j j

  • len(states) 是状态的总数。

  • alpha[t-1, :]:表示在时刻 t − 1 t-1 t1 各个状态的前向概率。

  • transition_probability[:, j]:表示从所有状态转移到状态 j j j 的概率。

  • np.dot(alpha[t-1, :], transition_probability[:, j]):计算在时刻 t t t 转移到状态 j j j 的总概率(前一时刻的前向概率和状态转移概率的点积)。

  • emission_probability[j, obs_seq[t]]:在状态 j j j 下观测到当前观测值 o t o_t ot 的概率。

  • 通过以上计算得到 alpha[t, j],即在时刻 t t t 处于状态 j j j 并观测到序列前 t t t 个观测值的概率。

  • 递推计算前向概率矩阵的后续行。对于每个时刻 t t t 和每个状态 j j j,前向概率的计算公式为:
    α [ t , j ] = ( ∑ i = 0 N − 1 α [ t − 1 , i ] × transition_probability [ i , j ] ) × emission_probability [ j , obs_seq [ t ] ] \alpha[t, j] = \left( \sum_{i=0}^{N-1} \alpha[t-1, i] \times \text{transition\_probability}[i, j] \right) \times \text{emission\_probability}[j, \text{obs\_seq}[t]] α[t,j]=(i=0N1α[t1,i]×transition_probability[i,j])×emission_probability[j,obs_seq[t]]
    这里使用了 NumPy 的 np.dot 函数计算前一行的前向概率和转移概率矩阵的点积。

前向算法递推计算的公式

对于每一个时刻 t t t 和每一个状态 j j j,前向概率 α [ t , j ] \alpha[t, j] α[t,j] 的计算公式为:

α [ t , j ] = ( ∑ i = 0 N − 1 α [ t − 1 , i ] × A i j ) × B j ( o t ) \alpha[t, j] = \left( \sum_{i=0}^{N-1} \alpha[t-1, i] \times A_{ij} \right) \times B_j(o_t) α[t,j]=(i=0N1α[t1,i]×Aij)×Bj(ot)

其中:

  • α [ t , j ] \alpha[t, j] α[t,j]:在时刻 t t t 处于状态 j j j 并观测到序列前 t t t 个观测值的概率。
  • α [ t − 1 , i ] \alpha[t-1, i] α[t1,i]:在时刻 t − 1 t-1 t1 处于状态 i i i 并观测到序列前 t − 1 t-1 t1 个观测值的概率。
  • A i j A_{ij} Aij:从状态 i i i 转移到状态 j j j 的概率。
  • B j ( o t ) B_j(o_t) Bj(ot):在状态 j j j 下观测到当前观测值 o t o_t ot 的概率。
  • N N N:状态的总数。

具体计算步骤举例

假设我们有如下参数:

  • start_probability = np.array([0.6, 0.4])
  • transition_probability = np.array([[0.7, 0.3], [0.4, 0.6]])
  • emission_probability = np.array([[0.1, 0.4, 0.5], [0.6, 0.3, 0.1]])
  • obs_seq = [0, 1, 2](对应于 ‘walk’, ‘shop’, ‘clean’)

初始化后:

α [ 0 , : ] = [ 0.06 , 0.24 ] \alpha[0, :] = [0.06, 0.24] α[0,:]=[0.06,0.24]

计算 t = 1 t = 1 t=1 时刻的前向概率:

  • 对于状态 ‘Rainy’ ( j = 0 j = 0 j=0):
    α [ 1 , 0 ] = ( α [ 0 , 0 ] × A 00 + α [ 0 , 1 ] × A 10 ) × B 0 ( o 1 ) \alpha[1, 0] = (\alpha[0, 0] \times A_{00} + \alpha[0, 1] \times A_{10}) \times B_0(o_1) α[1,0]=(α[0,0]×A00+α[0,1]×A10)×B0(o1)
    α [ 1 , 0 ] = ( 0.06 × 0.7 + 0.24 × 0.4 ) × 0.4 = 0.0336 \alpha[1, 0] = (0.06 \times 0.7 + 0.24 \times 0.4) \times 0.4 = 0.0336 α[1,0]=(0.06×0.7+0.24×0.4)×0.4=0.0336

  • 对于状态 ‘Sunny’ ( j = 1 j = 1 j=1):
    α [ 1 , 1 ] = ( α [ 0 , 0 ] × A 01 + α [ 0 , 1 ] × A 11 ) × B 1 ( o 1 ) \alpha[1, 1] = (\alpha[0, 0] \times A_{01} + \alpha[0, 1] \times A_{11}) \times B_1(o_1) α[1,1]=(α[0,0]×A01+α[0,1]×A11)×B1(o1)
    α [ 1 , 1 ] = ( 0.06 × 0.3 + 0.24 × 0.6 ) × 0.3 = 0.0468 \alpha[1, 1] = (0.06 \times 0.3 + 0.24 \times 0.6) \times 0.3 = 0.0468 α[1,1]=(0.06×0.3+0.24×0.6)×0.3=0.0468

以此类推,计算后续时刻 t = 2 t = 2 t=2 的前向概率。

计算总概率
total_prob = np.sum(alpha[-1, :])
  • 观测序列的总概率为前向概率矩阵最后一行的所有概率之和。
输出结果
print("Forward Probability Matrix:")
print(alpha)
print("\nTotal Probability of Observations:", total_prob)
  • 输出前向概率矩阵和观测序列的总概率。
代码执行结果示例

假设 alpha 矩阵计算结果如下:

[[0.06   0.24  ]
 [0.0336 0.0468]
 [0.0126 0.003 ]

total_prob0.0126 + 0.003 = 0.0156

综上所述,这段代码通过前向算法计算了给定观测序列在隐马尔可夫模型下的总概率,并输出了前向概率矩阵的中间结果。

2)解码问题

解码问题是指给定模型参数( π , A , B \pi, A, B π,A,B)和观察序列 O O O,找到最有可能的状态序列。常用维特比算法(Viterbi Algorithm)解决该问题。

维特比算法
δ t ( i ) = max ⁡ s 1 , s 2 , . . . , s t − 1 P ( s 1 , s 2 , . . . , s t = s i , o 1 , o 2 , . . . , o t ∣ λ ) \delta_t(i) = \max_{s_1, s_2, ..., s_{t-1}} P(s_1, s_2, ..., s_t = s_i, o_1, o_2, ..., o_t \mid \lambda) δt(i)=s1,s2,...,st1maxP(s1,s2,...,st=si,o1,o2,...,otλ)
其中, δ t ( i ) \delta_t(i) δt(i) 表示在时刻 t t t 到达状态 s i s_i si 且观察到 o 1 , o 2 , . . . , o t o_1, o_2, ..., o_t o1,o2,...,ot 的最大概率。

递推公式:
δ t + 1 ( j ) = max ⁡ i [ δ t ( i ) a i j ] b j ( o t + 1 ) \delta_{t+1}(j) = \max_i [\delta_t(i) a_{ij}] b_j(o_{t+1}) δt+1(j)=imax[δt(i)aij]bj(ot+1)

初始化:
δ 1 ( i ) = π i b i ( o 1 ) \delta_1(i) = \pi_i b_i(o_1) δ1(i)=πibi(o1)

回溯:
s t ∗ = arg ⁡ max ⁡ i [ δ T ( i ) ] s_t^* = \arg \max_i [\delta_T(i)] st=argimax[δT(i)]

代码示例
def viterbi_algorithm(A, B, pi, O):
    N = A.shape[0]  # 状态数
    T = len(O)      # 观察序列长度

    delta = np.zeros((T, N))
    psi = np.zeros((T, N), dtype=int)

    delta[0, :] = pi * B[:, O[0]]

    for t in range(1, T):
        for j in range(N):
            delta[t, j] = np.max(delta[t-1, :] * A[:, j]) * B[j, O[t]]
            psi[t, j] = np.argmax(delta[t-1, :] * A[:, j])

    states = np.zeros(T, dtype=int)
    states[T-1] = np.argmax(delta[T-1, :])
    for t in range(T-2, -1, -1):
        states[t] = psi[t+1, states[t+1]]

    return states

# 找到最有可能的状态序列
most_probable_states = viterbi_algorithm(A, B, pi, O)
print("最有可能的状态序列:", most_probable_states)
3)学习问题

学习问题是指给定观察序列 O O O,估计模型参数( π , A , B \pi, A, B π,A,B)。常用Baum-Welch算法(EM算法的一种)解决该问题。

Baum-Welch算法

  1. 初始化模型参数 π , A , B \pi, A, B π,A,B

  2. 计算前向概率 α \alpha α 和后向概率 β \beta β
    β t ( i ) = P ( o t + 1 , o t + 2 , . . . , o T ∣ s t = s i , λ ) \beta_t(i) = P(o_{t+1}, o_{t+2}, ..., o_T \mid s_t = s_i, \lambda) βt(i)=P(ot+1,ot+2,...,oTst=si,λ)

  3. 计算 ξ t ( i , j ) \xi_t(i, j) ξt(i,j) γ t ( i ) \gamma_t(i) γt(i)
    ξ t ( i , j ) = α t ( i ) a i j b j ( o t + 1 ) β t + 1 ( j ) P ( O ∣ λ ) \xi_t(i, j) = \frac{\alpha_t(i) a_{ij} b_j(o_{t+1}) \beta_{t+1}(j)}{P(O \mid \lambda)} ξt(i,j)=P(Oλ)αt(i)aijbj(ot+1)βt+1(j)
    γ t ( i ) = ∑ j = 1 N ξ t ( i , j ) \gamma_t(i) = \sum_{j=1}^{N} \xi_t(i, j) γt(i)=j=1Nξt(i,j)

  4. 更新参数 π , A , B \pi, A, B π,A,B
    π i = γ 1 ( i ) \pi_i = \gamma_1(i) πi=γ1(i)
    a i j = ∑ t = 1 T − 1 ξ t ( i , j ) ∑ t = 1 T − 1 γ t ( i ) a_{ij} = \frac{\sum_{t=1}^{T-1} \xi_t(i, j)}{\sum_{t=1}^{T-1} \gamma_t(i)} aij=t=1T1γt(i)t=1T1ξt(i,j)
    b j ( k ) = ∑ t = 1 , o t = k T γ t ( j ) ∑ t = 1 T γ t ( j ) b_j(k) = \frac{\sum_{t=1, o_t = k}^{T} \gamma_t(j)}{\sum_{t=1}^{T} \gamma_t(j)} bj(k)=t=1Tγt(j)t=1,ot=kTγt(j)

重复步骤2-4直到参数收敛。

代码示例
def baum_welch_algorithm(O, N, M, max_iter=100):
    T = len(O)
    A = np.random.rand(N, N)
    B = np.random.rand(N, M)
    pi = np.random.rand(N)

    A = A / A.sum(axis=1, keepdims=True)
    B = B / B.sum(axis=1, keepdims=True)
    pi = pi / pi.sum()

    for _ in range(max_iter):
        alpha = np.zeros((T, N))
        beta = np.zeros((T, N))



        alpha[0, :] = pi * B[:, O[0]]
        for t in range(1, T):
            for j in range(N):
                alpha[t, j] = np.sum(alpha[t-1, :] * A[:, j]) * B[j, O[t]]

        beta[T-1, :] = 1
        for t in range(T-2, -1, -1):
            for i in range(N):
                beta[t, i] = np.sum(A[i, :] * B[:, O[t+1]] * beta[t+1, :])

        xi = np.zeros((T-1, N, N))
        for t in range(T-1):
            denom = np.sum(alpha[t, :] * beta[t, :])
            for i in range(N):
                for j in range(N):
                    xi[t, i, j] = alpha[t, i] * A[i, j] * B[j, O[t+1]] * beta[t+1, j] / denom

        gamma = np.sum(xi, axis=2)

        pi = gamma[0, :]
        A = np.sum(xi, axis=0) / np.sum(gamma, axis=0, keepdims=True).T
        for j in range(N):
            for k in range(M):
                B[j, k] = np.sum(gamma[O == k, j]) / np.sum(gamma[:, j])

    return pi, A, B

# 示例数据:观察序列长度和观测空间大小
O = [0, 1, 0]
N = 2  # 状态数
M = 2  # 观察空间大小

# 估计模型参数
pi, A, B = baum_welch_algorithm(O, N, M)
print("估计的初始概率分布:", pi)
print("估计的状态转移概率矩阵:", A)
print("估计的观测概率矩阵:", B)

6.3 条件随机场

条件随机场(Conditional Random Field, CRF)是用于序列标注的概率图模型,与隐马尔可夫模型不同,CRF 是无向图模型,不假设观察序列的独立性。CRF 在自然语言处理等领域广泛应用。

公式推导

CRF 定义在给定观察序列 X = { x 1 , x 2 , . . . , x T } X = \{x_1, x_2, ..., x_T\} X={x1,x2,...,xT} 条件下,预测标签序列 Y = { y 1 , y 2 , . . . , y T } Y = \{y_1, y_2, ..., y_T\} Y={y1,y2,...,yT} 的条件概率分布。

CRF 的概率分布为:
P ( Y ∣ X ) = 1 Z ( X ) exp ⁡ ( ∑ t = 1 T ∑ k = 1 K λ k f k ( y t , y t − 1 , X , t ) ) P(Y \mid X) = \frac{1}{Z(X)} \exp \left( \sum_{t=1}^{T} \sum_{k=1}^{K} \lambda_k f_k(y_t, y_{t-1}, X, t) \right) P(YX)=Z(X)1exp(t=1Tk=1Kλkfk(yt,yt1,X,t))
其中, λ k \lambda_k λk 是模型参数, f k ( y t , y t − 1 , X , t ) f_k(y_t, y_{t-1}, X, t) fk(yt,yt1,X,t) 是特征函数, Z ( X ) Z(X) Z(X) 是归一化因子。

代码示例
import numpy as np
from sklearn_crfsuite import CRF

# 示例数据:观察序列和标签序列
X_train = [['I', 'love', 'Python'], ['CRF', 'is', 'great']]
y_train = [['O', 'O', 'B-Lang'], ['B-Tech', 'O', 'O']]

def word2features(sent, i):
    word = sent[i]
    features = {
        'bias': 1.0,
        'word.lower()': word.lower(),
        'word.isupper()': word.isupper(),
        'word.istitle()': word.istitle(),
        'word.isdigit()': word.isdigit(),
    }
    if i > 0:
        word1 = sent[i-1]
        features.update({
            '-1:word.lower()': word1.lower(),
            '-1:word.istitle()': word1.istitle(),
            '-1:word.isupper()': word1.isupper(),
        })
    else:
        features['BOS'] = True

    if i < len(sent)-1:
        word1 = sent[i+1]
        features.update({
            '+1:word.lower()': word1.lower(),
            '+1:word.istitle()': word1.istitle(),
            '+1:word.isupper()': word1.isupper(),
        })
    else:
        features['EOS'] = True

    return features

def sent2features(sent):
    return [word2features(sent, i) for i in range(len(sent))]

X_train_feats = [sent2features(s) for s in X_train]

# 创建并训练CRF模型
crf = CRF(algorithm='lbfgs', max_iterations=100)
crf.fit(X_train_feats, y_train)

# 预测
X_test = [['I', 'enjoy', 'CRF']]
X_test_feats = [sent2features(s) for s in X_test]
y_pred = crf.predict(X_test_feats)
print("预测结果:", y_pred)

7. 一些疑问

7.1 马尔可夫中的协方差有什么作用

在马尔可夫模型(尤其是隐马尔可夫模型,Hidden Markov Model, HMM)和其他相关模型中,协方差主要体现在以下几个方面:

1. 高斯隐马尔可夫模型(Gaussian Hidden Markov Model, GHMM)

在高斯隐马尔可夫模型中,状态的观测值被假定为服从高斯分布(正态分布)。每个状态 s i s_i si 对应一个高斯分布,其参数包括均值向量 μ i \mu_i μi 和协方差矩阵 Σ i \Sigma_i Σi。协方差矩阵描述了观测值在各个维度上的变异性及其相互关系。

观测概率

对于状态 i i i 下的观测值 o \mathbf{o} o,其观测概率 B i ( o ) B_i(\mathbf{o}) Bi(o) 表示为:

B i ( o ) = 1 ( 2 π ) d / 2 ∣ Σ i ∣ 1 / 2 exp ⁡ ( − 1 2 ( o − μ i ) ⊤ Σ i − 1 ( o − μ i ) ) B_i(\mathbf{o}) = \frac{1}{(2\pi)^{d/2} |\Sigma_i|^{1/2}} \exp \left( -\frac{1}{2} (\mathbf{o} - \mu_i)^\top \Sigma_i^{-1} (\mathbf{o} - \mu_i) \right) Bi(o)=(2π)d/2Σi1/21exp(21(oμi)Σi1(oμi))

其中:

  • d d d 是观测值的维度。
  • ∣ Σ i ∣ |\Sigma_i| Σi 是协方差矩阵的行列式。
  • Σ i − 1 \Sigma_i^{-1} Σi1 是协方差矩阵的逆。
2. 协方差矩阵的作用
描述观测数据的分布形状

协方差矩阵 Σ i \Sigma_i Σi 描述了在状态 i i i 下观测数据的分布形状和方向。协方差矩阵中的每个元素 σ j k \sigma_{jk} σjk 表示第 j j j 和第 k k k 维度之间的协方差。如果协方差矩阵是对角矩阵,则表示各个维度是相互独立的;如果不是对角矩阵,则表示各个维度之间存在相关性。

调整模型的灵活性

在高斯隐马尔可夫模型中,协方差矩阵决定了高斯分布的形状和方向。因此,协方差矩阵可以调整模型的灵活性,使其能够更好地适应具有不同相关性和变异性的观测数据。

计算观测概率

协方差矩阵在计算观测概率时起着至关重要的作用。通过协方差矩阵,可以确定观测值与高斯分布中心的距离,从而计算观测值的概率。

3. 实现示例

下面是一个简单的示例代码,演示如何在高斯隐马尔可夫模型中使用协方差矩阵计算观测概率。

import numpy as np
from scipy.stats import multivariate_normal

# 定义高斯分布的参数
mean = np.array([0, 0])  # 均值向量
cov = np.array([[1, 0.5], [0.5, 1]])  # 协方差矩阵

# 创建多变量正态分布
rv = multivariate_normal(mean, cov)

# 计算观测值的概率
observation = np.array([1, 1])
prob = rv.pdf(observation)

print("Observation Probability:", prob)

在这个示例中,我们定义了一个二维高斯分布,其均值向量为 [ 0 , 0 ] [0, 0] [0,0],协方差矩阵为 [ 1 0.5 0.5 1 ] \begin{bmatrix} 1 & 0.5 \\ 0.5 & 1 \end{bmatrix} [10.50.51]。然后,我们使用 scipy.stats 中的 multivariate_normal 函数来创建这个高斯分布,并计算观测值 [ 1 , 1 ] [1, 1] [1,1] 的概率。

在马尔可夫模型(特别是高斯隐马尔可夫模型)中,协方差矩阵主要用于描述观测数据的分布形状和方向。它不仅影响观测概率的计算,还决定了模型的灵活性和适应性。因此,在实际应用中,协方差矩阵是一个重要的参数,对模型性能有着重要影响。

8. 学习心得

本次学习了时间序列的基本概念,了解到时间序列可以通过分解模型来描述(加法模型、乘法模型和组合模型),分解模型将时间序列分解为不同的成分(趋势成分、季节成分、循环成分、不规则成分)。使用statsmodels.tsa.seasonal对零售销售数据进行组合模型分解,并画图展示时间序列数据分解成不通的成分。
移动平均法是一种时间序列分析的方法,通过取固定时间窗口内数据的平均值来平滑时间序列。它适用于平滑短期波动,突出长期趋势或周期。
指数平均法(Exponential Smoothing)是一种加权移动平均方法,越近期的数据权重越大。它能够较好地应对时间序列中的趋势和季节性变化。
在介绍自回归模型 (AR) 和移动平均模型 (MA) 后,引入了前二者的结合ARMA模型,用于描述平稳时间序列的数据。
ARIMA模型扩展了ARMA模型,适用于非平稳时间序列。ARIMA模型通过差分操作将非平稳时间序列转化为平稳序列,然后再应用ARMA模型。
组合投资策略是通过分散投资来降低风险并优化回报的金融方法。本次学习了三种常见的组合投资策略:马科维兹均值方差模型、最大夏普比率模型和风险平价模型。
马尔可夫模型基于马尔可夫性质,即未来状态只依赖于当前状态,而与过去状态无关。马尔可夫模型的实现可以分为三类主要问题:评估问题、解码问题和学习问题。通过人们对不通状态的天气(晴天、雨天)下通过前向算法递推计算人们的不同反应(行走、购物、清理)
(ps:学习本章的时候不太理解为什么会在本章节介绍马尔可夫模型???)

  • 17
    点赞
  • 19
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值