在 python/scikit/numpy 中替代 r 的指数平滑状态空间模型

Posted 2023-03-12

【中文标题】在 python/scikit/numpy 中替代 r 的指数平滑状态空间模型

【英文标题】：Alternative for r's Exponential smoothing state space model in python/scikit/numpy

【发布时间】：2016-05-13 13:01:05

【问题描述】：

在 R 中，我们有一个很好的预测模型，例如：

ets(y, model="ZZZ", damped=NULL, alpha=NULL, beta=NULL, gamma=NULL, 

phi=NULL, additive.only=FALSE, lambda=NULL, 

lower=c(rep(0.0001,3), 0.8), upper=c(rep(0.9999,3),0.98), 

opt.crit=c("lik","amse","mse","sigma","mae"), nmse=3, 

bounds=c("both","usual","admissible"), ic=c("aicc","aic","bic"),

restrict=TRUE, allow.multiplicative.trend=FALSE, use.initial.values=FALSE, ...)

在这种方法中，如果我们分配任何变量，它会自动获取季节类型、趋势和错误类型，例如model="ZZZ"/"AMA"/"MMZ"，并且会自动调整一些因素以获得准确的结果。

在 python 中，我们是否有任何类似于 ets 的内容 熊猫/numpy/scipy/scikit？

根据我的研究：pandas 中的Ewma 类似，但我们需要将所有参数硬编码为固定参数。 在 Holtwinter 中，我们需要为所有趋势和季节类型编写详细的方法。

所以我们有任何现成的函数来代替它 数据框作为输入并提供预测值，无需编写 我们自己有任何参数的内部函数吗？

任何微调回归模型 scikit/statsmodels？

【问题讨论】：

这里有一些关于 R 包采用的方法的参考：otexts.org/fpp/7/7robjhyndman.com/talks/ABS1.pdf。到目前为止我还没有找到实现完整状态空间框架的python包。

您可以通过rpy2 包使用Python 中的R 工具吗？

@mfripp 是的，我可以从 python 调用 R，但如果可以的话，我更喜欢直接使用 python！

statsmodels 状态空间模型看起来很有前途：statsmodels.org/dev/statespace.html 特别是自定义模型功能：statsmodels.org/dev/statespace.html#custom-state-space-models

这很有希望，但也不完整：github.com/mcskinner/ets

【参考方案1】：

经过一番搜索后，我没有发现任何似乎真正有希望作为 python 的 ets 替代品的东西。不过也有一些尝试：StatsModels 和 pycast's Forecasting methods，您可以检查它们是否适合您的需求。

您可以用来解决缺少的实现的一个选项是使用subprocess 模块从python 运行R 脚本。有一篇很好的文章告诉你如何做到这一点here。

为了以后做：

您需要创建一个 R 脚本（例如my_forecast.R），它将 计算（使用ets）并打印文件或stdout（使用cat()命令）上的预测，以便在脚本 运行。

您可以按如下方式从 python 脚本运行 R 脚本：

import subprocess

# You need to define the command that will run the Rscript from the subprocess
command = 'Rscript'
path2script = 'path/to/my_forecast.R'
cmd = [command, path2script]

# Option 1: If your script prints to a file
subprocess.run(cmd)
f = open('path/to/created/file', 'r')
(...Do stuff from here...)

# Option 2: If your script prints to stdout
forecasts = subprocess.check_output(cmd, universal_newlines=True)
(...Do stuff from here...)

您还可以向cmd 添加参数，Rscript 将使用这些参数作为命令行参数，如下所示：

args = [arg0, arg1, ...]

cmd = [command, path2script] + args 
Then pass cmd to the subprocess

---

编辑：

我找到了有关 Holt-Winters 预测的一系列示例文章：part1、part2 和 part3。除了这些文章中易于理解的分析之外，Gregory Trubetskoy（作者）还提供了他开发的代码：

初步趋势：

def initial_trend(series, slen):
    sum = 0.0
    for i in range(slen):
        sum += float(series[i+slen] - series[i]) / slen
    return sum / slen

# >>> initial_trend(series, 12)
# -0.7847222222222222

初始季节性组件：

def initial_seasonal_components(series, slen):
    seasonals = 
    season_averages = []
    n_seasons = int(len(series)/slen)
    # compute season averages
    for j in range(n_seasons):
        season_averages.append(sum(series[slen*j:slen*j+slen])/float(slen))
    # compute initial values
    for i in range(slen):
        sum_of_vals_over_avg = 0.0
        for j in range(n_seasons):
            sum_of_vals_over_avg += series[slen*j+i]-season_averages[j]
        seasonals[i] = sum_of_vals_over_avg/n_seasons
    return seasonals

# >>> initial_seasonal_components(series, 12)
# 0: -7.4305555555555545, 1: -15.097222222222221, 2: -7.263888888888888,
#  3: -5.097222222222222,  4: 3.402777777777778,   5: 8.069444444444445,  
#  6: 16.569444444444446,  7: 9.736111111111112,   8: -0.7638888888888887,
#  9: 1.902777777777778,  10: -3.263888888888889, 11: -0.7638888888888887

最后是算法：

def triple_exponential_smoothing(series, slen, alpha, beta, gamma, n_preds):
    result = []
    seasonals = initial_seasonal_components(series, slen)
    for i in range(len(series)+n_preds):
        if i == 0: # initial values
            smooth = series[0]
            trend = initial_trend(series, slen)
            result.append(series[0])
            continue
        if i >= len(series): # we are forecasting
            m = i - len(series) + 1
            result.append((smooth + m*trend) + seasonals[i%slen])
        else:
            val = series[i]
            last_smooth, smooth = smooth, alpha*(val-seasonals[i%slen]) + (1-alpha)*(smooth+trend)
            trend = beta * (smooth-last_smooth) + (1-beta)*trend
            seasonals[i%slen] = gamma*(val-smooth) + (1-gamma)*seasonals[i%slen]
            result.append(smooth+trend+seasonals[i%slen])
    return result

# # forecast 24 points (i.e. two seasons)
# >>> triple_exponential_smoothing(series, 12, 0.716, 0.029, 0.993, 24)
# [30, 20.34449316666667, 28.410051892109554, 30.438122252647577, 39.466817731253066, ...

您可以将它们放在一个文件中，例如：holtwinters.py 在具有以下结构的文件夹中：

forecast_folder
|
└── __init__.py
|
└── holtwinters.py

从这里开始，这是一个 python 模块，您可以将它放置在您想要的每个项目结构中，并在该项目中的任何位置使用它，只需导入它。

【讨论】：

我真的很想避免调用 R 子进程。我很欣赏详细的答案，但我设置了赏金，因为我希望在 python 中找到一些东西。

目前，我在答案的开头已经链接了一些尝试，@jseabold 是 StatsModels 的初始开发人员，所以如果你询问他，他可能会对此事有所了解!

我已经添加了一种现成的方法来在 python 中进行 Holt Winters 预测

我给你赏金，因为你付出了很多努力来回答这个问题，但这仍然不是我真正想要的。如果将来有人找到（或创建）了 Rob Hyndmans 的所有 30 个状态空间指数平滑模型的实现，我将很高兴再次提供奖励：otexts.org/fpp/7/7

@Zach 是我更需要的答案吗？

【参考方案2】：

ETS 现在作为 Statsmodels 最新的release 的一部分在 python 中可用。下面是一个简单的使用示例：

from statsmodels.tsa.api import ExponentialSmoothing

# make our own periodic data
# 3 weeks of hourly data
m = 24
days = 21
x = np.array(range(days * m)) / (m / 2 / np.pi)
st = 4 * np.sin(x)
lt = np.linspace(0, 15, days * m)
bt = 100
e = np.random.normal(scale=0.5, size=x.shape)

# make the ts we wish to forecast
y = lt + st + bt + e

# our guessed parameters
alpha = 0.4
beta = 0.2
gamma = 0.01

# initialise model
ets_model = ExponentialSmoothing(y, trend='add', seasonal='add', 
seasonal_periods=24)
ets_fit = ets_model.fit(smoothing_level=alpha, smoothing_slope=beta,
smoothing_seasonal=gamma)

# forecast p hours ahead
p_ahead = 48
yh = ets_fit.forecast(p_ahead)

# plot the y, y_smoothed and y_hat ts'
plt.plot(y, label='y')
plt.plot(ets_fit.fittedvalues, label='y_smooth')
plt.plot(range(days * m, days * m + p_ahead),yh, label='y_hat')

plt.legend()
plt.show()

这里有更多 examples 以笔记本形式出现。

最后，这里是source code，如果你想看看的话。