跳轉到

5.4 應用：Fibonacci 與馬可夫鏈#

本節目標#

  • 把 Fibonacci 遞迴改寫成矩陣冪次，用特徵值導出 Binet 公式——親眼看黃金比現身
  • 認識馬可夫鏈：轉移矩陣、機率向量、以及「穩態 = 特徵值 1 的特徵向量」
  • 用 np.linalg.eig 求穩態，再用迭代收斂實驗驗證
  • 走完「實際問題 → 矩陣 → 特徵值 → 預測長期行為」的完整套路

應用一：Fibonacci 數列 = 矩陣冪次#

Fibonacci 數列 \(0, 1, 1, 2, 3, 5, 8, 13, \dots\) 的規則人人都會背：

\[ F_{n+1} = F_n + F_{n-1} \]

遞迴要一步步算，第 100 項就得算 100 步。線性代數的魔法是： 把「兩個相鄰項」打包成一個向量，遞迴就變成矩陣乘法：

\[ \begin{bmatrix} F_{n+1} \\ F_n \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F_n \\ F_{n-1} \end{bmatrix}, \qquad \text{記 } A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \]

一路迭代下去：

\[ \begin{bmatrix} F_{n+1} \\ F_n \end{bmatrix} = A^n \begin{bmatrix} 1 \\ 0 \end{bmatrix} \]

**求第 100 項 = 求 \(A^{99}\) **——正好是 5.3 節剛練好的招式！

特徵值：黃金比現身#

\[ \det(A - \lambda I) = \lambda^2 - \lambda - 1 = 0 \quad\Longrightarrow\quad \lambda_{1,2} = \frac{1 \pm \sqrt{5}}{2} \]

看到 \(\frac{1+\sqrt{5}}{2} \approx 1.618\) 了嗎？**這就是黃金比 \(\varphi\) **！ \(F_{n+1} = F_n + F_{n-1}\) 這條樸素的遞迴規則裡，居然藏著古希臘的美學常數。 （這不是巧合： \(1 + \varphi = \varphi^2\) 正是遞迴關係本身的代數化身。）

Binet 公式推導#

由 5.3 節， \(A^n\) 的元素是 \(\varphi^n\) 與 \(\psi^n\) （ \(\psi = \frac{1-\sqrt5}{2} \approx -0.618\) ） 的線性組合；用初始條件 \(F_0 = 0\) 、 \(F_1 = 1\) 定出係數，得到閉式公式：

\[ F_n = \frac{\varphi^n - \psi^n}{\sqrt{5}} = \frac{1}{\sqrt{5}} \left[ \left(\frac{1+\sqrt5}{2}\right)^{n} - \left(\frac{1-\sqrt5}{2}\right)^{n} \right] \]

由於 \(|\psi| < 1\) ， \(n\) 大時 \(\psi^n \to 0\) ，所以：

\[ F_n \approx \frac{\varphi^n}{\sqrt{5}} \]

第 100 個 Fibonacci 數，不用迭代一步就能寫下。

import numpy as np

# 驗證三種演算法
n = 30
phi, psi = (1 + np.sqrt(5)) / 2, (1 - np.sqrt(5)) / 2

F_iter = int(round((phi**n - psi**n) / np.sqrt(5)))          # Binet 公式

A = np.array([[1.0, 1.0], [1.0, 0.0]])
F_power = int(np.linalg.matrix_power(A, n - 1)[0, 0])        # 矩陣冪次

vals, P = np.linalg.eig(A)
F_diag = (P @ np.diag(vals ** (n - 1)) @ np.linalg.inv(P) @ np.array([1.0, 0.0]))[0]

print(F_iter, F_power, F_diag)   # 832040 832040 832040.0 三者一致
print("F_n / F_{n-1} ->", (phi**n / phi**(n - 1)))            # 1.618... 收斂到黃金比

應用二：馬可夫鏈與穩態#

常客問題：某城市每人每個月只喝咖啡或茶，且習慣會變： 喝咖啡的人下個月有 \(70\%\) 繼續喝咖啡、 \(30\%\) 改喝茶； 喝茶的人下個月有 \(40\%\) 改喝咖啡、 \(60\%\) 繼續喝茶。 長期下來，市場佔比會變成什麼樣？

用狀態向量 \(\mathbf{x} = \begin{bmatrix} c \\ t \end{bmatrix}\) （ \(c + t = 1\) ，各佔比）描述人群，轉移規則寫成轉移矩陣：

\[ \mathbf{x}_{k+1} = M\mathbf{x}_k, \qquad M = \begin{bmatrix} 0.7 & 0.4 \\ 0.3 & 0.6 \end{bmatrix} \]

（注意每行加起來是 1——機率守恆。）這個「狀態只依賴上一步」的隨機模型 就是馬可夫鏈。

穩態 = 特徵值 1 的特徵向量#

「市場飽和、佔比不再變化」的數學描述是 \(\mathbf{x}^* = M\mathbf{x}^*\) —— 這正是特徵方程 \(M\mathbf{v} = \lambda\mathbf{v}\) 在 \(\lambda = 1\) 的特例！ 馬可夫轉移矩陣必定有特徵值 \(1\) （列隨機矩陣的性質），所以穩態一定存在：

\[ M\mathbf{x}^* = \mathbf{x}^* \quad\Longrightarrow\quad \begin{cases} 0.7c + 0.4t = c \\ 0.3c + 0.6t = t \end{cases} \ \Rightarrow\ 3c = 4t,\ \ c + t = 1 \quad\Longrightarrow\quad \mathbf{x}^* = \begin{bmatrix} 0.5714 \\ 0.4286 \end{bmatrix} \]

用 np.linalg.eig 求穩態 + 迭代驗證#

import numpy as np
import linalg_viz as viz

M = np.array([[0.7, 0.4],
              [0.3, 0.6]])

# --- 方法一：特徵向量法 ---
vals, vecs = np.linalg.eig(M)
print("特徵值:", vals)   # [ 1.  0.3]  ← 必有特徵值 1！
idx = np.argmin(np.abs(vals - 1))          # 找出最接近 1 的特徵值
x_star = np.real(vecs[:, idx])
x_star = x_star / x_star.sum()             # 標準化：機率和為 1
print("穩態:", np.round(x_star, 4))        # [0.5714 0.4286]

# --- 方法二：迭代驗證（模擬真實世界的收斂） ---
x = np.array([0.0, 1.0])                   # 假設一開始人人喝茶
history = [x.copy()]
for k in range(20):
    x = M @ x
    history.append(x.copy())

print("迭代 20 步後:", np.round(x, 4))     # 收斂到同一個穩態！

H = np.array(history)
viz.plot_lines(H[:, 0], H[:, 1], labels=["咖啡佔比", "茶佔比"])
viz.axhline(x_star[0], linestyle="--")     # 穩態畫成虛線
viz.title("馬可夫鏈收斂到穩態")
viz.show()

圖上兩條曲線一開始劇烈波動，之後貼著虛線走——不管初始狀態是 「人人喝茶」還是「各半」，長期都收斂到 \(\begin{bmatrix} 0.5714 \\ 0.4286 \end{bmatrix}\) 。 這就是「穩態由矩陣本身決定、與起點無關」的視覺證據。

為什麼一定收斂？特徵值給答案#

\(M\) 的特徵值是 \(1\) 和 \(0.3\) 。把狀態分解到特徵基底上：

\[ \mathbf{x}_k = M^k \mathbf{x}_0 = PD^kP^{-1}\mathbf{x}_0 = c_1 \cdot 1^k \mathbf{v}_1 + c_2 \cdot 0.3^k \mathbf{v}_2 \]

\(k \to \infty\) 時 \(0.3^k \to 0\) ，只剩特徵值 1 的分量—— 這就是收斂到穩態的完整解釋，而且 \(|0.3|\) 越小收斂越快。 （回頭看 5.3 的 \(A^n = PD^nP^{-1}\) ：應用的閉環在此合攏。）

換個模型：三態天氣#

import numpy as np

# 晴天/多雲/雨天，行 = 明天的條件機率
T = np.array([[0.6, 0.3, 0.1],
              [0.3, 0.4, 0.3],
              [0.1, 0.3, 0.6]])

vals, vecs = np.linalg.eig(T)
idx = np.argmin(np.abs(vals - 1))
steady = np.real(vecs[:, idx]); steady = steady / steady.sum()
print("穩態天氣分佈:", np.round(steady, 3))
# 約 [0.333 0.333 0.333]：長期三天氣均分

x = np.array([1.0, 0.0, 0.0])   # 今天大晴天
for k in range(30):
    x = T @ x
print("迭代 30 步:", np.round(x, 3))   # 同一個穩態

PageRank（第 8 章）就是同一套套路放大到數十億網頁： 「網頁重要性」是超連結轉移矩陣在特徵值 1 上的特徵向量。

小結與練習#

本節重點#

  1. 遞迴關係可打包成向量 + 矩陣： \(\mathbf{u}_{n+1} = A\mathbf{u}_n\) ，第 \(n\) 項 = \(A^n\) 作用在初值上
  2. Fibonacci 的特徵值是黃金比 \(\varphi = \frac{1+\sqrt5}{2}\) ，由此得 Binet 公式 \(F_n \approx \varphi^n / \sqrt5\)
  3. 馬可夫鏈：轉移矩陣每行和為 1；穩態 = 特徵值 1 的特徵向量（標準化和為 1）
  4. 收斂動力學由特徵值大小決定： \(|\lambda_2| < 1\) 保證收斂，且越小收斂越快
  5. 特徵值 1 的特徵向量法（直接解）與迭代法（模擬）必須給出相同答案——雙向驗證是好習慣

練習#

  1. 用 Binet 公式計算 \(F_{50}\) ，與 np.linalg.matrix_power 的結果比對。
  2. 常客問題中，若咖啡黨忠誠度升為 \(80\%\) （ \(M = \begin{bmatrix} 0.8 & 0.4 \\ 0.2 & 0.6 \end{bmatrix}\) ）， 手算新穩態，並用 np.linalg.eig 驗證。
  3. 畫出練習 2 的迭代收斂圖，比較新舊兩個矩陣哪個收斂更快 （提示：看 \(|\lambda_2|\) 的大小）。
  4. 三態天氣模型中，把雨天機率矩陣改成 \(\begin{bmatrix} 0.5 & 0.2 & 0.3 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.4 & 0.4 \end{bmatrix}\) ， 求新的穩態天氣分佈。
  5. 挑戰：證明列隨機矩陣（每行和為 1）必有特徵值 1 （提示：考慮 \(M^{\!\top}\mathbf{1} = \mathbf{1}\) ，即 \(M^{\!\top}\) 有特徵值 1， 再利用 \(M\) 與 \(M^{\!\top}\) 特徵值相同）。