6.4 最小平方法與線性變換：當方程組無解，找最好的近似解#

from sympy import Matrix, cos, sin, pi, simplify

theta = pi / 6   # 30 度
R = Matrix([[cos(theta), -sin(theta)],
            [sin(theta),  cos(theta)]])
print(simplify(R * Matrix([1, 0])))   # Matrix([[sqrt(3)/2], [1/2]])：正確！
print(simplify(R.T * R))              # 單位矩陣：旋轉不改變長度

M = Matrix([[cos(2*theta), sin(2*theta)],
            [sin(2*theta), -cos(2*theta)]])
print(simplify(M * Matrix([1, 0])))   # 鏡射 (1,0)：結果在 60 度方向

from sympy import Matrix, pi, cos, sin, simplify

theta = pi / 2
R = Matrix([[cos(theta), -sin(theta)],
            [sin(theta),  cos(theta)]])   # 逆時針轉 90 度
F = Matrix([[1, 0],
            [0, -1]])                     # 對 x 軸鏡射

# 先鏡射再旋轉
print(simplify(R * F))
# Matrix([[0, 1], [1, 0]])：沿直線 y=x 鏡射

# 先旋轉再鏡射
print(simplify(F * R))
# Matrix([[0, -1], [-1, 0]]，不同！

from sympy import Matrix, symbols, simplify

c0, c1 = symbols('c0 c1')
A = Matrix([[1, 0],
            [1, 1],
            [1, 2]])
b = Matrix([0, 1, 1])

# normal equation: A^T A c = A^T b
ATA = A.T * A
ATb = A.T * b
print(ATA)   # Matrix([[3, 3], [3, 5]])
print(ATb)   # Matrix([[2], [3]])

c = ATA.inv() * ATb
print(simplify(c))   # Matrix([[1/4], [3/8]]) → y = 1/4 + (3/8)x

import matplotlib.pyplot as plt
import numpy as np

xs = np.array([0, 1, 2], dtype=float)
ys = np.array([0, 1, 1], dtype=float)
c0v, c1v = 0.25, 0.375

plt.plot(xs, ys, 'ro', label='資料點')
t = np.linspace(-0.5, 2.5, 100)
plt.plot(t, c0v + c1v * t, 'b-', label=f"擬合線 y = {c0v} + {c1v}x")
for xi, yi in zip(xs, ys):
    plt.plot([xi, xi], [yi, c0v + c1v*xi], 'g--')   # 誤差線
plt.legend()
plt.grid(True)
plt.show()

# 誤差平方和
err = sum((c0v + c1v*x - y)**2 for x, y in zip(xs, ys))
print(err)   # 0.125

# 對照：用 numpy 的最小平方法
coef = np.polyfit(xs, ys, 1)
print(coef)   # [0.375 0.25]：斜率與截距完全一致！

print(A.T * A .solve(A.T * b))   # Matrix([[1/4], [3/8]])

from sympy import Matrix, GramSchmidt

v1 = Matrix([1, 1, 0])
v2 = Matrix([1, 0, 1])
v3 = Matrix([0, 1, 1])

ortho = GramSchmidt([v1, v2, v3])
for v in ortho:
    print(v.T)
# Matrix([[1, 1, 0]])
# Matrix([[1/2, -1/2, 1]])
# Matrix([[-2/3, 2/3, 2/3]])

# 驗證兩兩垂直：內積皆為 0
print(ortho[0].dot(ortho[1]))   # 0
print(ortho[0].dot(ortho[2]))   # 0
print(ortho[1].dot(ortho[2]))   # 0

unit = GramSchmidt([v1, v2, v3], orthonormal=True)
for v in unit:
    print(v.dot(v))   # 1, 1, 1：每個都是單位向量

import matplotlib.pyplot as plt
import numpy as np
from sympy import Matrix, Rational

def least_squares_line(points):
    """points: [(x, y), ...]，回傳 (c0, c1) 與誤差平方和（精確分數）"""
    A = Matrix([[1, Rational(x)] for x, _ in points])
    b = Matrix([Rational(y) for _, y in points])
    c = (A.T * A).solve(A.T * b)
    err = sum((c[0] + c[1]*Rational(x) - Rational(y))**2 for x, y in points)
    return c, err

# 自己造資料：故意不平共線
points = [(0, 1), (1, 2), (2, 2), (3, 4)]
(c0, c1), err = least_squares_line(points)
print(f"y = {c0} + {c1}x，誤差平方和 = {err}")
# y = 4/5 + 4x/5，誤差平方和 = 1/5

xs = [p[0] for p in points]
ys = [p[1] for p in points]
plt.plot(xs, ys, 'ro', label='資料點')
t = np.linspace(min(xs) - 0.5, max(xs) + 0.5, 100)
plt.plot(t, float(c0) + float(c1) * t, 'b-', label='最小平方擬合線')
plt.legend()
plt.grid(True)
plt.show()