MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题-优快云博客

多元自适应回归样条(Multivariate Adaptive Regression Splines, MARS)是Jerome Friedman于1991年提出的一种非参数回归技术。该方法专门用于建模预测变量集合与目标变量之间的复杂非线性关系，无需预先确定具体的函数形式。本文将深入探讨MARS算法的核心原理，并详细阐述其在时间序列预测任务中的应用策略与技术实现。

线性回归模型的局限性分析

在理解MARS算法的优势之前，有必要首先分析传统线性模型所面临的技术挑战。当提及"线性模型"时，通常指的是线性回归方法。

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#人工智能

线性回归模型存在三个根本性的技术限制。首先，该模型严格假设变量间存在线性关系，因此无法有效捕获数据中的曲线特征、阈值效应或斜率突变现象。其次，线性回归对多重共线性问题极为敏感，当预测变量之间存在高度相关性时，会导致回归系数估计的不稳定性和可靠性下降。第三，该方法在处理复杂特征交互时存在明显不足，需要研究者预先明确指定交互项，这在缺乏先验知识的情况下往往难以实现。

现实世界中的大多数问题都涉及非线性关系，这使得传统线性回归的应用受到显著限制。

多项式回归的技术挑战

为了解决线性假设的局限性，一种常见的方法是引入多项式项来处理非线性关系。当自变量X与因变量y之间存在曲线关系时，可以通过添加多项式项来扩展基础线性模型。

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#回归_02

采用全局高次多项式来解决线性假设问题往往会引发新的技术困难。全局多项式具有系统性影响，即曲线任何局部的调整都会影响整体拟合效果。高次多项式容易产生过拟合现象，在数据边界处可能出现极端振荡。此外，多项式方法缺乏自适应性，需要研究者手动选择多项式的次数，这种选择往往缺乏客观标准。

MARS算法的核心技术原理

MARS算法采用了一种根本不同的策略来处理非线性建模问题。该算法不是在整个数据集上构建单一的复杂多项式，而是采用分而治之的方法。MARS通过节点(断点)将预测变量空间划分为多个局部区域，在每个区域内拟合简单的低次多项式(通常为一次线性函数)。算法允许每个变量在不同位置设置独立的节点，并且仅在对模型性能有显著提升时才添加变量间的交互项。

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#分段回归_03

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#分段回归_04

以具体示例说明，对于某个变量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算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#时间序列_05

不同建模方法的比较分析

为了更好地理解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算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#回归_06

从比较结果可以看出，线性回归方法强制使用单一直线拟合，这导致对非线性趋势的欠拟合问题。多项式回归虽然通过全局曲线弯曲提供了更大的灵活性，但存在过拟合风险和外推性能差的问题。相比之下MARS方法通过自动选择节点位置实现分段线性拟合，在有效控制模型复杂度的同时实现对数据的局部自适应建模。

MARS算法的工作机制

MARS算法的核心工作原理可以概括为自动识别数据行为变化的关键位置(节点)，在这些位置周围构建分段线性基函数，然后通过模型简化过程保留最重要的变化模式。

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#分段回归_07

在MARS算法中，节点的概念类似于机械中的铰链结构，在节点前后，拟合曲线可以呈现不同的斜率特征。基函数是MARS算法的核心构建要素，对于位置c处的每个节点，算法会创建两个对应的基函数，这种设计使得MARS能够独立建模节点两侧的不同行为模式。两个基函数的数学表达式为：BF1 = max(0,x−c) 和 BF2 = max(0,c−x)。

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#机器学习_08

MARS算法的拟合过程包含两个关键阶段。前向传递阶段是模型构建的增长过程，算法从单一的平坦基线开始，系统性地搜索所有变量和可能的节点位置，识别出能够最大化改善模型拟合效果的基函数对，这个过程持续进行直至满足预设的停止条件。后向传递阶段也称为剪枝过程，在此阶段MARS使用广义交叉验证(GCV)准则逐个移除对模型贡献最小的基函数，最终得到一个结构简化且具有良好泛化能力的模型。

MARS在时间序列建模中的应用

虽然MARS可以应用于时间序列分析，但需要注意的是，它本质上并非专门的时间序列模型，不具备ARIMA或LSTM等模型固有的时间依赖性概念。因此在将MARS应用于时间序列预测时，需要对数据进行适当的预处理，使其能够有效捕获时间模式。

由于MARS算法处理的是表格化数据格式(X，y)，首要任务是将时间序列问题转换为监督学习框架。

MARS算法理论和Python代码实现：用分段回归解决非线性时间序列预测问题_#机器学习_09

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()