劉應興的部落格

先前談過一個簡單的族群成長模型是假設族群大小 N(t) 是非隨機的，雖然可以允許淨成長率 a(t) 本身是一個隨機過程，又假設它對 t 可積，但它與真實情況總有些不符，例如這模型的 N(t) 不能限制在整數值；又如假設 a(t) 可積，若 a(t) 非隨機倒還好，但若認為 a(t) 應具有隨機性，可積的要求可說是一個重大的限制；而一方面認為 a(t) 是隨機的，另方面又把 N(t) 當做非隨機的，也有些怪異難諧。本文則將 N(t) 本身當做一個隨機過程，把它當做連續時間馬可夫過程 (continuous time Markov process) 的一個例子。

一個馬可夫鏈 Xn (Markov chain, 離散時間，n = 0, 1, 2, ...) 或馬可夫過程 Xt (連續時間，t ≧ 0)，基本上就是指隨機過程 Xt 滿足

P[Xt in A | Fs] = P[Xt in A | Xs]  for any (measurable) A and any s < t

簡單的（初級課程的）說法是：給定時間 s 及以前訊息，[Xt in A] 這事件發生的條件機率，只和時間 s 時的狀態 Xs 有關。數學上嚴謹講法，σ-體 Fs 代表時間 s 及以前的所有相關訊息，A 是狀態空間 (state space) S 的一個可測子集，所謂狀態空間就是所有隨機變數 Xt 的值域所在範圍，以生死過程而言就是 {0,1,2,...}。σ-體 Ft 可以定義為

Ft = σ({Xs, s ≦ t})

也就是由所有 Xs, s≦t 所生成的 σ-體，或說是總 σ-體中使 Xt, t≦s 可測的最小子 σ-體。前述馬可夫過程定義條件也稱為馬可夫性質 (Markov property)，在這條件中並沒有要求由 Xs 轉移到 Xt 的機率結構是與時間（特別是起點 s）無關的，雖然通常考慮的馬可夫過程其轉移機率 (transition probability) 是固定，與時間無關，稱為時間均勻的 (time homogeneous)。有時候像這裡談的，狀態空間是離散型的，不管時間是離散或連續都稱之為馬可夫鏈；而狀態空間是連續型的有兩種，一是離散時間的 Harris chain，另一是連續時間的，例如 Wiener 過程 (Brownian motion)。

我們考慮狀態空間為 S = {0, 1, 2, ...} 的連續時間馬可夫過程，Poisson 過程是一個例子，此過程在極短時間僅允許由 X(t) = n 轉移到 X(t+h) = n+1

P[X(t+h) = n+1 | X(t) = n] = λ(t)h + o(h)
P[X(t+h) = n | X(t) = n] = 1-λ(t)h + o(h)

不允許向左轉移（機率 0），X(t+h)-X(t) ≧ 2 的機率也可被忽略（機率 o(h)）。另有 Galton-Watson 過程和分支過程 (branching process) 倒是允許多個子代，但它們都考慮離散時間，也就是一代一代分得很明確。但在族群成長、等候問題中，族群大小或隊列長度卻比較像是依循生死過程 (birth and death process)

P[N(t+h) = n+1 | N(t) = n] = λ(n,t)h + o(h)
P[N(t+h) = n-1 | N(t) = n] = μ(n,t)h + o(h)
P[N(t+h) = n | N(t) = n] = 1-λ(n,t)h-μ(n,t)h + o(h)

當然，如果考慮雙胞胎、多胞胎的情況，或許由 N(t) = n 可以不限於轉移至 n-1, n, n+1 三狀態，而是可以轉移至更多狀態。

從上述生死過程的三個轉移機率設定式，令 P(n,t) = P[N(t)=n]，則得

P(n,t+h) = P(n,t)(1-λ(n,t)h-μ(n,t)h+o(h))
              + P(n-1,t)(λ(n-1,t)+o(h)) + P(n+1,t)(μ(n+1)h+o(h))

但 n = 0 時得

P(0,t+h) = P(0,t)(1-λ(0,t)h-μ(0,t)h+o(h)) + P(1,t)(μ(1)h+o(h))

計算 (P(n,t+h)-P(n,t))/h 並令 h → 0+，則得微分遞迴式（D 為對連續型變數微分）

      DP(n,t) = -(λ(n,t)+μ(n,t))P(n,t) + λ(n-1,t)P(n-1,t) + μ(n+1,t)P(n+1,t)
      DP(0.t) = -(λ(0,t)+μ(0,t))P(0,t) + μ(1,t)P(1,t)

在存 birth 的 Poisson process，DP(n,t) 只有 P(n,t) 和 P(n-1,t) 項，因此可遞迴解微分方程式；但此處 DP(n,t) 同時涉及 P(n,t) 與 P(n±1,t) ，這是與純出生過程不同的地方。

回到生死過程原來的假設，我們可以說

N(t+h) = N(t) + 1  機率 λ(N(t),t)h + o(h)
            = N(t) - 1   機率 μ(N(t),t)h + o(h)
            = N(t)        機率 1 - λ(N(t),t)h - μ(N(t),t)h + o(h)

如果 N(t) 代表一個地區的人口，在出生率、死亡率不隨時間而變的情況，我們可以假設

μ(n,t) = μn,   λ(n,t) = λn + θ

其中 θ 代表淨遷入，則

E[N(t+h) | N(t)] = N(t) + θh + N(t)(λ-μ)h + o(h)

因此，期望人口數滿足

M(t+h) = M(t) + θh + M(t)(λ-μ)h + o(h)

於是得微分方程

M'(t) = (λ-μ)M(t) + θ

設 M(0) = M0，解得

M(t) = M0 e^{(λ-μ)t} + (θ/(λ-μ))[e^{(λ-μ)t} -1]

如果沒有人口遷移，則簡化為指數成長模型；有人口遷移時，除了常數項，基本上也是指數成長。上列結果是假設 λ ≠ μ，若 λ = μ，則

M(t) = M0 + θt

因為自然增長為 0，而每年人口淨遷移數固定，總人口（期望值）成直線型變化。

基於上列人口期望數的結果，我們發現在出生死亡皆為固定比率時，除卻遷移因素，結果是指數型成長；很自然地我們猜測，如果將出生率做邏輯斯模型限制

λ(n,t) = λn(１-n/M*) + θ

是否可得邏輯斯成長的期望人口規模？依一般式 M(t) 的微分方程應是

M'(t) = E[λ(N(t),t) - μ(N(t),t)]

所以在出生率將隨 N(t) 增長而遞減時，我們得到的微分方程式和 V(t) = Var[N(t)] 有關：

M'(t) = (λ-μ)M(t) - (λ/M*)(M^2(t)+V(t)) + θ

但 V(t) 顯然與 λ(n,t) 及 μ(n,t) 的設定是有關的，想得到此模型的期望人口數成長函數仍有困難。

若 N(t) 是一等候線在線人數，來者按 Poisson 過程到達，到達率 λ；服務時間是指數時間，平均服務時間是 1/μ。如果只有一個服務窗口，這是 M/M/1 的等候系統；如果有 s 個服務窗口，就是 M/M/s 系統。所以 s = 1 時，

μ(n,t) = μ,   λ(n,t) = λ

對一般的 s,  λ(n,t) 不變，而

μ(n,t) = μn,  if n ≦ s;  = μs,  if n > s  

在線人數瞬間以 λ(n,t) h 的機率進來一個，又以 μ(n,t) h 的機率離去一個，典型的生死過程。

回到一般的生死過程，令 Tn 表示由狀態 n 轉移至狀態 n+1 所需的時間。首先，看 T0，

P[T0 > t+h] = P[T0 > t]P[T0 > t+h | T0 > t]

那麼 P[T0 > t+h | T0 > t] 等於多少？條件 T0 > t 蘊含 N(t) = 0，因此我們猜測

P[T0 > t+h | T0 > t] = P[N(t0+h) = 0 | N(t0) = 0]

可是這是事實嗎？由於馬可夫性質，給予 t 及以前的訊息都歸結到 t 時的狀態，因此

P[T0 > t+h | T0 > t] = P[T0 > t+h | N(t) = 0]

將事件 [N(t+h) = 0] 表示為

[N(t+h) = 0] = [T0 > t+h]∪[T0 ≦ t+h, N(t+h) = 0]

因此

P[T0 > t+h | N(t) = 0]  = P[N(t+h) = 0 | N(t) = 0] - P[T0 ≦ t+h, N(t+h) = 0 | N(t) = 0]

第二項表示在 (t,t+h) 中間某一點由狀態 0 轉移至狀態 1，又在另一點由狀態 1 轉移回狀態 0，所以最大也就是 O(h^2) 等級，故

        P[T0 > t+h | T0 > t] = P[N(t0+h) = 0 | N(t0) = 0] +o(h) = 1-λ(0,t)h+o(h)

如果以上沒弄錯，那麼

DP[T0 > t] = - λ(0,t) P[T0 > t]

結果得 P[T0>t] = e^{-∫_[0,t] λ(0,t') dt'}，若 λ(0,t) = λ0, 與 t 無關，則 T0 是一簡單指數分布變量，

E[T0] = 1/λ0

如果 λ(0,t) 不是常數，而是 λ0(t) 呢？

E[T0] = ∫_[0,∞) P[T0>t] dt = ∫_[0,∞) e^{-∫_[0.t] λ(0,u) du} dt

無 λ(0,t) 之明確算式，不能計算出 E[T0] 的具體結果。

現在考慮 n > 0 時 Tn 的期望值。假設 λ(n,t) = λ_n, μ(n,t) = μ_n。從狀態 n 下一瞬間可能移轉至狀態 n-1 或 n+1，也可能暫時停留，依 T0 的討論，我們可以說轉出去平均需要的時間是

由狀態 n 轉出（至 n-1 或 n+1）期望時間 = 1/(λ_n+μ_n)

我們可以認為由 n 轉移到 n+1 所需時間是一指數隨機變數 X，P[X>x] = e^(-λ_n x)；由 n 轉移至 n-1 所需時間是另一指數變量 Y，P[Y>y] = e^(-μ_n y)，X 與 Y 獨立。則 min(X,Y) 也是一指數變量，

P[min(X,Y) > z] = e^{-(λ_n+μ_n)z}

但 P[X<Y] = λ_n/(λ_n+μ_n)，P[X＞Y] =μ_n/(λ_n+μ_n)。而

E[Tn | X < Y] = 1/(λ_n+μ_n)
E[Tn | X > Y] = 1/(λ_n+μ_n) + E[T(n-1)] + E[Tn]

故合併條件期望值得

E[Tn] = 1/(λ_n+μ_n) + (μ_n/(λ_n+μ_n))(E[T(n-1)]+E[Tn])

或等價於

E[Tn] = 1/λ_n + (μ_n/λ_n)E[T(n-1)]

由 E[T0] = 1/λ0 開始，依上列遞迴式可得所有 E[Tn]：

        n    n       / n
E[Tn] = Σ    Π  μ_j / Π λ_j
       i=0 j=i+1   / j=i

式中 i = n 時分子連乘 j 是 n+1 起步，代表無任何一項，取值 1。

類似，用條件變異數方法及指數時間，可以導出 Tn 變異數的遞迴式。注意 Tn, n=0,1,2,... 是相互獨立的，因此我們也可計算由狀態 m 到 n > m 所需時間的期望值和變異數。至此讀者難免有個疑問：

你只談 n 到 n+1 的時間，那麼由 n 到 n-1 呢？

如 Tn, n > 0 的討論，令 Sn 代表由 n 初次轉移到 n-1 的時間，

E[Sn | X > Y] = 1/(λ_n+μ_n)
E[Sn | X > Y] = 1/(λ_n+μ_n) + E[S(n+1)] + E[Sn]

故得

E[Sn] = 1/(λ_n+μ_n) + (λ_n/(λ_n+μ_n))(E[S(n+1)]+E[Sn])

但不像 Tn 有個 E[T0] = 1/λ0 當起始值，這產生了困難。

在時間均勻假設下，令 P(t) 是連續時間馬可夫過程經過 t 時間的轉移機率「矩陣」（此處行列數並非必定有限，如生死過程），也就是說，其元素為

P(i,j,t) ≡ P[X(t) = j | X(0) = i]

則，對任意 s > 0, t > 0, 均得

P(i,j,s+t) = Σ_k P(i,k,t)P(k,j,s) = Σ_k P(i,k,s)P(k,j,t)

取 s = h = 很小的數，則

P(i,j,t+h) = Σ_k P(i,k,t)P(k,j,h) = Σ_k P(i,k,h)P(k,j,t)

在生死過程，中間一式得

　　  　DP(i,j,t) = -(λ(j)+μ(j))P(i,j,t) + λ(j-1)P(i,j-1,t) + μ(j+1)P(i,j+1,t)

這類似前面的狀態機率遞迴微分方程式

 DP(n,t) = -(λ(n)+μ(n))P(n,t) + λ(n-1）P(n-1,t) + μ(n+1)P(n+1,t)
 DP(0.t) = -(λ(0)+μ(0))P(0,t) + μ(1)P(1,t)

轉移機率的第二個等式，即第三式，在生死過程得

DP(i,j,t) = - (λ(i)+μ(i))P(i,j,t) + μ(i)P(i-1,j,t) + λ(i)P(i+1,j,t) 

在馬可夫過程（鏈）, 由 s 時間轉移矩陣和 t 時間轉移矩陣合成 s+t 時間轉移矩陣稱 Chapman-Kolmogorov 等式，導出之遞迴微分方程式，终點遞迴的（前一組）稱Kolmogorov 前向 (forward) 方程式；起點遞迴的（後一組）稱 Kolmogorov 後向 (backward) 方程式。

設 λ(k,t) = λk, μ(k,t) = μk, 即固定出生率、死亡率的情形，則 P(n,t) 的微分差分方程系統如下：

 DP(n,t) = -n(λ+μ)P(n,t) + (n-1)λP(n-1,t) + (n+1)μP(n+1,t)
 DP(0.t) =  μP(1,t)

令 g(z,t) 為時間 t 時機率分布 P(n,t) = P[N(t)=n], n = 0, 1, 2, ... 的機率母函數 E[z^N(t)]；給予 t = 0 時一個初始分布 P(n,0)，其機率母函數為 g(z,0)。將 P(n,t) 的微分式乘以 z^n 再對所有非負整數 n 加總，則由上列方程組可得

Dt g(z,t) = (λz-μ)(z-1) Dz g(z,t)

式中 Dt 表示對 t 做偏微，Dz 則是對 z 做偏微。由於 g(z,t) 是機率母函數，g(1,t) = 1 對所有 t；又加上初始機率分布之母函數 g(z,0)，如能解得上列偏微分方程，也就得到任意 t 時的族群大小機率分布 P(n,t)。

在 u = g(z,t) 這曲面的等高線上，其斜率

dz/dt = -(Dt g(z,t))/(Dz g(z,t)) = -(λz-μ)(z-1)

其解：

若 λ = μ, 則 t = 1/[λ(z-1)] + C
若 λ ≠ μ, 則 t = log[(z-1)/(z-μ/λ)] + C

式中 log 為自然對數。這表示 g(z,t) 可以表示為：

g(z,t) = ψ(t-1/[λ(z-1)]),              若 λ = μ;
g(z,t) = ψ(t-log|(z-1)/(z-μ/λ)|,   若 λ ≠ μ.

找出符合初始 g(z,0) 及 g(1,t) = 1 並且 g(z,t) 對所有正 z 值恆正的解，展開成 z 的冪級數，其係數即 P(n,t)。

如族群成長，可能是一個成長模式或衰退模式，前者 N(t) 將機率發散至無窮；後者最終將是族群滅絕。如果一個生死過程最終將導致均衡，意謂 P(n,t) 當 t→∞ 時將有個極限。注意 P(n,t) 是與時間 0 時狀態有關的；若看 t→∞ 的情形，起始狀態並不影響，P(i,n,t) 也將和 P(n,t) 有相同極限機率。這裡假設時間平穩，則依遞迴微分方程組，微分為 0 時得

λ(0)P(0) = μ(1)P(1)
(λ(n)+μ(n))P(n) = μ(n+1)P(n+1)+λ(n-1)P(n-1)

可導出

P(0) = lim_{t→∞} P(0,t)
     = 1/(1+Σ_k Π_{j=1~k}(λ(j-1)/μ(j)))
P(n) = P(0)Π_{j=1~n}(λ(j-1)/μ(j))

此平穩機率存在的充要條件是

Σ_k Π_{j=1~k}(λ(j-1)/μ(j)) < ∞

在 M/M/s 等侯系統，也就是

Σ λ^n/(sμ)^n < ∞

這條件相當於到達率 λ 比服務完成率 sμ 小，即 λ < sμ。在族群成長模型，固定出生率 λ，死亡率 μ 及固定淨遷移情況，則是 λ < μ，出生率小於死亡率；若出生率不小於死亡率，則由期望數 M(t) 也很容易看出族群是成長的（當 λ = μ 時則是線性的，θ > 0 成長，θ < 0 衰退，θ = 0  時平均數不動但變異數成直線上升。）

更詳細的結果可參閱隨機過程的書，或網路搜尋，如 Wiki。