如何向量化OLS样本外误差计算循环？规避全局变量方案

初始实现代码 #

我经常遇到数组子集化循环处理的问题，以下是OLS样本外误差计算的初始实现：

import numpy as np

def get_betas(x, y):
    # Calculate parts of OLS formula
    xy = x.T @ y
    xx = x.T @ x
    xx_inv = np.linalg.inv(xx)
    
    # Calculate beta
    beta = xx_inv @ xy
    
    return beta

def calc_oos_error(x, y, k):
    errors = []
    
    for t in range(k, len(y)):
        beta = get_betas(x[:t], y[:t])
        error = y[t] - x[t] @ beta
        errors.append(error)
        
    return np.mean(errors), np.std(errors)

x = np.random.normal(size = 1000)
y = 10 * x + np.random.normal(scale = 1e-2, size = len(x))

x = x.reshape(-1, 1)

a, b = calc_oos_error(x, y, 10)

朴素向量化尝试的错误 #

尝试用np.vectorize直接向量化函数时会报错：

# Vectorize naive way
calc_oos_error_vec = np.vectorize(calc_oos_error)

# Try to run
calc_oos_error_vec(x = x, y = y, k = np.arange(1, 20))

错误信息：

TypeError: object of type 'numpy.float64' has no len().

问题根源是np.vectorize会把k数组的每个元素单独传入函数，导致函数内部的len(y)调用出错（此时y被当成单个元素而非数组）。

现有方案的缺陷 #

一种临时方案是使用全局变量引用x和y，但这种方式逻辑混乱，极易引发错误：

def calc(t):
    beta = get_betas(x[:t], y[:t])
    error = y[t] - x[t] @ beta
    return error

calc = np.vectorize(calc)    

def calc_oos_error_vec(x, y, k):
    errors = calc(np.arange(k, len(y)))
    return np.mean(errors), np.std(errors)

优化的向量化解决方案 #

利用单变量OLS的解析公式，通过累积求和实现完全向量化，避免循环和全局变量：

import numpy as np

def calc_oos_error_vec(x, y, k):
    x_flat = x.flatten()
    n = len(y)
    # 计算累积统计量：前t个样本的和、平方和、交叉乘积和
    cum_x = np.cumsum(x_flat)
    cum_y = np.cumsum(y)
    cum_xy = np.cumsum(x_flat * y)
    cum_x2 = np.cumsum(x_flat ** 2)
    
    # 确定需要计算的t范围（从k到n-1）
    t_vals = np.arange(k, n)
    sample_sizes = t_vals  # 每个t对应前t个样本（0到t-1），样本量为t
    
    # 基于累积统计量计算每个t对应的OLS系数beta
    numerator = cum_xy[t_vals - 1] - (cum_x[t_vals - 1] * cum_y[t_vals - 1]) / sample_sizes
    denominator = cum_x2[t_vals - 1] - (cum_x[t_vals - 1] ** 2) / sample_sizes
    beta = numerator / denominator
    
    # 计算每个t对应的样本外误差
    errors = y[t_vals] - x_flat[t_vals] * beta
    
    return np.mean(errors), np.std(errors)

# 测试用例
x = np.random.normal(size=1000)
y = 10 * x + np.random.normal(scale=1e-2, size=len(x))
x = x.reshape(-1, 1)

# 对比原方法和向量化方法的结果
mean_old, std_old = calc_oos_error(x, y, 10)
mean_vec, std_vec = calc_oos_error_vec(x, y, 10)

print(f"原方法均值：{mean_old:.6f}，标准差：{std_old:.6f}")
print(f"向量化方法均值：{mean_vec:.6f}，标准差：{std_vec:.6f}")

方案说明 #

1. 利用单变量OLS解析解：单变量线性回归的系数可以通过样本均值、协方差直接计算，无需矩阵求逆，简化计算逻辑
2. 累积求和替代切片循环：通过np.cumsum一次性计算所有前缀的统计量，避免每次循环切片数组
3. 完全向量化运算：所有计算步骤均为numpy数组操作，效率远高于循环，且无全局变量依赖

内容的提问来源于stack exchange，提问作者Charles0349