基於馬爾科夫鏈的流網絡分析

從 集智百科
跳到: 導覽搜尋

一個平衡的流網絡可以自然地轉變成馬爾科夫鏈,而馬爾科夫鏈是數學家們研究得非常深入而透徹的數學對象,因此這為我們更好地分析流網絡提供了便利條件。

下面给出基本计算简表。定义:流矩阵F,对应的马尔科夫链矩阵M


 m_{ij}=\frac{f_{ij}}{\sum_{j=1}^{N+1}f_{ij}}
 ,   基础矩阵: 
 U=I+M+M^2+\cdot\cdot\cdot=(I-M)^{-1}
 

i到j的总流:
 t_{i,j}=T_0\frac{u_{0,i}u_{i,j}}{u_{i,i}}
 ,    i到j的首达流:
 \phi_{i,j}=T_0 \frac{u_{0,i}u_{i,j}}{u_{i,i}u_{j,j}}
 

其中T0为从源输入到整个网络的流量,也就是IS(参看流动网络)。根据这两个公式,可以计算任意网络任意节点对的总流和首达流。

i到j的平均时间:
 k_{i,j}=\frac{(MU^2)_{ij}}{(MU)_{ij}}
 ,     i到j的首达平均时间:
 l_{i,j}=\frac{u_{jj}(M_{-j}U_{-j}^2)_{ij}}{u_{ij}}=\frac{(MU^2)_{ij}}{u_{ij}}-\frac{(MU^2)_{jj}}{u_{jj}}
 

其中M_{-j}表示将M的第i行全置为0,U_{-j}=I+M_{-j}+M_{-j}^{2}+\cdot\cdot\cdot=(I-M)^{-1}


目錄

從流網絡到馬爾科夫鏈的轉換

圖1:一個示例流網絡

對於平衡流動來說,我們總可以將整個流系統看作是一個多粒子沿着網絡流動的系統。例如,對於圖1所示網絡,當某個處於1號節點的粒子往外流動的時候,它可能會跳到2號節點,也可能會跳到3號節點。可以想到,由於1到2號節點的流量比到3號的節點大,所以這個粒子更有可能跳到2號節點。我們不妨假設粒子跳轉到某個節點的概率是與這條邊上的流量成正比的。也就是說,粒子會以5/8=50/80的概率從1跳到2,而以3/8=30/80的概率從1到3。同樣的道理,當粒子跳入到任何一個節點後,它都會繼續沿着邊以概率跳轉。

我們可以寫出這個網絡所對應的流矩陣(參見流動網絡)如下:

示例網絡的流量矩陣

其中,矩陣外圍的數字為對節點的編碼。0對應源,6對應匯。陰影區域對應網絡的實線部分,為我們關注的節點主體。 整個平衡態的流網絡可以看作是一個馬爾可夫鏈,其中每個節點相當於是粒子所處的可能狀態,節點之間的跳轉可以看作是概率轉移。具體的,對於任何的處於平衡的流網絡F,我們定義矩陣M為概率轉移矩陣[1],其中M中的第i行第j列的元素就是粒子處於i狀態下轉移到j狀態的概率,記為:


m_{ij}=Pr\{to~~j|after~~visit~~i\}=\frac{f_{ij}}{\sum_{j=1}^{N+1}f_{ij}}

當粒子進入匯節點N+1,它就不會再往外跳了,也就是說匯所在的行全為0,因此,狀態N+1就相當於馬爾可夫鏈的吸收態。注意,根據矩陣M的定義,它自然滿足概率歸一化條件(不包括匯節點):


\sum_{j=1}^{N+1}m_{ij}=\frac{\sum_{j=1}^{N+1}f_{ij}}{\sum_{j=1}^{N+1}f_{ij}}=1,~~~~\forall i<N+1

針對圖1的示例網絡,轉化過來的馬爾科夫矩陣為:

Markovmatrix.png

其中陰影部分為我們主要關心的馬爾科夫轉移概率。

逆向馬爾可夫鏈

另外一個比較有趣的事實是,根據流矩陣F,我們還可以定義另外一個馬爾科夫鏈M',稱為逆向馬爾科夫鏈[2]。它刻畫的是一個逆向的因果過程,也就是說如果我在圖1中2號節點觀察到了一個新到的粒子,那麼它有可能從1號節點過來也可能從3號節點過來。顯然1到2的流量要比3到2的流量大得多,因此粒子從1跳過來的可能性會更大。那麼,既然整個網絡是平穩的,我們不難想到,在2號節點觀測到一個粒子的條件下,它是從1號節點跳過來的概率是5/6=50/60,它從3號節點跳過來的概率是1/6=10/60。

更一般地,我們定義逆向馬爾科夫鏈M',其中任意的元素m'ij定義為:觀察到粒子來到了節點i的條件下,它可能是從j節點流過來的概率。它按照下式計算:


m'_{ij}=Pr\{from~~j|observed~~particle~~on~~i\}=\frac{f_{ji}}{\sum_{j=0}^{N}f_{ji}}

同樣地道理,它也是歸一化的,即:


\sum_{j=0}^{N}m'_{ji}=1,~~~~\forall i>0

這個逆向馬爾可夫鏈又可以看作是逆向流矩陣FT(即原始流矩陣F的轉置)所對應的正向馬爾科夫矩陣,這是因為:


m'_{ij}=\frac{f_{ji}}{\sum_{j=0}^{N}f_{ji}}=\frac{F^T_{\{ij\}}}{\sum_{j=0}^{N}F^T_{\{ij\}}}=m_{ij}(F^T)

其中逆向流矩陣就相當於把原始流網絡所有的有向邊都掉轉方向,並保持所有的流量不變。由於網絡是平衡的,所以逆向的流網絡必然也是平衡的。

這似乎給我們描繪了一個有趣的圖景:在一個平衡的流網絡中,所有的粒子都順着流動從源到匯。而與此同時,如果我們從匯觀察到粒子,並不斷地詢問這個粒子是從哪個節點跳轉過來的,我們就得到了一個反向的信息流動,就彷彿是有假想的粒子從匯流到了源一樣。

路徑分析

下面,我們沿着粒子在流網絡中流動,分析粒子所能經過的所有可能路徑。讓我們以圖1所示網絡為例,考察從1節點到4節點的路徑。我們可以將所有從1到4的路徑分成2類:從1到4的首達路徑,以及包含了從4到4的循環路徑。而對於首達路徑,又可以分成2類:直接路徑、首達流循環路徑[3]

首達路徑

首先,我們可以列出從1到4的所有直接路徑,所謂的直接路徑是指粒子從1流到4,路徑上的所有節點只被訪問一次(沒有循環)的所有可能路徑,對於這個簡單的網絡來說,我們可以直接列出這些路徑是:


1\rightarrow 2\rightarrow 4


1\rightarrow 3\rightarrow 2\rightarrow 4

其次,從1到4的首達流循環路徑是指粒子從1流到4,中間允許任意次地循環重複經過某個節點,但是一旦到達4之後就不再循環重複的所有可能路徑。也就是說,4節點只能唯一地出現在路徑的尾端。對於圖1所示的網絡,我們知道1-->2-->5-->3-->2-->4是一條首達循環路徑,而1-->2(-->5-->3-->2)(-->5-->3-->2)-->4也是一條首達循環路徑,其中括號的部分被重複了兩次。當然,括號部分中重複3次、4次…… 都是不同的首達循環路徑。所以,我們可以把這一系列從5-->3-->2之間循環任意次的路徑簡記成


1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m\rightarrow 2\rightarrow 4

其中,指數的意思為括號部分的路徑字符可以重複m次,m為任意的整數(可以取無窮大),因此有無窮條從1到4的首達循環路徑。

值得指出的是,在這裡,我們相當於定義了兩個路徑的乘法,如果路徑P1=i1-->...-->i,路徑P2=j1,....,j,則P1乘以P2定義為:


P_1 \cdot P_2=i_1\rightarrow i_2 \rightarrow \cdot\cdot\cdot i \rightarrow j_1\rightarrow j_2\rightarrow \cdot\cdot\cdot j

可以驗證上述的路徑平方就是我們定義的路徑的乘法。

同樣的道理,還有一系列不同的從1到4的首達循環路徑:

綜上所述,我們實際上可以把所有的從1到4的首達路徑(包含了不循環的和循環的)統一表示成:


1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m\rightarrow 4,~~~~ 
1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2\rightarrow 4

其中m為大於等於0的任意整數。當m=0的時候,就是直接路徑。

循環路徑

第二大類路徑是那些到了4號節點又沿着網絡的流動重複到達4號節點的流動路徑,它們就是循環路徑,例如:1-->2-->4-->5-->3-->2-->4,又如:1-->3-->2-->5-->3-->4-->5-->3-->2-->4。利用我們上面引入的記號,我們可以把這兩個路徑進行推廣統一寫為:


[1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m\rightarrow 4]\cdot[(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^t


[1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2\rightarrow 4]\cdot [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^t

其中m,m',t都是大於等於0的任意整數。我們看到,上式兩種路徑就是在所有可能的首達路徑後面「後綴」了一個從4到4的首達路徑的任意次方。如果我們再定義路徑的加法為將兩個相加的路徑並列在一起(相當於集合的並集),則我們可以簡潔統一地把從1到4的所有可能路徑表達為:


\sum_{m,m',t}{[1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m \rightarrow 4]\cdot [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^t+ [1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2\rightarrow 4]\cdot [(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^{t}}

可以驗證路徑的加法和乘法同樣滿足結合律,所以上式又可以寫作:


\sum_{m,m',t}{\{[1\rightarrow 2 (\rightarrow 5\rightarrow 3\rightarrow 2)^m ] + [1\rightarrow 3 (\rightarrow 2\rightarrow 5\rightarrow 3)^m\rightarrow 2]\}\cdot [(\rightarrow 4\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^{t}}

(Template:EquationRef)

對於這種簡單的網絡來說,要表達出所有從1到4的路徑都已經如此複雜,當我們考慮更複雜的網絡的時候,機械地列出所有路徑就幾乎是不可能的了。所以,我們必須找出一套更合理的方法。

路徑矩陣

我們可以用矩陣來表示路徑,這樣,通過矩陣的乘積可以表示路徑的合成。首先,根據原始網絡的鄰接矩陣,我們定義一個矩陣P:


P_{ij}=\left\{\begin{array}{ll} i\rightarrow j & \mbox {if i,j connected}, \\
 0  &\mbox {else} \end{array}\right.

也就是,如果i和j相鄰,則在(i,j)的位置上放置i-->j,表示從i到j的路徑。該矩陣就是路徑矩陣。

例如,圖1所示網絡的路徑矩陣可以寫為:


P_1=\begin{pmatrix}
0       &0\rightarrow 1&0             &0             &0             &0             &0             & \\
0       &0             &1\rightarrow 2&1\rightarrow 3&0             &0             &0             & \\
0       &0             &0             &0             &2\rightarrow 4&2\rightarrow 5&2\rightarrow 6& \\
0       &0             &3\rightarrow 2&0             &0             &0             &3\rightarrow 6& \\
0       &0             &0             &0             &0             &4\rightarrow 5&4\rightarrow 6& \\
0       &0             &0             &5\rightarrow 3&0             &0             &5\rightarrow 6& \\
0       &0             &0             &0             &0             &0             &0             & \\
\end{pmatrix}

下面我們來考慮路徑矩陣的乘積。設路徑矩陣P和Q乘積(它們都是n*n的方陣)得到R,則


R_{ij}=\sum_{k}P_{ik}\cdot Q_{kj}

也就是說,路徑矩陣的乘積與普通矩陣一樣。注意等式右側的乘法對應路徑的乘法,也就是如果P_{ik}=i\rightarrow k, Q_{kj}=k\rightarrow j,則:P_{ik}\cdot Q_{kj}=i\rightarrow k\rightarrow j,即把兩個路徑拼接在一起。

這樣P1自乘一次就得到:


\begin{array}{ll}
P_1^2=\begin{pmatrix}
0       &0\rightarrow 1&0             &0             &0             &0             &0             & \\
0       &0             &1\rightarrow 2&1\rightarrow 3&0             &0             &0             & \\
0       &0             &0             &0             &2\rightarrow 4&2\rightarrow 5&2\rightarrow 6& \\
0       &0             &3\rightarrow 2&0             &0             &0             &3\rightarrow 6& \\
0       &0             &0             &0             &0             &4\rightarrow 5&4\rightarrow 6& \\
0       &0             &0             &5\rightarrow 3&0             &0             &5\rightarrow 6& \\
0       &0             &0             &0             &0             &0             &0             & \\
\end{pmatrix}\cdot
\begin{pmatrix}
0       &0\rightarrow 1&0             &0             &0             &0             &0             & \\
0       &0             &1\rightarrow 2&1\rightarrow 3&0             &0             &0             & \\
0       &0             &0             &0             &2\rightarrow 4&2\rightarrow 5&2\rightarrow 6& \\
0       &0             &3\rightarrow 2&0             &0             &0             &3\rightarrow 6& \\
0       &0             &0             &0             &0             &4\rightarrow 5&4\rightarrow 6& \\
0       &0             &0             &5\rightarrow 3&0             &0             &5\rightarrow 6& \\
0       &0             &0             &0             &0             &0             &0             & \\
\end{pmatrix}&\\
=

\begin{pmatrix}
0       &0             &0\rightarrow 1\rightarrow 2&0\rightarrow 1\rightarrow 3    &0             &0             &0         & \\
0       &0             &0     &1\rightarrow 3\rightarrow 2&1\rightarrow 2\rightarrow 4&1\rightarrow 2\rightarrow 5&(1\rightarrow 2\rightarrow 6)+(1\rightarrow 3\rightarrow 6)& \\
0       &0             &0     &2\rightarrow 5\rightarrow 3&0&2\rightarrow 4\rightarrow 5&(2\rightarrow 4\rightarrow 6)+(2\rightarrow 5\rightarrow 6)& \\
0       &0             &0     &0                          &3\rightarrow 2\rightarrow 4&3\rightarrow 2\rightarrow 5&3\rightarrow 2\rightarrow 6& \\
0       &0             &0     &4\rightarrow 5\rightarrow 3&0                          &0                          &4\rightarrow 5\rightarrow 6& \\
0       &0             &5\rightarrow 3\rightarrow 2&0     &0                          &0                          &5\rightarrow 3\rightarrow 6& \\
0       &0             &0     &0                          &0                          &0 \\
\end{pmatrix}
&\end{array}

其中i,j元素表示從i到j長度為2的所有路徑。同樣的道理,我們可以計算P1m,它的任意元素就表示從i到j長度為m的所有可能路徑。從這點來看,矩陣的乘積與路徑的合成有着天然的對應。

全路徑

接下來我們就要利用路徑矩陣這個工具計算從任意節點i到任意節點j的所有可能路徑。我們知道Pm的第i,j個元素包括了從i到j的長度為m的可能路徑,那麼要得到從i到j的所有可能路徑(無論路徑有多長),我們只要把所有的P1、P2……Pm……加起來就可以了。所以:


Paths\{i\rightarrow j\}=(\sum_{m=1}^{\infty}P^m)_{ij}

例如,對於圖1的網絡,1到4的所有可能路徑就是:


Paths\{1\rightarrow 4\}=(\sum_{m=1}^{\infty}P_1^m)_{1,4}

如果我們把P1代入到上式,可以驗證它跟我們在eq1中得到的全路徑的表達式是一樣的。

除此之外,我們也可以寫出所有從i到j的首達路徑(j在整個路徑中只能出現在結尾)為:


First\{i\rightarrow j\}=\frac{(\sum_{m=1}^{\infty}P^m)_{i,j}}{(\sum_{m=1}^{\infty}P^m)_{j,j}}

這裡,兩個路徑相除定義為乘法的逆運算。即,如果R\cdot Q=P,則定義\frac{P}{Q}=R。例如,若P=1\rightarrow 2\rightarrow 3\rightarrow 4, Q=3\rightarrow 4, 則\frac{P}{Q}=1\rightarrow 2

我們可以寫出圖1中所有從1到4的首達路徑為:


First\{1\rightarrow 4\}=\frac{(\sum_{m=1}^{\infty}P_1^m)_{1,4}}{(\sum_{m=1}^{\infty}P_1^m)_{4,4}}

根據eq1,我們知道所有的首達路徑拼接上(乘上)任意一條從4到4的路徑就構成了一條非首達路徑,那麼我們將所有的路徑除以一條從4到4的路徑就可以得到首達路徑了。也就是在eq1中將[(\rightarrow 5\rightarrow 3\rightarrow 2)^{m'+1}\rightarrow 4]^{t}這個因子除去。

流量分析

總流量

我們知道,每一條從i到j的路徑就代表了粒子從i到j的一種可能流動。當有大量的粒子在體系中流動的時候,每條路徑上都會分配一定的流量。下面,我們就要求出從i到j經過所有可能路徑的粒子總流量t_{i,j}。這個流量可以通過這樣的假想實驗來定義:假設原始的粒子沒有顏色,而粒子一旦經過節點i就被染成了紅色,並且一直保持不褪色。那麼,我從j節點處統計:每一時刻流入j節點的紅色粒子數有多少?這個數量就相當於從i到j的總流量:t_{i,j}

如果從i到j的所有路徑中不存在着循環,如下圖所示的樹狀結構:

Simpletreeill.PNG

我們很容易計算出從1到4的總流量為:


t_{1,4}=50\times \frac{3}{5}\times \frac{2}{3}=20=T_1 (M+M^2)_{1,4}

即1節點的總流量乘以經過一步概率轉移從1到4的概率,再加上總流量乘以經過2步概率轉移從1到4的概率。然而,當網絡中存在着循環結構的時候,計算稍微麻煩,但是我們不妨將上面通過馬爾科夫鏈計算的公式推廣到有環的結構上,也就是對M的m次方求和一直進行下去,直到無窮大。例如,對於如圖1所示的流網絡,從1到4的總流量為:


t_{1,4}=T_1 (\sum_{k=1}^{\infty}M^k)_{1,4}

這裡,T1為從源到1的流量,等於80。我們知道M1,4表示粒子直接從1到4的轉移概率,這樣T1M1,4就是粒子經過一步轉移從1到4的流量。同樣,(M^2)_{1,4}表示粒子經過兩步轉移從1到4的轉移概率,那麼T_1(M^2)_{1,4}就是粒子經過2步轉移從1到4的流量。……,我們將經過1步轉移、2步轉移、……從1到4的粒子流總量加起來自然就是從1到4的粒子流量。

基礎矩陣

然而,上式不方便計算,因為存在着矩陣的無窮級數和。注意到,M是一個行列式小於1的方陣,因此,我們可以用下列矩陣級數的求和公式:


(\sum_{m=1}^{\infty}M^m)_{1,4}=\frac{M}{I-M}

這裡I是與M同樣維度的單位陣。矩陣的無窮等比數列之和與數的無窮等比數列和的計算一模一樣。我們定義:



U=\frac{1}{I-M}=(I-M)^{-1}

基礎矩陣(Fundamental Matrix),這沿襲了投入產出分析的用法。(I-M)^{-1}表示矩陣I-M的逆。在這個例子中,U可以寫出來:


\left(
\begin{array}{ccccccc}
 1 & 1 & \frac{3}{4} & \frac{7}{16} & \frac{1}{4} & \frac{1}{2} & 1 \\
 0 & 1 & \frac{3}{4} & \frac{7}{16} & \frac{1}{4} & \frac{1}{2} & 1 \\
 0 & 0 & \frac{42}{41} & \frac{7}{82} & \frac{14}{41} & \frac{28}{41} & 1 \\
 0 & 0 & \frac{12}{41} & \frac{42}{41} & \frac{4}{41} & \frac{8}{41} & 1 \\
 0 & 0 & \frac{3}{164} & \frac{21}{328} & \frac{165}{164} & \frac{21}{41} & 1 \\
 0 & 0 & \frac{3}{82} & \frac{21}{164} & \frac{1}{82} & \frac{42}{41} & 1 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{array}
\right)

總流量

於是,將數字代入,從1到4的總流就是:


t_{1,4}=T_1 (MU)_{1,4}=80\times \frac{1}{5}=16

下面,我們再來計算從2到4的總流量。按照上面的類似計算,我們可以得到:


t_{2,4}=T_2 (MU)_{2,4}=60\times \frac{14}{41}=20.4878

然而,從圖1中,我們看到從2到4的路徑只有一條2-->4,所以從1到4的流量直觀地看應該是20,而我們上面計算的結果卻比20大。為什麼呢?仔細分析會發現,原因出在,對於2號節點的總流量T2包含着從2到4又循環回2的流量(從3到2的那部分流量)。而這部分流量又在U矩陣中被重複計算了。所以,要把它刨除掉。即只計算從0節點經過1流過來流量50,該流量就稱為從源0到2的首達流。如果所有粒子都是無色的,而只有經過2號節點的粒子被染成紅色,並且永不褪色,那麼,從0到2的首達流就是:


\phi_{0,2}=50+30\times \frac{2}{7}=\frac{410}{7}

不難看出,無色粒子可能從1直接到2,也可能經過3到2。至於首達流的更一般的計算將放到下一節討論。


t_{2,4}=\phi_{0,2} (MU)_{2,4}=\frac{410}{7}\times \frac{14}{41}=20

這樣的計算才可以得到正確的結果。因此,任意點i到j的總流量應該按下式計算:


t_{ij}=\phi_{0,i} (MU)_{ij}


根據上述公式,從i到i自己的總流量也就應該是:
t_{ii}=\phi_{0,i} (MU)_{ii}
然而,這個式子並不完全。因為,它並不包括流\phi_{0i}自身。習慣上,我們定義從i到i的總流應該包括\phi_{0i}。 因此,寫完全應該是:


t_{ii}=\phi_{0,i} ((MU)_{ii}+1)

注意到,根據U的定義,MU=U-I,所以


t_{ii}=\phi_{0,i} (U)_{ii}


最後,綜合上述所有討論,從i到j的總流量應該按照下式計算:


t_{ij}=\phi_{0,i}(U)_{ij}

即使對於i不等於j的情況上式也成立,因為U的非對角元與MU的非對角元相同。

首達流

上面的討論中曾引出了首達流的概念,下面我們就來討論一般地從i到j的首達流應該怎樣計算。首先,我們定義從i到j的首達流為所有從i出發首次到達j的粒子總流量。還是假設所有的粒子經過i節點後被染成紅色,並且同時假設粒子一旦經過j後,我們就把紅色塗掉,那麼每一時刻到達j的紅色粒子總數就是從i到j的首達流。這樣,所有經過i到達j的粒子就可以分成2類:'首達流循環流,它們的總和就是總流動。

例如,對於圖1來說,從1到4的總流動可以分為首達流和循環流:


t_{1,4}=\phi_{1,4}+\psi_{1,4}

我們已經知道t1,4的計算了,如果我們會計算ψ1,4的話,就很容易求出首達流了。所謂的從1到4的循環流一定是那些曾經至少已經經過一次4號節點的粒子所構成的流動。而我們知道,所有已經經過4號節點的粒子一定是那些首次到達了4號節點的粒子又重新從4出發,沿着所有可能路徑再流回到4節點的粒子。這些粒子的路徑就是如eq1式中所描述的,並且t不為0。同樣的,根據流量與路徑的關係,我們可以計算出:


\psi_{1,4}=\phi_{1,4}(\sum_{m=1}^{\infty}M^m)_{4,4}

也就是1到4的首達流再經過所有可能路徑回到4的總流量。這樣,我們就有:


t_{1,4}=\phi_{1,4}+\psi_{1,4}=\phi_{1,4}+\phi_{1,4}(\sum_{m=1}^{\infty}M^m)_{4,4}=\phi_{1,4}(1+(\sum_{m=1}^{\infty}M^m)_{4,4})=\phi_{1,4}u_{4,4}

其中,u4,4表示(U)4,4 這樣就可以求出1到4的首達流:


\phi_{1,4}=\frac{t_{1,4}}{u_{4,4}}

我們定義,當i=j的時候:


\phi_{i,i}=T_i

而我們已經知道了t1,4的表達式,因此1到4的首達流可以計算出,等於16*165/164=660/41。

上述方法適用於更一般的網絡,因此,我們可以得到計算任一點i到j首達流的正確計算公式為:


\phi_{i,j}=\frac{t_{i,j}}{u_{j,j}}

小結

總結來看,我們這裡定義了從i到j的兩種流動:i到j的總流和i到j的首達流。如果將經過i節點的粒子染色,則總流就是每個時刻j節點接收到的紅色粒子總數,而首達流就是所有那些第一次經過j節點的紅色粒子總數。它們的計算方式如下:

總流:


t_{i,j}=\phi_{0,i}u_{i,j}=\frac{t_{0,i}}{u_{i,i}}u_{i,j}=T_0\frac{u_{0,i}u_{i,j}}{u_{i,i}}

首達流:


\phi_{i,j}=\frac{t_{i,j}}{u_{j,j}}=T_0 \frac{u_{0,i}u_{i,j}}{u_{i,i}u_{j,j}}

其中T0為從源輸入到整個網絡的流量,也就是IS(參看流動網絡)。根據這兩個公式,可以計算任意網絡任意節點對的總流和首達流。

我們看到,無論是總流還是首達流,都要用到U矩陣,而


U=I+M+M^2+\cdot\cdot\cdot+M^{\infty}

相當於包含了所有從i到j的一步轉移概率、2步轉移概率、3步轉移概率……。但注意,uij並不是概率,因為u可能大於1。即便它不大於1,也要注意到M、M2, ……矩陣每一行都是歸一化的,所以把它們加起來就不再滿足歸一化條件(無論對於行還是列)。所以,uij並不表示從i到j的概率,而是表示i對j的總影響(Impact)。這是一個很含糊的概念,但目前也只能這樣來解釋U了。

而當uij乘以了一個流量φi,則所得的矩陣的每一個元素就有了明確的意義:表示從i到j的總流量。這個流量雖然從定義上看是一個瞬時t,i到j的粒子總數,但是根據計算式,它實際上是把t時刻的0步轉移、t-1時刻的一步轉移、t-2時刻的2步轉移,……,所有時間步的粒子數加起來。由於流網絡是平衡的,所以馬爾科夫鏈是平穩的,所以我們可以把瞬時空間上的流量轉變成所有歷史的求和,這一點是最讓人吃驚的事實。

平均時間(距離)

前面研究的是流量,下面我們研究與時間相關的主題。大量的粒子從網絡中流動,我們假設單位時間粒子沿着有向連邊跳轉。這樣,粒子需要m個時間步才能走完一條m長的路徑。下面,我們將研究從網絡上任意一點i到任意一點j的時間(距離)。

首先,由於網絡中存在着環路,因此,從i流向j的粒子有可能經歷無限長的時間步。為了避免這個問題,我們只考慮首達時間,也就是粒子從i出發第一次到達j的時間。但是,由於粒子可能沿不同的路徑到達j,而這些路徑的長度又不一樣,所以沒有確定的首達時間值。於是,我們只能用這些粒子的首達時間值來作為我們研究的目標,也就是平均首達時間。有趣的是,對於平衡的流網絡來說,平均首達時間的計算存在着簡潔而明確的表達式和計算方法。

任意點i到匯的平均時間

首先,我們來求從任意節點i(不包括匯)出發,到達匯N+1的粒子們的平均首達時間。由於到達N+1的粒子不再流向任何其它節點,因此這個平均首達時間其實就是粒子從i到達N+1的平均時間。下面,我們來計算從i到匯的平均時間。

我們不妨假設從任意節點i到達匯N+1的平均時間為li,那麼,我們可以列出如下等式:


l_i=1+\sum_{j=0}^{N}{m_{ij}l_j}, ~~~~\forall i\neq N+1

(Template:EquationRef)

下面來解釋這個式子。我們知道任何粒子要想跳到匯都有兩種可能性:直接跳到N+1,或者先跳到任意一個其他的非匯節點j,然後再跳到匯N+1。對於第一種情況,它的路徑長度確定性的是1,發生概率是m_{i,N+1}。而對於第二種情況,假設粒子跳到了j,再由j跳到N+1,那麼我們知道粒子從j跳到N+1的平均首達時間是lj,而粒子跳轉到j的概率是mi,j,因此粒子跳到任意的非匯節點,再跳到匯的平均路徑長度就是\sum_{j=0}^{N}{m_{ij}(1+l_j)},其中和號里的式子加1的原因是每條路經長度都要包含從i跳到j的這一步。於是,把這兩種可能性合起來就有:


l_i=m_{i,N+1}\cdot 1+\sum_{j=0}^{N}{m_{ij}(1+l_j)}=m_{i,N+1}+\sum_{j=0}^{N}m_{ij}+\sum_{j=0}^{N}m_{ij}l_j=1+\sum_{j=0}^{N}m_{ij}l_j, ~~~~\forall i\neq N+1

最後用到了馬爾科夫鏈的歸一化條件:\sum_{j=0}^{N+1}m_{ij}=1。注意到eq2實際上是一個線性方程組,未知數有N+1個,方程個數也有N+1個,可以直接求解。事實上,如果我們定義:向量L=(l_0,l_1,\cdot\cdot\cdot,l_N)^T,則eq2可以寫成矩陣的形式:



L=\mu+M_{-(N+1)} L


其中\mu=(1,1,\cdot\cdot\cdot,1)^T,1構成的N+1維列向量。M_{-(N+1)}表示原始馬爾科夫矩陣去掉最後一行和最後一列。則,最終L可以解出:



L=(I-M_{-(N+1)}^{-1})\cdot \mu=U_{-(N+1)}\cdot \mu=(\sum_{j=0}^{N}u_{ij})_{(N+1)\times 1}

(Template:EquationRef)


其中矩陣U-(N+1)為U矩陣去掉最後一行和最後一列。最後一個等式是根據馬爾科夫鏈刪除節點的性質中的定理1得到的。這樣L中的第i個元素就是i節點到匯的平均時間。

對於圖1所示網絡,我們可以計算出L=(63/16, 47/16, 175/82, 66/41, 525/328, 197/164),其中從源到匯的平均路首達時間是63/16,也就是平均每個粒子經過3.9375步跳出整個網絡。它剛好是\sum_{j=0}^{N}u_{0j}=\frac{C_0}{T_0}+1=\frac{TST}{IS}+1。這後兩個等式需要參考流網絡的異速標度律.

從源到任意點i的平均首達時間

我們無法簡單套用公式eq3來計算從源到任意點i的平均首達時間[2],這是因為如果將i節點去除,而非N+1,那麼整個網絡就有可能不連通,使得最後的求解出現問題。但是,我們仍然可以用類似的推理方法:我們假設0到任意i的平均首達時間為li。從源0到任意節點i的路徑總可以分成兩種情況,要麼從0直接到i,要麼從0到某個j再到i,於是假設0到j的距離為lj,則第二種情況下的平均首達時間就是:\sum_{j=1}^{N+1} m'_{ji}(l_j+1)。注意,在這裡,馬爾科夫轉移概率是逆向馬爾科夫轉移概率m'ji而不再是m了,這是因為我們已經確定性地知道了粒子肯定到達i的情況下,它可能從j節點過來的概率。於是,在這種情況下,概率歸一化的對象是給定在i觀察到粒子的情況下,從所有可能的中間節點過來的概率。所以,如果使用mji的話,則不能滿足概率歸一的條件,所以只能用m'ij。最後,再綜合第一種情況,就有:


l_i=m'_{0,i}+\sum_{j=1}^{N+1} m'_{ji}(l_j+1)=m'_{0,i}+\sum_{j=1}^{N+1} m'_{ji}+\sum_{j=1}^{N+1} m'_{ji}l_j=1+\sum_{j=1}^{N+1}m'_{ji}l_j,~~~~\forall i\neq 0

同樣可以把它寫作一個矩陣方程:


L=\mu+M'_{-0}L

於是,求解得:


L=(I-M'_{-0})^{-1}\mu=U'_{-0}\mu

注意,與eq3相比,我們把M換成了逆向馬爾科夫矩陣M'。也就是說,求從源0到任意點i的首達時間問題相當於是把整個流網絡倒轉過來,把源當成匯,把匯當成源,然後求任意節點i到匯的首達時間問題。

對於圖1所示網絡,可以計算出: L'=(1, 365/164, 193/82, 529/164, 285/82, 63/16) 注意,第一個元素對應的是0到1的平均首達時間,……,以此類推。那麼0到匯的平均首達時間是63/16,這與上面的計算結果是一致的。

任意節點i到任意節點j的平均首達時間

更一般地,我們應如何求出任意點i到j的平均首達時間呢?首先,我們需要明確這個平均首達時間到底是什麼意思。還是假設所有流經i點的粒子都被染上了紅色,並且粒子一經過j點就被抹掉紅色。這樣,j點接收的粒子就是首次到達j點的粒子。同時,假設每個粒子都有一個內置的秒錶。它一旦被染成紅色,秒錶就開始計時,每個時間步走動一次。於是我們在j點讀出粒子的計時,看它從i出發走了多少時間步過來。由於每個粒子的時間不一樣,於是,我們就可以求出一個均值。這就是我們要求的平均首達時間。

如何計算這個平均首達時間呢?上面介紹的方法只適用於源或匯,而不能簡單地拓展到非源或匯節點。下面,我們將介紹一種非常簡單實用的方法,我稱其為級數法。首先,根據平均首達時間的定義,我們知道從i到j的平均首達時間可以寫為:


l_{i,j}=\sum_{k=1}^{\infty}kp^{k}_{i,j}

這裡,pki,j是粒子經過k步從i首達j的概率。這個概率怎麼求呢?一種直觀的想法就是用(Mk)i,j來計算這個數值。但是,這是錯誤的,因為,pki,j這個概率是指所有那些從i到j首達的粒子中,剛好經過了k步的粒子比例。而(Mk)i,j表達的則是,在所有從i出發經過k步轉移後的粒子中,那些剛好到達節點j的粒子比例。這兩者顯然不一樣,因為從i出發的粒子經過k步後不見得都到達j。

讓我們只考察到達j的那些粒子,每一時刻,這個流量是固定的,有一些粒子是從i出發首達的粒子,而其中有一部分剛好經歷了k步。那麼,從i到j的首達粒子就可以用首達流來計算。剛好經過了k步的首達粒子如何計算呢?我們知道Φ0,i(Mk)i,j是所有經過k步轉移,從i到j的粒子數。但是這些粒子並不是首達粒子,有一些曾到過j又從j到j的粒子數。為了避免這種重複,我們不妨在M上做一些手腳:我們把M中j節點對應的行都設為0,即,定義:


(M_{-j})_{p,q}=\left\{\begin{array}{ll} m_{p,q} & \mbox {if } p\neq j, \\
 0 & \mbox {if } p = j\end{array}\right.

這樣,概率pki,j就可以這樣計算:


p^{k}_{i,j}=\frac{\phi_{0i}(M_{-j}^k)_{i,j}}{\phi_{i,j}}

於是平均首達時間就變成了求下列無窮級數和:


l_{i,j}=\sum_{k=0}^{\infty}k\frac{\phi_{0i}(M_{-j}^k)_{i,j}}{\phi_{i,j}}=\frac{\phi_{0i}(\sum_{k=1}^{\infty}kM_{-j}^k)_{ij}}{\phi_{i,j}}

那麼,下面的任務就是如何求解無窮級數\sum_{k=1}^{\infty}k(M_{-j})^k的和。為此,我們設此和為X,則在等式兩邊右乘以M-j得到:


XM_{-j}=\sum_{k=1}^{\infty}kM_{-j}^{k+1}=M_{-j}^2+2M_{-j}^3+\cdot\cdot\cdot

與等式X=M_{-j}+2M_{-j}^2+3M_{-j}^3+\cdot\cdot\cdot相減得到:


XM_{-j}-X=-(M_{-j}+M_{-j}^2+M_{-j}^3+\cdot\cdot\cdot)=-M_{-j}U_{-j}

於是,解出:


X=M_{-j}U_{-j}^2

其中,U_{-j}=\sum_{k=0}^{\infty}M_{-j}^k=(I-M_{-j})^{-1} 所以,i到j的平均首達時間就是:



l_{i,j}=\frac{\phi_{0i}(M_{-j}U_{-j}^2)_{ij}}{\phi_{i,j}}=\frac{\phi_{0i}(M_{-j}U_{-j}^2)_{ij}}{\frac{\phi_{0i}u_{ij}}{u_{jj}}}=\frac{u_{jj}(M_{-j}U_{-j}^2)_{ij}}{u_{ij}}

又根據馬爾科夫鏈刪除節點的性質,我們可以把上式簡化為:



l_{i,j}=\frac{u_{jj}(M_{-j}U_{-j}^2)_{ij}}{u_{ij}}=\frac{(MU^2)_{ij}}{u_{ij}}-\frac{(MU^2)_{jj}}{u_{jj}}


於是我們得到了從任意節點i到j的平均首達時間計算公式。將這個公式套用到圖1所示的網絡上,我們可以得到從任意點i到j的平均首達時間,這就構成了平均首達時間矩陣:

L=
\left(
\begin{array}{ccccccc}
 0 & 1 & 88/41 & 373/164 & 7218/2255 & 557/164 &
   63/16 \\
 \infty  & 0 & 47/41 & 209/164 & 4963/2255 &
   393/164 & 47/16 \\
 \infty  & \infty  & 0 & 9/4 & 58/55 & 5/4 & 175/82
   \\
 \infty  & \infty  & 1 & 0 & 113/55 & 9/4 & 66/41 \\
 \infty  & \infty  & 3 & 2 & 0 & 1 & 525/328 \\
 \infty  & \infty  & 2 & 1 & 168/55 & 0 & 197/164 \\
 \infty  & \infty  & \infty  & \infty  & \infty  & \infty  & 0
\end{array}
\right)

我們看到,這與前面計算的從0到所有點,以及從所有點到匯的距離是一致的,這進一步驗證了這個方法的正確性。

平均時間

根據上述級數法求解平均首達時間的啟發,我們還可以求出從任意點i到達任意點j的平均時間而非平均首達時間。這裡,平均時間按照如下方式定義:所有經過i的節點染成紅色,並且永不退色,而且,經過i的粒子就開啟一個計時器,則從j點處統計紅色粒子,計算它們計時器所經歷的時間。

有趣的是,按照上述方法,這個平均時間比平均首達時間更容易計算。還是按照平均時間的定義:


k_{i,j}=\sum_{k=1}^{\infty}kp^{k}_{i,j}

這裡p^{k}_{i,j}的計算不用再考慮首達問題,於是

p^{k}_{i,j}=\frac{\phi_{0i}(M^k)_{i,j}}{T'_{i,j}}

這裡,T'_{i,j}應為從i到j的總流。理應代入前面得到的關於T_{i,j}的表達式。這大體上正確,只有在i=j的時候不正確。因為在T_{i,i}中包含了從0到i的首達流\phi_{0,i},而在平均時間的計算中,這部分流量不應包含,故而T'_{i,i}=T_{i,i}-1。所以:


T'_{i,j}=\phi_{0,i}(MU)_{ij}

從而時間:


p_{i,j}^k=\frac{(M^k)_{i,j}}{(MU)_{ij}}

所以,平均時間的計算公式為:



k_{i,j}=\frac{(MU^2)_{ij}}{(MU)_{ij}}


根據這個公式,計算圖1所示網絡,得到彼此平均時間的矩陣為:


K=
\left(
\begin{array}{ccccccc}
 \infty & 1 & 365/164 & 193/82 & 529/164 & 285/82 &
   63/16 \\
 \infty  & \infty & 201/164 & 111/82 & 365/164 & 203/82
   & 47/16 \\
 \infty  & \infty  & 273/82 & 191/82 & 177/164 &
   109/82 & 175/82 \\
 \infty  & \infty  & 177/164 & 273/82 & 341/164 &
   191/82 & 66/41 \\
 \infty  & \infty  & 505/164 & 341/164 & 669/164 &
   177/164 & 525/328 \\
 \infty  & \infty  & 341/164 & 177/164 & 505/164 &
   273/82 & 197/164 \\
 \infty  & \infty  & \infty  & \infty  & \infty  & \infty  & \infty
\end{array}
\right)

有趣地是,對角線元素表達的是節點i到i的平均回歸時間,它可以用於節點i的中心度刻畫。

k_{i,i}=(\infty,\infty,273/82,273/82,669/164,273/82,\infty)

該數值越小,說明這個節點在整個網絡中的流量循環也就越快。

另外,我們還可以通過比較平均時間矩陣和平均首達時間矩陣,會發現任意兩點之間的平均時間都比平均首達時間要大一些,通過做這兩個矩陣的差,我們得到:


K-L=
\left(
\begin{array}{ccccccc}
        &  0 & 13/164 & 13/164 & 223/9020 & 13/164 & 0 \\
         &  & 13/164 & 13/164 & 223/9020 & 13/164 & 0 \\
   &   & 273/82 & 13/164 & 223/9020 &13/164 & 0 \\
   &   & 13/164 & 273/82 & 223/9020 &13/164 & 0 \\
   &   & 13/164 & 13/164 & 669/164 &13/164 & 0 \\
   &   & 13/164 & 13/164 & 223/9020 &273/82 & 0 \\
   &   &   &   &   &   & 
\end{array}
\right)

觀察這個矩陣,我們發現,所有有值的列除了對角線擁有相同的數值。這表明任何一點到達某一個點的平均時間與首達平均時間之差均為常數,並且出於同一環上的不同節點具有相同的數值,例如2,3,5。

參考文獻

歷史上,經濟學家最早給出了U矩陣(也稱為Leontief矩陣)的定義,請參照投入產出分析。後來生態學家Finn將這種分析方法引入到[1]了生態流網絡中用於分析能量流網絡。Patten、Ulanowicz等人將這套方法系統化,提出了Environ理論[4][5][6]。隨後,這些人又將這套方法進一步擴展以討論能量流的循環與存儲問題[7]以及能量流的節律等問題[8]

將流網絡建模成馬爾科夫鏈不僅僅可以討論流量、時間等問題,還可以討論如電導經過某點中轉的平均首達時間等問題。[9]


  1. 1.0 1.1 Finn, John T.. "Measures of Ecosystem Structure and Function Derived from Analysis of Flows". Journal of Theoretical Biology 56: 363-380.
  2. 2.0 2.1 Levine, Stephen. "Several Measures of Trophic Structure Applicable to Complex Food Webs". Journal of Theoretical Biology 83: 195-207.
  3. Higashi, Masahiko. "Network trophic dynamics: the modes of energy utilization in ecosystems". Ecological Modelling 66: 1-42.
  4. Patten, Bernard C. (1982). "Environs: Relativistic elementary particles for ecology". The American Naturalist 119: 179-219.
  5. Patten, Bernard C. (1985). "Energy Cycling in the Ecosystem". Ecological Modelling 28: 1-71.
  6. Ulanowicz, Roberte E.; Wolff, Wilfried F. (1990). "Ecosystem Flow Networks: Loaded Dice?". Mathematical Biosciences 103: 45-68.
  7. Patten, Bernard C.; Higashi, Thomas P.; Burns (1990). "Trophic Dynamics in Ecosystem Networks: Significance of Cycles and Storage". Ecological Modelling 51: 41-28.
  8. Higashi, Masahiko; Burns, Bernard C.; Patten (1993). "Network trophic dynamics: the tempo of energy movement and availability in ecosystems". Ecological Modelling 66: 43-64.
  9. Hilfer, A. (1988). "Probabilistic interpretation of the Einstein relation". Physical Review A 37: 578-581.
個人工具
名字空間
動作
導覽
工具箱