Python实现迭代重加权最小二乘法（IRLS）

import numpy as np
import matplotlib.pyplot as plt

def irls(x, y, degree, iterations, epsilon):
    # 初始化参数
    theta = np.zeros(degree + 1)
    m = len(x)
    
    # 添加偏置项
    X = np.column_stack([np.power(x, i) for i in range(degree + 1)])
    
    # 迭代更新参数
    for _ in range(iterations):
        # 计算预测值
        y_pred = np.dot(X, theta)
        
        # 计算残差
        residuals = y - y_pred
        
        # 计算权重矩阵
        weights = np.diag(1 / (np.abs(residuals) + epsilon))
        
        # 更新参数
        theta = np.linalg.inv(X.T @ weights @ X) @ X.T @ weights @ y
    
    return theta

# 生成样本数据
x = np.linspace(0, 10, 100)
y = 2 * x + 1 + np.random.randn(100)

# 使用IRLS拟合曲线
degree = 1  # 多项式的次数
iterations = 10  # 迭代次数
epsilon = 0.01  # 避免除零错误的小值
theta = irls(x, y, degree, iterations, epsilon)

# 绘制拟合曲线
x_range = np.linspace(0, 10, 100)
y_pred = np.dot(np.column_stack([np.power(x_range, i) for i in range(degree + 1)]), theta)
plt.scatter(x, y, label='Data')
plt.plot(x_range, y_pred, color='red', label='Fitted curve')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

在上述代码中，我们定义了一个irls函数，该函数接受输入的x和y数据，以及多项式的次数、迭代次数和epsilon（用于避免除零错误的小值）作为参数。函数内部首先初始化参数theta为全零向量，并将x的幂次方作为特征矩阵X。然后，使用迭代的方式更新参数theta，直到达到指定的迭代次数。在每次迭代中，计算预测值y_pred、残差residuals和权重矩阵weights，并使用加权最小二乘法更新参数theta。最后，返回更新后的参数theta