机器学习基础系列——线性回归

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# 真实参数
true_w = np.array([2.0, -3.5])
true_b = 4.2

# 生成数据
n_samples = 100
X = np.random.randn(n_samples, 2)
y = X @ true_w + true_b + np.random.randn(n_samples) * 0.5

print(f"X shape: {X.shape}")  # (100, 2)
print(f"y shape: {y.shape}")  # (100,)

def linear_regression_normal(X, y):
    """使用正规方程求解"""
    # 添加偏置项
    X_b = np.c_[np.ones((len(X), 1)), X]
    
    # 正规方程
    theta = np.linalg.inv(X_b.T @ X_b) @ X_b.T @ y
    
    return theta[0], theta[1:]  # b, w

b, w = linear_regression_normal(X, y)
print(f"估计的 w: {w}")
print(f"估计的 b: {b:.4f}")
print(f"真实的 w: {true_w}")
print(f"真实的 b: {true_b}")

class LinearRegressionGD:
    def __init__(self, learning_rate=0.01, n_iterations=1000):
        self.lr = learning_rate
        self.n_iterations = n_iterations
        self.w = None
        self.b = None
        self.losses = []
    
    def fit(self, X, y):
        n_samples, n_features = X.shape
        
        # 初始化参数
        self.w = np.zeros(n_features)
        self.b = 0
        
        # 梯度下降
        for _ in range(self.n_iterations):
            # 预测
            y_pred = X @ self.w + self.b
            
            # 计算损失
            loss = np.mean((y - y_pred) ** 2) / 2
            self.losses.append(loss)
            
            # 计算梯度
            dw = -(1 / n_samples) * X.T @ (y - y_pred)
            db = -(1 / n_samples) * np.sum(y - y_pred)
            
            # 更新参数
            self.w -= self.lr * dw
            self.b -= self.lr * db
        
        return self
    
    def predict(self, X):
        return X @ self.w + self.b

# 训练
model = LinearRegressionGD(learning_rate=0.1, n_iterations=1000)
model.fit(X, y)

print(f"GD 估计的 w: {model.w}")
print(f"GD 估计的 b: {model.b:.4f}")

# 绘制损失曲线
plt.figure(figsize=(10, 4))
plt.plot(model.losses)
plt.xlabel('迭代次数')
plt.ylabel('损失')
plt.title('训练损失曲线')
plt.show()

class StandardScaler:
    def fit(self, X):
        self.mean = np.mean(X, axis=0)
        self.std = np.std(X, axis=0)
        return self
    
    def transform(self, X):
        return (X - self.mean) / self.std
    
    def fit_transform(self, X):
        return self.fit(X).transform(X)

# 使用
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 标准化后训练
model_scaled = LinearRegressionGD(learning_rate=0.1, n_iterations=100)
model_scaled.fit(X_scaled, y)

class MinMaxScaler:
    def fit(self, X):
        self.min = np.min(X, axis=0)
        self.max = np.max(X, axis=0)
        return self
    
    def transform(self, X):
        return (X - self.min) / (self.max - self.min)

class PolynomialFeatures:
    def __init__(self, degree=2):
        self.degree = degree
    
    def transform(self, X):
        n_samples = len(X)
        features = [np.ones(n_samples)]
        
        for d in range(1, self.degree + 1):
            features.append(X ** d)
        
        return np.column_stack(features)

# 生成非线性数据
X_nonlinear = np.linspace(-3, 3, 100)
y_nonlinear = 0.5 * X_nonlinear**2 + X_nonlinear + 2 + np.random.randn(100) * 0.5

# 多项式特征
poly = PolynomialFeatures(degree=2)
X_poly = poly.transform(X_nonlinear)

# 正规方程求解
theta = np.linalg.lstsq(X_poly, y_nonlinear, rcond=None)[0]
print(f"多项式系数: {theta}")

# 预测
y_pred = X_poly @ theta

plt.scatter(X_nonlinear, y_nonlinear, alpha=0.5, label='数据')
plt.plot(X_nonlinear, y_pred, 'r-', label='拟合曲线')
plt.legend()
plt.show()

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

# 划分数据
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# 方法1：基本线性回归
model = LinearRegression()
model.fit(X_train, y_train)

print(f"权重: {model.coef_}")
print(f"偏置: {model.intercept_:.4f}")

# 评估
y_pred = model.predict(X_test)
print(f"MSE: {mean_squared_error(y_test, y_pred):.4f}")
print(f"R²: {r2_score(y_test, y_pred):.4f}")

# 方法2：Pipeline
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('regressor', LinearRegression())
])

pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)
print(f"Pipeline R²: {r2_score(y_test, y_pred):.4f}")

def evaluate_regression(y_true, y_pred):
    mse = np.mean((y_true - y_pred) ** 2)
    rmse = np.sqrt(mse)
    mae = np.mean(np.abs(y_true - y_pred))
    
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    r2 = 1 - ss_res / ss_tot
    
    return {
        'MSE': mse,
        'RMSE': rmse,
        'MAE': mae,
        'R²': r2
    }

metrics = evaluate_regression(y_test, y_pred)
for name, value in metrics.items():
    print(f"{name}: {value:.4f}")

residuals = y_test - y_pred

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# 残差分布
axes[0].hist(residuals, bins=20, edgecolor='black')
axes[0].axvline(0, color='r', linestyle='--')
axes[0].set_xlabel('残差')
axes[0].set_ylabel('频数')
axes[0].set_title('残差分布')

# 残差 vs 预测值
axes[1].scatter(y_pred, residuals, alpha=0.5)
axes[1].axhline(0, color='r', linestyle='--')
axes[1].set_xlabel('预测值')
axes[1].set_ylabel('残差')
axes[1].set_title('残差 vs 预测值')

plt.tight_layout()
plt.show()