AI数据分析实战：构建端到端的预测系统——从数据清洗到模型部署

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

# 生成模拟数据：100家门店，10个品类，过去2年每日数据
np.random.seed(42)
dates = pd.date_range('2024-01-01', '2025-12-31', freq='D')
stores = range(1, 101)
families = ['FAMILY_' + str(i) for i in range(1, 11)]

index = pd.MultiIndex.from_product([dates, stores, families], names=['date', 'store', 'family'])
df = pd.DataFrame(index=index).reset_index()

# 基础销量（带趋势、季节性、随机噪声）
df['sales'] = 100 + np.sin(2*np.pi*df['date'].dt.dayofyear/365)*30 \
              + df['date'].dt.dayofyear/365 * 20 \
              + np.random.normal(0, 10, len(df))

# 促销标志（随机10%的天数有促销）
df['onpromotion'] = np.random.choice([0,1], size=len(df), p=[0.9,0.1])

# 节假日标志（随机20个节假日）
holidays = np.random.choice(dates, size=20, replace=False)
df['holiday'] = df['date'].isin(holidays).astype(int)

# 温度（随机）
df['temp'] = 20 + 5*np.sin(2*np.pi*df['date'].dt.dayofyear/365) + np.random.normal(0,2,len(df))

# 刻意引入缺失值
df.loc[np.random.choice(df.index, size=5000, replace=False), 'sales'] = np.nan
df.loc[np.random.choice(df.index, size=2000, replace=False), 'temp'] = np.nan

print(df.shape)  # (730000, 6)
df.head()

from scipy import interpolate

# 按门店+品类分组，进行线性插值（考虑时间顺序）
df = df.sort_values(['store', 'family', 'date'])
df['sales'] = df.groupby(['store', 'family'])['sales'].transform(
    lambda x: x.interpolate(method='linear', limit_direction='both')
)

# 温度填充：按门店+月份中位数
df['temp'] = df.groupby(['store', df['date'].dt.month])['temp'].transform(
    lambda x: x.fillna(x.median())
)

# 异常值检测（按门店+品类+月份，滚动窗口）
def detect_outliers(group, window=30, n_sigmas=3):
    rolling_median = group.rolling(window, center=True).median()
    rolling_mad = group.rolling(window, center=True).apply(lambda x: np.median(np.abs(x - np.median(x))))
    lower = rolling_median - n_sigmas * rolling_mad
    upper = rolling_median + n_sigmas * rolling_mad
    return (group < lower) | (group > upper)

df['outlier'] = df.groupby(['store', 'family', df['date'].dt.month])['sales'].transform(
    lambda x: detect_outliers(x, window=30, n_sigmas=3)
)
# 将异常值替换为滚动中位数
df.loc[df['outlier'], 'sales'] = df.groupby(['store', 'family'])['sales'].transform(
    lambda x: x.rolling(30, center=True).median()
)

import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.seasonal import seasonal_decompose

daily_sales = df.groupby('date')['sales'].sum()
decomp = seasonal_decompose(daily_sales, model='additive', period=7)
decomp.plot()
plt.show()

# 按门店+品类分组，生成滞后和滚动特征
features = df.copy()
features = features.sort_values(['store', 'family', 'date'])

for lag in [1, 7, 14, 28]:
    features[f'sales_lag_{lag}'] = features.groupby(['store', 'family'])['sales'].shift(lag)

for window in [7, 30]:
    for stat in ['mean', 'std', 'min', 'max']:
        features[f'sales_roll_{window}_{stat}'] = features.groupby(['store', 'family'])['sales'].transform(
            lambda x: x.rolling(window, min_periods=1).agg(stat).shift(1)  # shift(1)避免包含当前
        )

# 时间特征
features['dayofweek'] = features['date'].dt.dayofweek
features['month'] = features['date'].dt.month
features['dayofyear'] = features['date'].dt.dayofyear
features['weekofyear'] = features['date'].dt.isocalendar().week
features['quarter'] = features['date'].dt.quarter
features['is_weekend'] = (features['dayofweek'] >= 5).astype(int)

# 促销相关：过去7天促销次数
features['promo_count_7d'] = features.groupby(['store', 'family'])['onpromotion'].transform(
    lambda x: x.rolling(7, min_periods=1).sum().shift(1)
)

# 温度滞后
features['temp_lag_1'] = features.groupby(['store', 'family'])['temp'].shift(1)

from sklearn.model_selection import KFold

def target_encode(train, target, col, n_folds=5):
    kf = KFold(n_splits=n_folds, shuffle=True, random_state=42)
    encoded = np.zeros(len(train))
    for train_idx, val_idx in kf.split(train):
        x_train = train.iloc[train_idx]
        x_val = train.iloc[val_idx]
        encoding_map = x_train.groupby(col)[target].mean()
        encoded[val_idx] = x_val[col].map(encoding_map)
    # 全局均值填充未见类别
    global_mean = train[target].mean()
    encoded = pd.Series(encoded).fillna(global_mean)
    return encoded

# 对store和family分别编码
features['store_te'] = target_encode(features, 'sales', 'store')
features['family_te'] = target_encode(features, 'sales', 'family')

train = features[features['date'] < '2025-07-01']
valid = features[(features['date'] >= '2025-07-01') & (features['date'] < '2025-10-01')]
test = features[features['date'] >= '2025-10-01']

# 特征列（剔除ID类、日期、目标）
exclude = ['date', 'store', 'family', 'sales', 'outlier']
X_train = train.drop(columns=exclude)
y_train = train['sales']
X_valid = valid.drop(columns=exclude)
y_valid = valid['sales']
X_test = test.drop(columns=exclude)
y_test = test['sales']

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

lr = LinearRegression()
lr.fit(X_train, y_train)
pred_lr = lr.predict(X_valid)
print(f"Linear Regression - MAE: {mean_absolute_error(y_valid, pred_lr):.2f}, RMSE: {np.sqrt(mean_squared_error(y_valid, pred_lr)):.2f}")

import xgboost as xgb
import lightgbm as lgb

xgb_model = xgb.XGBRegressor(n_estimators=1000, learning_rate=0.05, subsample=0.8, colsample_bytree=0.8, random_state=42)
xgb_model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], early_stopping_rounds=50, verbose=False)
pred_xgb = xgb_model.predict(X_valid)

lgb_model = lgb.LGBMRegressor(n_estimators=1000, learning_rate=0.05, subsample=0.8, colsample_bytree=0.8, random_state=42)
lgb_model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], callbacks=[lgb.early_stopping(50), lgb.log_evaluation(0)])
pred_lgb = lgb_model.predict(X_valid)

print(f"XGBoost - MAE: {mean_absolute_error(y_valid, pred_xgb):.2f}, RMSE: {np.sqrt(mean_squared_error(y_valid, pred_xgb)):.2f}")
print(f"LightGBM - MAE: {mean_absolute_error(y_valid, pred_lgb):.2f}, RMSE: {np.sqrt(mean_squared_error(y_valid, pred_lgb)):.2f}")

import optuna

def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 500, 2000),
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
        'subsample': trial.suggest_float('subsample', 0.6, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'reg_alpha': trial.suggest_float('reg_alpha', 1e-3, 10.0, log=True),
        'reg_lambda': trial.suggest_float('reg_lambda', 1e-3, 10.0, log=True),
        'min_child_weight': trial.suggest_int('min_child_weight', 1, 10),
        'random_state': 42
    }
    model = xgb.XGBRegressor(**params)
    model.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], early_stopping_rounds=50, verbose=False)
    pred = model.predict(X_valid)
    return np.sqrt(mean_squared_error(y_valid, pred))

study = optuna.create_study(direction='minimize', sampler=optuna.samplers.TPESampler(seed=42))
study.optimize(objective, n_trials=50, show_progress_bar=True)

best_params = study.best_params
print("Best params:", best_params)

best_xgb = xgb.XGBRegressor(**best_params, n_estimators=best_params['n_estimators'])
best_xgb.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], early_stopping_rounds=50, verbose=False)
pred_test = best_xgb.predict(X_test)
print(f"Test MAE: {mean_absolute_error(y_test, pred_test):.2f}, RMSE: {np.sqrt(mean_squared_error(y_test, pred_test)):.2f}")

import shap

explainer = shap.TreeExplainer(best_xgb)
shap_values = explainer.shap_values(X_valid.sample(1000, random_state=42))  # 抽样加快

# 汇总图
shap.summary_plot(shap_values, X_valid.sample(1000), plot_type="bar", show=False)
shap.summary_plot(shap_values, X_valid.sample(1000), show=False)

# 单个样本的力力图
shap.force_plot(explainer.expected_value, shap_values[0], X_valid.iloc[0], matplotlib=True, show=False)

import pickle

# 保存模型
with open('xgb_model.pkl', 'wb') as f:
    pickle.dump(best_xgb, f)

# 保存特征列名（用于输入校验）
feature_cols = X_train.columns.tolist()
with open('feature_cols.pkl', 'wb') as f:
    pickle.dump(feature_cols, f)

# app.py
from fastapi import FastAPI, HTTPException
from pydantic import BaseModel
import pandas as pd
import pickle
import numpy as np
from typing import List

app = FastAPI(title="Sales Prediction API")

# 加载模型和特征列
with open('xgb_model.pkl', 'rb') as f:
    model = pickle.load(f)
with open('feature_cols.pkl', 'rb') as f:
    feature_cols = pickle.load(f)

class PredictionInput(BaseModel):
    date: str
    store: int
    family: str
    onpromotion: int
    holiday: int
    temp: float
    # 注意：滞后特征等需要由调用方提供，或者我们在API内部根据历史数据计算（为简化，此处要求传入所有特征）

class PredictionResponse(BaseModel):
    predicted_sales: float

@app.post("/predict", response_model=PredictionResponse)
def predict(input_data: PredictionInput):
    try:
        # 将输入转为DataFrame，并补充所有特征（此处简化，仅做演示）
        # 实际中，我们需要根据传入的原始数据计算滞后和滚动特征，这需要历史数据存储。
        # 这里假设调用方已经计算好所有特征，我们只做校验。
        input_dict = input_data.dict()
        # 构建特征向量（需要与训练时顺序一致）
        # 实际场景中，应该有一个特征工程函数，这里略
        # 因为我们没有计算所有特征，仅演示结构
        # 所以我们硬编码一个示例特征向量（仅用于演示）
        # 生产环境需完整实现特征计算管道
        # 此处不展开，以防误导
        pass
    except Exception as e:
        raise HTTPException(status_code=400, detail=str(e))

# 实际部署时，建议使用 uvicorn app:app --host 0.0.0.0 --port 8000