基于Numpy实现时间序列学生化残差——函数编写问题咨询

Hey there! Let’s tackle this studentized residuals calculation problem you’re facing. I’ve dealt with similar numerical quirks before, so let’s break down where things might be going wrong and build a robust Python function together.

You’re using the correct formula for studentized residuals:

R_student(i) = e(i)/(std(i)*sqrt(1-h(i)))

But let’s clarify each component and avoid pitfalls that cause abnormal results:

• e(i): The raw residual for the i-th observation (actual y minus predicted y from the full model)
• std(i): The standard error of the model after removing the i-th observation (not the full model’s standard error!)
• h(i): The i-th diagonal element of the hat matrix H = X(XᵀX)⁻¹Xᵀ — this measures how influential the i-th observation is

The most common issues here are:

1. Numerical instability from direct matrix inversion: Calculating (XᵀX)⁻¹ directly can blow up if your design matrix X has multicollinearity or is nearly singular.
2. Incorrect degrees of freedom for std(i): When computing the leave-one-out MSE, you need to adjust the denominator correctly.
3. Forgetting the intercept in X: Your design matrix must include a column of 1s for the intercept term (as you noted in your X structure) — skipping this will make the hat matrix calculation wrong.

Instead of refitting the model n times (once per observation removed), we can use mathematical shortcuts to compute everything efficiently and stably. Here’s a function using NumPy:

import numpy as np

def calculate_studentized_residuals(X, y):
    # Ensure X includes the intercept column (first column of ones)
    n, p = X.shape
    if not np.allclose(X[:, 0], np.ones(n)):
        raise ValueError("X must include an intercept column (first column of ones)")
    
    # Fit the full linear model to get coefficients and raw residuals
    beta_hat = np.linalg.lstsq(X, y, rcond=None)[0]
    y_pred = X @ beta_hat
    raw_residuals = y - y_pred
    
    # Calculate hat matrix diagonal elements (h(i)) using pseudoinverse for stability
    XTX = X.T @ X
    XTX_pinv = np.linalg.pinv(XTX)  # Handles singular/near-singular XTX
    hat_diag = np.diag(X @ XTX_pinv @ X.T)
    
    # Sum of squared residuals for the full model
    sse_full = np.sum(raw_residuals ** 2)
    
    # Compute studentized residuals for each observation
    studentized_res = []
    for i in range(n):
        # Skip if h(i) is 1 (perfect leverage, which breaks the formula)
        if hat_diag[i] >= 1.0:
            studentized_res.append(np.nan)
            continue
        
        # Leave-one-out SSE (no need to refit the model!)
        sse_loo = sse_full - (raw_residuals[i] ** 2) / (1 - hat_diag[i])
        # Leave-one-out MSE: degrees of freedom = n - p - 1 (full model has p params, leave-one-out loses 1 observation)
        mse_loo = sse_loo / (n - p - 1)
        std_i = np.sqrt(mse_loo)
        
        # Compute the studentized residual
        r_student = raw_residuals[i] / (std_i * np.sqrt(1 - hat_diag[i]))
        studentized_res.append(r_student)
    
    return np.array(studentized_res)

Let’s run a quick test with sample data to verify:

# Sample data: X has intercept column, y is a simple linear relationship with noise
X = np.column_stack((np.ones(5), np.array([1, 2, 3, 4, 5])))
y = np.array([2, 4, 5, 4, 6])  # Added a slight outlier at the end

residuals = calculate_studentized_residuals(X, y)
print("Studentized Residuals:\n", residuals)

You should see that the outlier (last observation) has a larger absolute studentized residual, which is exactly what we expect.

• Pseudoinverse instead of direct inversion: np.linalg.pinv handles cases where XᵀX is singular or nearly singular, preventing numerical errors.
• Leave-one-out shortcut: We use the full model’s residuals to compute the leave-one-out SSE without refitting, which is faster and less error-prone.
• Leverage check: We catch cases where h(i) is 1 (a perfect leverage point) and return NaN instead of dividing by zero.

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