一元線性回歸

# -*- coding: utf-8 -*-
"""
一元線性回歸（Linear Regression）正式教材用範例（繁中超詳細註解）
====================================================================

本程式示範兩種求解方式：
(1) 封閉解（Closed-form / 最小平方法）：一次算出最佳 b0, b1
(2) 梯度下降（Gradient Descent）：模擬「訓練」過程，逐步更新 b0, b1

核心模型：
    y_hat = b0 + b1*x

核心 Loss（MSE，均方誤差）：
    MSE = (1/n) * Σ (y - y_hat)^2

核心評估（R^2，決定係數）：
    R^2 = 1 - SS_res/SS_tot

------------------------------------------------------------
Windows 輸出編碼的保險
------------------------------------------------------------
- Windows 主控台可能使用 cp950
- 某些特殊符號可能無法輸出，造成 UnicodeEncodeError
- 本程式採取兩個策略：
  (1) 嘗試把 stdout 改成 UTF-8（errors="replace" 避免炸裂）
  (2) 輸出避免使用高風險符號（例如改用 ~ 表示約略）
"""

from __future__ import annotations

import sys
import numpy as np

# -----------------------------
# A) 輸出編碼保險：盡量用 UTF-8，避免主控台編碼炸裂
# -----------------------------
try:
    if hasattr(sys.stdout, "reconfigure"):
        sys.stdout.reconfigure(encoding="utf-8", errors="replace")
except Exception:
    pass

# -----------------------------
# B) matplotlib：可選畫圖
# -----------------------------
try:
    import matplotlib.pyplot as plt
    HAS_MATPLOTLIB = True
except ImportError:
    HAS_MATPLOTLIB = False


def fit_closed_form(x: np.ndarray, y: np.ndarray) -> tuple[float, float]:
    """
    封閉解（最小平方法）：
        b1 = Σ((x - x̄)(y - ȳ)) / Σ((x - x̄)^2)
        b0 = ȳ - b1*x̄
    """
    x_mean = x.mean()
    y_mean = y.mean()

    numerator = np.sum((x - x_mean) * (y - y_mean))
    denominator = np.sum((x - x_mean) ** 2)

    b1 = numerator / denominator
    b0 = y_mean - b1 * x_mean
    return b0, b1


def fit_gradient_descent(
    x: np.ndarray,
    y: np.ndarray,
    lr: float = 0.02,
    epochs: int = 2000
) -> tuple[float, float, list[float]]:
    """
    梯度下降訓練（Gradient Descent）

    梯度（結論）：
        d(MSE)/d(b0) = (-2/n) * Σ( y - (b0 + b1*x) )
        d(MSE)/d(b1) = (-2/n) * Σ( x * (y - (b0 + b1*x)) )

    更新：
        b0 = b0 - lr * grad_b0
        b1 = b1 - lr * grad_b1
    """
    n = len(x)
    b0 = 0.0
    b1 = 0.0
    losses: list[float] = []

    for epoch in range(1, epochs + 1):
        y_pred = b0 + b1 * x
        mse = float(np.mean((y - y_pred) ** 2))
        losses.append(mse)

        error = y - y_pred
        grad_b0 = (-2 / n) * float(np.sum(error))
        grad_b1 = (-2 / n) * float(np.sum(x * error))

        b0 = b0 - lr * grad_b0
        b1 = b1 - lr * grad_b1

        if epoch % 400 == 0:
            print(f"[GD] epoch={epoch:4d}  MSE={mse:.4f}  b0={b0:.4f}  b1={b1:.4f}")

    return b0, b1, losses


def r2_score(y_true: np.ndarray, y_pred: np.ndarray) -> float:
    ss_res = float(np.sum((y_true - y_pred) ** 2))
    ss_tot = float(np.sum((y_true - y_true.mean()) ** 2))
    return 1 - (ss_res / ss_tot)


def main() -> None:
    x_list = [1, 2, 3, 4, 5, 6]
    y_list = [52, 55, 61, 66, 72, 75]
    x = np.array(x_list, dtype=float)
    y = np.array(y_list, dtype=float)

    b0_cf, b1_cf = fit_closed_form(x, y)
    y_pred_cf = b0_cf + b1_cf * x
    r2_cf = r2_score(y, y_pred_cf)

    print("\n=== Closed-form (Least Squares) ===")
    print(f"b0={b0_cf:.4f}, b1={b1_cf:.4f}, R^2={r2_cf:.4f}")
    print(f"Model: y = {b0_cf:.4f} + {b1_cf:.4f} * x")

    b0_gd, b1_gd, losses = fit_gradient_descent(x, y, lr=0.02, epochs=2000)
    y_pred_gd = b0_gd + b1_gd * x
    r2_gd = r2_score(y, y_pred_gd)

    print("\n=== Gradient Descent Training ===")
    print(f"b0={b0_gd:.4f}, b1={b1_gd:.4f}, R^2={r2_gd:.4f}")
    print(f"Model: y = {b0_gd:.4f} + {b1_gd:.4f} * x")

    x_new = 7
    y_new_pred = b0_cf + b1_cf * x_new
    print("\nPrediction (using closed-form):")
    print(f"x={x_new} -> y~{y_new_pred:.2f}")

    if HAS_MATPLOTLIB:
        plt.scatter(x, y, label="Data (x,y)")
        plt.plot(x, y_pred_cf, label="Closed-form line")
        plt.plot(x, y_pred_gd, label="GD line")
        plt.title("Simple Linear Regression Demo")
        plt.xlabel("x")
        plt.ylabel("y")
        plt.legend()
        plt.show()

        plt.plot(np.arange(1, len(losses) + 1), losses, label="GD Loss (MSE)")
        plt.title("Gradient Descent Training Loss")
        plt.xlabel("epoch")
        plt.ylabel("MSE")
        plt.legend()
        plt.show()
    else:
        print("\n(matplotlib not installed, skip plotting)  Install: pip install matplotlib")


if __name__ == "__main__":
    main()