经典回归模型汇总 | Regression | SLie's Blog|琴弦之轮

线性回归_{Linear Regression}

待更：https://blog.csdn.net/longgb123/article/details/79079434
【计量经济学导论】03. 矩阵形式的线性回归模型_回归模型估计解矩阵-CSDN博客

广义线性回归

加权回归

递归最小二乘

递归最小二乘法.pdf (dezeming.top)
最优化算法-递推最小二乘法 - xiconxi - 博客园 (cnblogs.com)
自适应滤波：递归最小二乘 - LeeLIn。 - 博客园 (cnblogs.com)

import numpy as np
import matplotlib.pyplot as plt

class RLS:
    def __init__(self, beta=1.0, dim=200):
        self.beta = beta # 遗忘因子
        self.dim = dim # 数据维度/特征数

        self.theta = np.zeros(self.dim)  # 初始化参数均为 0

        self.P = np.eye(self.dim)  # 预测误差协方差矩阵的逆
        self.K = np.zeros(self.dim)  # 增益向量

    def update(self, X, y):
        XT_P_X = np.dot(X.T, np.dot(self.P, X)) + self.beta
        if XT_P_X != 0: # 避免除零错误
            self.K = np.dot(self.P, X) / XT_P_X
        else:
            self.K = np.zeros_like(X)  

        self.theta += self.K * (y - np.dot(X.T, self.theta))
        
        self.P = np.dot((np.eye(self.dim) - self.K.reshape(-1,1) @  X.reshape(1,-1)), self.P) / self.beta
        # self.K @ X.T  [易错] 这个代码得到的结果是一个值，而不是一个矩阵！ 应该要用 reshape

    def estimate(self, X, y):
        y_hat = X.T @ self.theta
        error = (y - y_hat) ** 2
        return error


def data_generate(num_samples, dim):
    '''根据样本量和特征数 生成数据'''

    _theta = np.random.normal(0, 1, dim)
    X = np.random.normal(0, 1, (num_samples, dim))
    eta = np.sqrt(0.01) * np.random.randn(num_samples)
    Y = np.dot(X, _theta) + eta

    return X, Y, _theta


if __name__ == '__main__':
    
    num_experiment = 100    # 实验次数
    beta = 1                # 遗忘因子
    num_features = 200      # 特征数
    num_samples = 1000     # 样本数/迭代次数

    db_avg_err = []         # 存储每次迭代的平均误差(dB)

    for _ in range(num_experiment): # 实验次数
        x_n, y_n, _theta = data_generate(num_samples, num_features)
        rls = RLS(beta=beta, dim=num_features)
        
        # 在线迭代更新
        errors = []
        for x, y in zip(x_n, y_n):
            rls.update(x, y)
            err = rls.estimate(x, y)
            errors.append(err)

        db_avg_err.append(errors)

    db_avg_err = np.mean(db_avg_err, axis=0)
    db_avg_err = 10 * np.log10(db_avg_err)

    plt.plot(range(num_samples), db_avg_err, label='Average Error (dB)')
    
    plt.xlabel('Iterations')
    plt.ylabel('Error (dB)')
    plt.title(f'RLS Convergence  ($\\beta={beta}$)')
    plt.legend()
    plt.show()

    # 简单查看拟合效果    
    print(_theta[:5])
    print(rls.theta[:5])