# Python实战：用ARIMA模型预测餐厅销量的5个关键步骤（附完整代码） 最近和一位经营连锁餐厅的朋友聊天，他正为库存管理和人员排班头疼。食材备多了怕浪费，备少了又影响顾客体验；员工排班全凭经验，忙时人手不够，闲时又人力过剩。他问我，有没有什么数据驱动的办法，能提前知道未来几天大概会有多少客人？这让我想起了时间序列预测这个经典工具。对于餐厅、零售这类有明显周期性和趋势性的业务，ARIMA模型就像一位经验丰富的“老掌柜”，能从过去的数据里嗅出未来的味道。今天，我们就抛开复杂的数学公式，直接上手Python，用五个清晰的关键步骤，手把手带你构建一个能预测餐厅销量的ARIMA模型。无论你是想优化自家小店运营的数据分析师，还是希望将预测能力融入产品的开发者，这套流程都能让你快速落地，看到实实在在的预测结果。 ## 1. 理解你的数据：时间序列预测的起点 在动手写代码之前，我们得先和手里的数据“交个朋友”。对于餐厅销量数据，它通常是一个按时间顺序排列的序列，比如过去365天每天的营业额。时间序列分析的核心假设是“未来会延续过去的模式”，但这个模式必须满足一定的条件，模型才能捕捉到。 首先，我们需要明确两个核心概念：**平稳性**和**白噪声**。你可以把平稳的时间序列想象成一条在固定轨道上平稳运行的过山车，它围绕一个固定的均值上下波动，且波动的幅度（方差）大致恒定。而非平稳序列则像一条没有固定轨道的过山车，可能持续上升（趋势）或波动越来越大。ARIMA模型要求输入的数据是平稳的，或者能通过差分变得平稳。 > 提示：为什么强调平稳性？因为只有平稳序列的统计特性（如均值、方差）不随时间变化，模型基于历史数据学习到的规律才能稳定地应用于未来预测。否则，模型会“迷失”在不断变化的趋势里。 那么，餐厅销量数据天然就是平稳的吗？大概率不是。它通常包含： - **趋势性**：比如随着店铺口碑传播，整体销量呈缓慢上升趋势。 - **季节性**：每周的周末销量更高，每年节假日会有销售高峰。 - **周期性**：可能还存在以月或季度为单位的更长周期波动。 - **随机性**：由天气、临时促销、突发事件等不可控因素引起的波动。 我们的第一步，就是通过可视化和统计检验，来诊断数据的这些特性。下面是一个快速查看数据基本情况的代码示例： ```python import pandas as pd import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.stattools import adfuller # 假设数据已加载为DataFrame，包含‘date’和‘sales’两列 # df = pd.read_csv('restaurant_sales.csv', parse_dates=['date'], index_col='date') print("数据概览:") print(df.info()) print("\n前5行数据:") print(df.head()) # 绘制原始序列时序图 plt.figure(figsize=(14, 6)) plt.subplot(2, 1, 1) plt.plot(df.index, df['sales'], label='原始销量') plt.title('餐厅原始销量时序图') plt.xlabel('日期') plt.ylabel('销量') plt.legend() plt.grid(True) # 绘制自相关图(ACF) plt.subplot(2, 1, 2) plot_acf(df['sales'], lags=40, ax=plt.gca()) # 查看40期滞后 plt.title('原始序列自相关图(ACF)') plt.tight_layout() plt.show() ``` 通过时序图，你能直观看到数据是否存在明显的上升/下降趋势或周期性波动。自相关图则更量化：如果自相关系数随着滞后阶数增加缓慢衰减（不快速趋近于0），则强烈暗示序列非平稳。 ## 2. 让数据“安静”下来：平稳性检验与差分处理 看完图表有了感性认识，接下来需要用统计检验来定量判断平稳性。最常用的方法是**Augmented Dickey-Fuller (ADF) 单位根检验**。它的原假设是“序列存在单位根，即非平稳”。如果检验得到的p值小于显著性水平（通常为0.05），我们就拒绝原假设，认为序列是平稳的。 对于非平稳序列，ARIMA模型中的 **“I” (Integrated)** 部分就派上用场了，即通过差分运算来消除趋势或季节性。一阶差分就是计算当前值与前一个值的差值（`df['sales'].diff(1)`），这通常能消除线性趋势。如果一阶差分后仍不平稳，可能需要进行二阶差分，或考虑季节性差分。 ```python # 进行ADF平稳性检验 def adf_test(timeseries): print('ADF检验结果:') result = adfuller(timeseries.dropna(), autolag='AIC') # 自动选择最佳滞后阶数 output = pd.Series(result[0:4], index=['ADF统计量', 'p-value', '滞后阶数', '观测数']) for key, value in result[4].items(): output[f'临界值({key})'] = value print(output.to_string()) if result[1] <= 0.05: print("-> p值 <= 0.05，拒绝原假设，序列平稳。") else: print("-> p值 > 0.05，无法拒绝原假设，序列非平稳。") return result[1] p_value_original = adf_test(df['sales']) # 如果非平稳，进行一阶差分 if p_value_original > 0.05: df['sales_diff_1'] = df['sales'].diff(1) print("\n对原始序列进行一阶差分后...") p_value_diff = adf_test(df['sales_diff_1']) # 可视化差分后的序列 plt.figure(figsize=(14, 5)) plt.plot(df.index[1:], df['sales_diff_1'].dropna(), color='green') plt.title('一阶差分后序列时序图') plt.xlabel('日期') plt.ylabel('销量差分值') plt.grid(True) plt.show() ``` 除了ADF检验，**白噪声检验**也至关重要。白噪声序列是完全随机的，前后数据点毫无关联，对这种序列做预测没有意义。我们建模的前提是，差分后的平稳序列必须是**非白噪声**序列，即存在可被模型捕捉的相关性。可以使用Ljung-Box检验来判断。 ```python from statsmodels.stats.diagnostic import acorr_ljungbox # 对（差分后的）平稳序列进行白噪声检验 lb_result = acorr_ljungbox(平稳序列, lags=[10], return_df=True) # 检验前10阶 print("Ljung-Box白噪声检验结果 (滞后10阶):") print(lb_result) if lb_result['lb_pvalue'].iloc[0] < 0.05: print("-> p值 < 0.05，拒绝原假设，序列为非白噪声，可以进一步建模。") else: print("-> 序列可能为白噪声，需重新审视数据或问题。") ``` ## 3. 为模型“配钥匙”：确定ARIMA的(p,d,q)参数 经过第二步，我们得到了一个平稳、非白噪声的序列，并确定了差分阶数 **d**（例如，一阶差分则d=1）。接下来是最关键也最具技巧性的一步：确定自回归阶数 **p** 和移动平均阶数 **q**。这就像为模型这把锁找到最合适的钥匙。 这里主要依靠两把“利器”：**自相关函数图** 和 **偏自相关函数图**。 - **ACF图**：描述当前观测值与过去各期观测值之间的总体相关性。它拖尾（逐渐衰减至0）或截尾（在某一阶后突然接近0）的模式，能提示q值。 - **PACF图**：在控制了中间滞后项的影响后，当前观测值与某一特定滞后观测值之间的“纯”相关性。它的截尾模式能提示p值。 | 模型识别准则 | ACF表现 | PACF表现 | 建议的(p, q) | | :--- | :--- | :--- | :--- | | **AR(p)模型** | 拖尾（指数衰减或正弦波动） | **p阶后截尾** | p = 截尾处阶数， q=0 | | **MA(q)模型** | **q阶后截尾** | 拖尾 | p=0, q = 截尾处阶数 | | **ARMA(p,q)模型** | 拖尾 | 拖尾 | p, q均需从图中初步判断，再结合定阶准则 | ```python # 绘制平稳序列的ACF和PACF图，辅助定阶 fig, axes = plt.subplots(1, 2, figsize=(14, 4)) # ACF图 plot_acf(平稳序列, lags=40, ax=axes[0]) axes[0].set_title('平稳序列自相关图(ACF)') # PACF图 plot_pacf(平稳序列, lags=40, ax=axes[1], method='ywm') # 使用Yule-Walker方程估计 axes[1].set_title('平稳序列偏自相关图(PACF)') plt.tight_layout() plt.show() ``` 看图是门艺术，存在主观性。因此，我们通常结合**信息准则**进行定量定阶，最常用的是 **AIC** 和 **BIC**。其原理是，在保证模型拟合优度的前提下，选择使AIC或BIC值最小的(p, q)组合，以平衡模型的复杂度和拟合能力。我们可以通过网格搜索来实现： ```python import itertools import warnings from statsmodels.tsa.arima.model import ARIMA warnings.filterwarnings("ignore") # 忽略模型拟合中的警告 # 确定p, q的搜索范围（根据ACF/PACF截尾处初步判断） p_range = range(0, 4) # 例如，PACF在2阶后截尾，可设0-3 q_range = range(0, 4) # 例如，ACF在1阶后截尾，可设0-3 d = 1 # 之前确定的一阶差分 best_aic = float('inf') best_bic = float('inf') best_order_aic = None best_order_bic = None results_list = [] for p, q in itertools.product(p_range, q_range): try: model = ARIMA(df['sales'], order=(p, d, q)) model_fit = model.fit() aic = model_fit.aic bic = model_fit.bic results_list.append([p, q, aic, bic]) if aic < best_aic: best_aic = aic best_order_aic = (p, d, q) if bic < best_bic: best_bic = bic best_order_bic = (p, d, q) except Exception as e: continue # 跳过无法拟合的参数组合 # 将结果转为DataFrame方便查看 results_df = pd.DataFrame(results_list, columns=['p', 'q', 'AIC', 'BIC']) print("不同(p,q)组合的AIC/BIC值:") print(results_df.sort_values('AIC').head()) print(f"\n根据AIC准则，最优参数为: {best_order_aic}") print(f"根据BIC准则（惩罚模型复杂度更强），最优参数为: {best_order_bic}") ``` 通常，AIC和BIC给出的最优阶数会接近。如果差异较大，可以优先选择更简洁的模型（即p、q更小的），或者用样本外预测效果来最终裁决。 ## 4. 构建模型与诊断：不仅仅是拟合，更要验证 确定了(p,d,q)参数，我们就可以正式构建ARIMA模型并拟合数据了。但拟合完成并不意味着大功告成，我们必须对模型进行严格的诊断，核心是检验**残差序列**。一个理想的模型，其残差应该近似为一个**均值为0、方差恒定、且无自相关性的白噪声序列**。如果残差还有模式可循，说明模型没有完全捕捉数据中的信息，有改进空间。 诊断主要从三个方面进行： 1. **残差序列的ACF/PACF图**：检查残差在各阶滞后上是否还存在显著的自相关。理想情况是所有滞后阶的自相关系数都落在置信区间内（图中蓝色阴影区域）。 2. **残差的正态性检验**：通过QQ图或统计检验（如Shapiro-Wilk检验）查看残差是否近似服从正态分布。这对预测区间的构建很重要。 3. **Ljung-Box检验（再次使用）**：对残差序列进行白噪声检验，确认其是否已无统计上显著的自相关。 ```python # 使用选定的最优参数拟合模型（这里以AIC最优为例） final_order = best_order_aic model = ARIMA(df['sales'], order=final_order) model_fitted = model.fit() print("="*50) print("模型拟合报告:") print(model_fitted.summary()) print("="*50) # 1. 获取残差 residuals = model_fitted.resid[模型拟合的起始索引:] # 注意排除前几期因差分和滞后缺失的数据 # 2. 绘制残差诊断图 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 2.1 残差时序图 axes[0, 0].plot(residuals.index, residuals) axes[0, 0].axhline(y=0, color='r', linestyle='--', alpha=0.5) axes[0, 0].set_title('残差序列时序图') axes[0, 0].set_xlabel('日期') axes[0, 0].set_ylabel('残差') # 2.2 残差ACF图 plot_acf(residuals, lags=40, ax=axes[0, 1]) axes[0, 1].set_title('残差自相关图(ACF)') # 2.3 残差直方图与正态分布曲线 from scipy.stats import norm import numpy as np axes[1, 0].hist(residuals, bins=30, density=True, alpha=0.7, edgecolor='black') mu, std = norm.fit(residuals) xmin, xmax = axes[1, 0].get_xlim() x = np.linspace(xmin, xmax, 100) p = norm.pdf(x, mu, std) axes[1, 0].plot(x, p, 'k', linewidth=2, label=f'拟合正态分布\n(μ={mu:.2f}, σ={std:.2f})') axes[1, 0].set_title('残差分布直方图') axes[1, 0].legend() # 2.4 残差QQ图 from statsmodels.graphics.gofplots import qqplot qqplot(residuals, line='s', ax=axes[1, 1]) # 's'表示标准线 axes[1, 1].set_title('残差QQ图（检验正态性）') plt.tight_layout() plt.show() # 3. 对残差进行Ljung-Box白噪声检验 lb_test_result = acorr_ljungbox(residuals, lags=[10, 20], return_df=True) print("\n残差序列Ljung-Box检验结果:") print(lb_test_result) if (lb_test_result['lb_pvalue'] > 0.05).all(): print("-> 所有滞后阶的p值均大于0.05，无法拒绝残差为白噪声的原假设，模型诊断通过。") else: print("-> 部分滞后阶p值小于0.05，残差可能仍存在自相关，模型有待改进。") ``` 如果诊断发现残差不是白噪声，可能需要回到第三步，尝试调整(p, q)参数，或者考虑更复杂的模型，如**季节性ARIMA**来捕捉周度、月度等固定周期模式。 ## 5. 预测未来与结果解读：将数字转化为洞察 模型通过诊断后，我们就可以放心地用它进行预测了。`statsmodels`的`get_forecast()`或`predict()`方法可以方便地生成未来多期的点预测和置信区间。点预测是模型认为最可能发生的值，而置信区间（通常是95%）给出了预测值可能波动的范围，这在实际业务中往往比单一预测值更有参考价值。 ```python # 进行未来7天的预测 forecast_steps = 7 forecast_result = model_fitted.get_forecast(steps=forecast_steps) # 获取点预测值、预测区间上下限和标准误 forecast_mean = forecast_result.predicted_mean forecast_ci = forecast_result.conf_int(alpha=0.05) # 95%置信区间 # 创建包含历史数据和预测结果的DataFrame forecast_index = pd.date_range(start=df.index[-1] + pd.Timedelta(days=1), periods=forecast_steps, freq='D') forecast_df = pd.DataFrame({ '预测销量': forecast_mean.values, '预测下限_95%': forecast_ci.iloc[:, 0].values, '预测上限_95%': forecast_ci.iloc[:, 1].values }, index=forecast_index) print("未来7天销量预测:") print(forecast_df) # 可视化历史数据、拟合值及未来预测 plt.figure(figsize=(14, 7)) # 绘制历史数据 plt.plot(df.index, df['sales'], label='历史实际销量', color='blue', alpha=0.7) # 绘制样本内拟合值（可选，查看模型拟合效果） fitted_values = model_fitted.fittedvalues plt.plot(fitted_values.index, fitted_values, label='模型拟合值', color='green', alpha=0.9, linestyle='--') # 绘制未来预测 plt.plot(forecast_df.index, forecast_df['预测销量'], label='未来预测销量', color='red', marker='o') plt.fill_between(forecast_df.index, forecast_df['预测下限_95%'], forecast_df['预测上限_95%'], color='red', alpha=0.2, label='95% 置信区间') plt.title('餐厅销量：历史、拟合与未来预测') plt.xlabel('日期') plt.ylabel('销量') plt.legend(loc='best') plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() ``` 拿到预测结果后，如何向业务方解读？这里有几个要点： - **点预测值**：可以告诉店长，“模型预测下周一销量大约在XXX份”。这是安排备货和基础人力的核心参考。 - **置信区间**：必须强调预测的不确定性。“虽然最可能值是XXX，但有95%的把握认为销量会在[下限, 上限]这个范围内。”这为制定弹性预案（如准备安全库存）提供了依据。 - **结合业务常识**：如果模型预测周末销量反而下降，这明显违背常识，就需要回头检查数据或模型是否遗漏了重要的季节性因素（比如未使用SARIMA模型）。永远要让数据科学为业务逻辑服务，而不是替代它。 最后，别忘了模型需要**持续更新**。餐厅的销售模式会变（比如换了菜单、附近开了新店），用太久远的历史数据训练模型，预测未来可能会失灵。一个实用的做法是，定期（比如每季度）用最新的数据重新训练模型，或者建立一套自动化监控流程，当预测误差持续增大时触发模型重训。在实际项目中，我习惯将整个流程脚本化，从数据拉取、预处理、模型训练到预测输出，形成一个端到端的Pipeline，这样每次更新只需运行一个脚本，大大提升了落地和维护的效率。