MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题
MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题
多元自适应回归样条（Multivariate Adaptive Regression Splines, MARS）是Jerome Friedman于1991年提出的一种非参数回归技术。该方法专门用于建模预测变量集合与目标变量之间的复杂非线性关系，无需预先确定具体的函数形式。本文将深入探讨MARS算法的核心原理，并详细阐述其在时间序列预测任务中的应用策略与技术实现。
线性回归模型的局限性分析
在理解MARS算法的优势之前，有必要首先分析传统线性模型所面临的技术挑战。当提及"线性模型"时，通常指的是线性回归方法。
线性回归模型存在三个根本性的技术限制。首先，该模型严格假设变量间存在线性关系，因此无法有效捕获数据中的曲线特征、阈值效应或斜率突变现象。其次，线性回归对多重共线性问题极为敏感，当预测变量之间存在高度相关性时，会导致回归系数估计的不稳定性和可靠性下降。第三，该方法在处理复杂特征交互时存在明显不足，需要研究者预先明确指定交互项，这在缺乏先验知识的情况下往往难以实现。
现实世界中的大多数问题都涉及非线性关系，这使得传统线性回归的应用受到显著限制。
多项式回归的技术挑战
为了解决线性假设的局限性，一种常见的方法是引入多项式项来处理非线性关系。当自变量X与因变量y之间存在曲线关系时，可以通过添加多项式项来扩展基础线性模型。
采用全局高次多项式来解决线性假设问题往往会引发新的技术困难。全局多项式具有系统性影响，即曲线任何局部的调整都会影响整体拟合效果。高次多项式容易产生过拟合现象，在数据边界处可能出现极端振荡。此外，多项式方法缺乏自适应性，需要研究者手动选择多项式的次数，这种选择往往缺乏客观标准。
MARS算法的核心技术原理
MARS算法采用了一种根本不同的策略来处理非线性建模问题。该算法不是在整个数据集上构建单一的复杂多项式，而是采用分而治之的方法。MARS通过节点（断点）将预测变量空间划分为多个局部区域，在每个区域内拟合简单的低次多项式（通常为一次线性函数）。算法允许每个变量在不同位置设置独立的节点，并且仅在对模型性能有显著提升时才添加变量间的交互项。
以具体示例说明，对于某个变量X，MARS可能识别出在X < 3的区域内存在一种斜率关系，在3 ≤ X < 5的区间内存在另一种斜率关系，而在X ≥ 5的区域内又呈现第三种斜率模式。
import numpy as np
import matplotlib.pyplot as plt
# 定义MARS风格的分段函数：
# y = 5 + 2 * max(0, X - 3) - 1.5 * max(0, 5 - X)
def mars_like(x):
x = np.asarray(x)
return 5 + 2*np.maximum(0, x - 3) - 1.5*np.maximum(0, 5 - x)
# 生成x值并计算y
X = np.linspace(0, 10, 400)
Y = mars_like(X)
# 绘图
plt.figure(figsize=(9, 5))
plt.plot(X, Y, label=r"$y = 5 + 2\max(0, x-3) - 1.5\max(0, 5-x)$")
plt.axvline(3, linestyle="--", alpha=0.6)
plt.axvline(5, linestyle="--", alpha=0.6)
# 每个区域的斜率注释
plt.text(1.2, mars_like(1.2)+0.2, "Region 1: x < 3\nslope = +0.0", fontsize=10)
plt.text(3.4, mars_like(3.4)+0.2, "Region 2: 3 ≤ x < 5\nslope = +2.0", fontsize=10)
plt.text(6.0, mars_like(6.0)+0.2, "Region 3: x ≥ 5\nslope = +2.0 - 1.5 = +0.5", fontsize=10)
plt.title("Piecewise Function (MARS-style) with Knots at x=3 and x=5")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.tight_layout()
plt.show()
不同建模方法的比较分析
为了更好地理解MARS算法的技术优势，下面通过对比实验来展示线性回归、多项式回归和MARS方法在处理非线性数据时的不同表现。
import matplotlib.pyplot as plt
import numpy as np
# 生成示例数据
np.random.seed(0)
X = np.linspace(0, 10, 100)
y_true = np.piecewise(X, [X < 3, (X >= 3) & (X < 7), X >= 7],
[lambda x: 2 + 0.5*x,
lambda x: 4 + 1.5*(x - 3),
lambda x: 10 - 0.8*(x - 7)])
y_noisy = y_true + np.random.normal(scale=0.8, size=len(X))
# 拟合简单线性回归
coef = np.polyfit(X, y_noisy, 1)
y_linear = np.polyval(coef, X)
# 拟合多项式回归（3次）
coef_poly = np.polyfit(X, y_noisy, 3)
y_poly = np.polyval(coef_poly, X)
# 模拟MARS拟合（分段线性，已知节点用于说明）
knots = [3, 7]
y_mars = np.piecewise(X, [X < knots[0], (X >= knots[0]) & (X < knots[1]), X >= knots[1]],
[lambda x: np.polyval(np.polyfit(X[X<3], y_noisy[X<3], 1), x),
lambda x: np.polyval(np.polyfit(X[(X>=3)&(X<7)], y_noisy[(X>=3)&(X<7)], 1), x),
lambda x: np.polyval(np.polyfit(X[X>=7], y_noisy[X>=7], 1), x)])
# 创建子图
fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)
# 线性
axes[0].scatter(X, y_noisy, color='gray', alpha=0.6, label="Data")
axes[0].plot(X, y_linear, color='red', label="Linear Fit")
axes[0].set_title("Linear Regression")
axes[0].legend()
# 多项式
axes[1].scatter(X, y_noisy, color='gray', alpha=0.6, label="Data")
axes[1].plot(X, y_poly, color='blue', label="Polynomial Fit (deg=3)")
axes[1].set_title("Polynomial Regression")
axes[1].legend()
# MARS
axes[2].scatter(X, y_noisy, color='gray', alpha=0.6, label="Data")
axes[2].plot(X, y_mars, color='green', label="MARS Fit (piecewise)")
for knot in knots:
axes[2].axvline(knot, color='black', linestyle='--', alpha=0.5)
axes[2].set_title("MARS (Piecewise Linear)")
axes[2].legend()
fig.suptitle("From Linear to Polynomial to MARS", fontsize=16)
plt.tight_layout()
plt.show()
从比较结果可以看出，线性回归方法强制使用单一直线拟合，这导致对非线性趋势的欠拟合问题。多项式回归虽然通过全局曲线弯曲提供了更大的灵活性，但存在过拟合风险和外推性能差的问题。相比之下MARS方法通过自动选择节点位置实现分段线性拟合，在有效控制模型复杂度的同时实现对数据的局部自适应建模。
MARS算法的工作机制
MARS算法的核心工作原理可以概括为自动识别数据行为变化的关键位置（节点），在这些位置周围构建分段线性基函数，然后通过模型简化过程保留最重要的变化模式。
在MARS算法中，节点的概念类似于机械中的铰链结构，在节点前后，拟合曲线可以呈现不同的斜率特征。基函数是MARS算法的核心构建要素，对于位置c处的每个节点，算法会创建两个对应的基函数，这种设计使得MARS能够独立建模节点两侧的不同行为模式。两个基函数的数学表达式为：BF1 = max(0,x−c) 和 BF2 = max(0,c−x)。
MARS算法的拟合过程包含两个关键阶段。前向传递阶段是模型构建的增长过程，算法从单一的平坦基线开始，系统性地搜索所有变量和可能的节点位置，识别出能够最大化改善模型拟合效果的基函数对，这个过程持续进行直至满足预设的停止条件。后向传递阶段也称为剪枝过程，在此阶段MARS使用广义交叉验证（GCV）准则逐个移除对模型贡献最小的基函数，最终得到一个结构简化且具有良好泛化能力的模型。
MARS在时间序列建模中的应用
虽然MARS可以应用于时间序列分析，但需要注意的是，它本质上并非专门的时间序列模型，不具备ARIMA或LSTM等模型固有的时间依赖性概念。因此在将MARS应用于时间序列预测时，需要对数据进行适当的预处理，使其能够有效捕获时间模式。
由于MARS算法处理的是表格化数据格式（X，y），首要任务是将时间序列问题转换为监督学习框架。
Python实现与技术细节
以下是MARS算法在时间序列预测中的完整技术实现：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import List, Tuple, Optional
np.random.seed(123)
# 1) 生成合成非线性季节性时间序列
n = 360 # 12个月 * 30天作为示例长度
t = np.arange(n)
# 真实过程：趋势 + 月度季节性 + 慢速季节性机制 + 阈值非线性
trend = 0.03 * t
monthly = 8 * np.sin(2 * np.pi * t / 30) # 30天周期
season90 = 3 * np.sign(np.sin(2 * np.pi * t / 90)) # 大约每45天机制翻转
threshold_boost = np.where(monthly > 6, 4, 0) # 阈值非线性
y_true = 20 + trend + monthly + season90 + threshold_boost
y = y_true + np.random.normal(scale=1.2, size=n)
# 2) 框架化为监督学习
def make_lag_features(series: np.ndarray, max_lag: int) -> pd.DataFrame:
s = pd.Series(series)
data = {f"lag_{k}": s.shift(k) for k in range(1, max_lag + 1)}
return pd.DataFrame(data)
max_lag = 10
df = make_lag_features(y, max_lag)
# 日历/季节性特征
df["sin_30"] = np.sin(2 * np.pi * np.arange(len(df)) / 30)
df["cos_30"] = np.cos(2 * np.pi * np.arange(len(df)) / 30)
df["sin_90"] = np.sin(2 * np.pi * np.arange(len(df)) / 90)
df["cos_90"] = np.cos(2 * np.pi * np.arange(len(df)) / 90)
df["time_idx"] = np.arange(len(df))
df["target"] = y
df = df.dropna().reset_index(drop=True)
# 按时间顺序分割
train_frac = 0.8
split_idx = int(len(df) * train_frac)
X_train = df.drop(columns=["target"]).iloc[:split_idx].to_numpy()
y_train = df["target"].iloc[:split_idx].to_numpy()
X_test = df.drop(columns=["target"]).iloc[split_idx:].to_numpy()
y_test = df["target"].iloc[split_idx:].to_numpy()
feature_names = df.drop(columns=["target"]).columns.tolist()
# 3) 轻量级MARS类模型（铰链基 + 线性项），前向选择 + 后向剪枝
def hinge(vec: np.ndarray, c: float, direction: int) -> np.ndarray:
# direction: +1 => max(0, x - c), -1 => max(0, c - x)
if direction == 1:
return np.maximum(0.0, vec - c)
else:
return np.maximum(0.0, c - vec)
@dataclass
class BasisSpec:
kind: str # "hinge" 或 "linear"
feat_idx: int
knot: Optional[float] # 线性为None
direction: Optional[int] # 铰链为+1/-1，线性为None
@dataclass
class MarsLiteTS:
basis: List[BasisSpec]
beta: np.ndarray
def design_matrix(X: np.ndarray, basis_list: List[BasisSpec]) -> np.ndarray:
parts = [np.ones((X.shape[0], 1))] # 截距
for b in basis_list:
col = X[:, b.feat_idx]
if b.kind == "linear":
vec = col
else:
vec = hinge(col, b.knot, b.direction)
parts.append(vec.reshape(-1, 1))
return np.hstack(parts)
def gcv(y_true: np.ndarray, y_hat: np.ndarray, m_terms: int, penalty: float = 2.0) -> float:
n = len(y_true)
rss = np.sum((y_true - y_hat) ** 2)
# 有效参数数量：1（截距）+ m_terms，按惩罚缩放
c_eff = 1 + penalty * m_terms
denom = max(1e-3, (1 - c_eff / n))
return rss / (n * denom ** 2)
def fit_mars_lite_ts(X: np.ndarray, y: np.ndarray, max_terms: int = 20, q_knots: int = 15,
allow_linear: bool = True, penalty: float = 2.0, min_improve: float = 1e-3) -> MarsLiteTS:
n, p = X.shape
# 从每个特征的分位数获得候选节点
quantiles = np.linspace(0.05, 0.95, q_knots)
cand_knots = [np.quantile(X[:, j], quantiles) for j in range(p)]
basis_list: List[BasisSpec] = []
# 仅截距基线
D = design_matrix(X, basis_list)
beta, *_ = np.linalg.lstsq(D, y, rcond=None)
yhat = D @ beta
best_score = gcv(y, yhat, m_terms=0, penalty=penalty)
# 前向选择
for _ in range(max_terms):
best_add = None
best_add_score = best_score
# 尝试线性项
if allow_linear:
for j in range(p):
spec = BasisSpec(kind="linear", feat_idx=j, knot=None, direction=None)
D_try = design_matrix(X, basis_list + [spec])
beta_try, *_ = np.linalg.lstsq(D_try, y, rcond=None)
yhat_try = D_try @ beta_try
score = gcv(y, yhat_try, m_terms=len(basis_list) + 1, penalty=penalty)
if score < best_add_score - 1e-12:
best_add_score = score
best_add = (spec, beta_try, yhat_try)
# 尝试铰链项
for j in range(p):
for c in cand_knots[j]:
for d in (1, -1):
spec = BasisSpec(kind="hinge", feat_idx=j, knot=float(c), direction=d)
D_try = design_matrix(X, basis_list + [spec])
beta_try, *_ = np.linalg.lstsq(D_try, y, rcond=None)
yhat_try = D_try @ beta_try
score = gcv(y, yhat_try, m_terms=len(basis_list) + 1, penalty=penalty)
if score < best_add_score - 1e-12:
best_add_score = score
best_add = (spec, beta_try, yhat_try)
if best_add is None or (best_score - best_add_score) < min_improve:
break
# 接受添加
spec, beta, yhat = best_add
basis_list.append(spec)
best_score = best_add_score
# 后向剪枝（贪心）
improved = True
while improved and len(basis_list) > 0:
improved = False
current_best = best_score
drop_idx = None
for j in range(len(basis_list)):
trial = [b for i, b in enumerate(basis_list) if i != j]
D_try = design_matrix(X, trial)
beta_try, *_ = np.linalg.lstsq(D_try, y, rcond=None)
yhat_try = D_try @ beta_try
score = gcv(y, yhat_try, m_terms=len(trial), penalty=penalty)
if score < current_best - 1e-12:
current_best = score
drop_idx = j
best_beta_try = beta_try
best_yhat_try = yhat_try
if drop_idx is not None:
basis_list.pop(drop_idx)
beta = best_beta_try
yhat = best_yhat_try
best_score = current_best
improved = True
# 最终重新拟合
D_final = design_matrix(X, basis_list)
beta_final, *_ = np.linalg.lstsq(D_final, y, rcond=None)
return MarsLiteTS(basis=basis_list, beta=beta_final)
def predict(model: MarsLiteTS, X: np.ndarray) -> np.ndarray:
D = design_matrix(X, model.basis)
return D @ model.beta
# 训练模型
model = fit_mars_lite_ts(X_train, y_train, max_terms=25, q_knots=12, allow_linear=True, penalty=2.0, min_improve=1e-3)
# 评估
y_pred = predict(model, X_test)
rmse = float(np.sqrt(np.mean((y_test - y_pred) ** 2)))
mape = float(np.mean(np.abs((y_test - y_pred) / np.maximum(1e-6, np.abs(y_test)))) * 100)
# 4) 绘图
plt.figure(figsize=(12, 6))
plt.plot(y, label="Actual")
plt.axvline(split_idx + max_lag, linestyle="--", label="Train/Test split")
# 将预测对齐到原始索引
start = split_idx + max_lag
x_axis = np.arange(start, start + len(y_pred))
plt.plot(x_axis, y_pred, label="MARS-like forecast")
plt.title(f"MARS-like Time Series Forecast | RMSE: {rmse:.2f}, MAPE: {mape:.1f}%")
plt.xlabel("Time")
plt.ylabel("Value")
plt.legend()
plt.tight_layout()
plt.show()
MARS算法的适用场景分析
MARS算法在时间序列建模中具有独特的技术优势，特别适用于以下几类应用场景。
当研究者怀疑数据中存在阈值效应时，MARS算法能够有效识别和建模这类非连续性变化。例如，在销售预测中，当温度超过30°C时销售量可能出现显著跃升，这种阈值驱动的行为模式正是MARS算法擅长处理的问题类型。
对于存在制度转换的时间序列数据，MARS算法能够自动检测变化点并在转换点后适应新的趋势斜率。这种能力使得MARS在处理结构性变化的经济时间序列或具有季节性转换的自然现象数据时表现优异。
当系统受到多个非线性外部驱动因素影响时，MARS算法通过其分段建模能力能够有效捕获不同变量在不同取值范围内的差异化影响模式，从而提供比传统线性模型更准确的预测结果。
作者：Okan Yenigün
本文分享自 DeepHub IMBA 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ，欢迎热爱写作的你一起参与！
腾讯云开发者
