有限離散型二維條件分配相容性演算法之研究 - 政大學術集成

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

國 立 政 治 大 學 應 用 數 學 系

數 學 教 學 碩 士 在 職 專 班

碩士學位論文

有限離散型二維條件分配相容性

演算法之研究

On the algorithms for the compatibility

of bivariate finite conditional distributions

碩專班學生： 劉軒志

撰

指導教授：

姜志銘 博士

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

目錄

中 文 摘 要 - - - i

Abstract --- --- ---ii

1.

緒論---1

2.

理論背景---3

3.

演算法與程式設計---9

3.1 IBD 演算法---9

3.2 Rank One 演算法 ---17

3.3 對應矩陣 A 及 B 的聯合機率矩陣

JAB

之演算法 ---24

3.4 程式設計 ---31

4.

結論 ---42

參考文獻 ---48

附錄---49

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

附錄 1 程式使用說明--- 50

附錄 2 Main 主程式碼--- 57

附錄 3 檢查兩條件矩陣子程式碼---61

附錄 4 IBD 子程式碼---63

附錄 5 Rank One 子程式碼---66

附錄 6 隨機矩陣程式測試碼---69

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

中文摘要

給定兩個條件機率分配，判斷他們是否相容？是否有唯一的聯合機率分配？ 以及相容時，如何找出所有可能的聯合機率分配？是研究相容性相當重要的課

題。本文針對有限離散型二維條件機率分配，以 Arnold and Press(1989) 最先提

出的比值矩陣法，及由 Song , Li, Chen, Jiang, and Kuo (2010) 所提出的檢驗法為

架構，提出新演算法且利用此演算法來設計程式，使程式能判斷兩條件機率分配

是否相容，以及相容後可求出對應的所有聯合機率分配。本文亦依據新演算法並

應用 MATLAB 軟體設計程式，讓使用者可以很快地對上述三個問題得到答案。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

Abstract

When two conditional distributions are given, the following three important

questions are likely to be raised. Are they compatible? Is the corresponding joint

distribution unique if they are compatible? How do you find all the corresponding joint

distributions if they are compatible? In this thesis, basing on ratio matrix method given

first by Arnold and Press (1989), and on the method for checking compatibility

existence, for checking uniqueness, and for finding all possible joint distributions

provided by Song, Li, Chen, Jiang, and Kuo (2010), we provide a new algorithm to

answer these questions. Using this new algorithm, we also provide a MATLAB

computer program so that any user could get the answer quickly for the above three

questions.

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

1．緒論

1.1 前言

給定兩個離散條件機率分配，判斷他們是否相容？是否有唯一的聯合機率分 配？以及不唯一時，如何找出所有可能的聯合機率分配？是研究相容性相當重要 的課題。已有很多文獻提到解決這些問題的理論，但實務上相關的計算或演算法 的部分較少提及，本文研究的內容主要在研究演算法並設計程式。

1.2 研究動機與目的

本文針對二維離散型相容性問題為主體，以 Arnold and Press (1989)最先提出

的比值矩陣法，及由 Song, Li, Chen, Jiang, and Kuo (2010)所提的檢驗法為架構，研

究並提供新演算法，一方面讓後續研究者或實際使用者可以此演算法為基礎，針 對離散型相容性問題，以個人所擁有之軟體設計程式加以應用。另一方面本文亦 依據此演算法，利用 MATLAB 數學軟體，設計提供解決離散型相容性問題的程式， 讓一般實際使用者可直接使用本程式探討相容性問題。使用者輸入兩矩陣分別為 A 及 B 至本文程式，程式會自動判別此兩矩陣是否皆為條件矩陣；如果是，則再判 別是否相容；若相容，則可求出對應 A 及 B 的聯合機率矩陣。若唯一時，則可得 出唯一對應 A 及 B 的聯合機率矩陣，若不唯一，則可得出對應 A 及 B 的所有聯合 機率矩陣。本文研究主要在二維離散型相容性問題，希望給後續研究有關二維離 散型相容性問題者，可以提供程式應用，或以此演算法設計程式，在實際應用層 面對研究或探討相容性問題有所助益。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

1.3 架構

本文共分為四個章節。第一章緒論：包含前言、研究動機與目的及架構。第 二章理論背景：介紹有關二維離散型相容性文獻內容、理論背景包含定義、定理、 名詞解釋。第三章演算法與程式設計：首先介紹演算法，包括不可約化的塊狀對

角矩陣 (irreducible block diagonal matrix ，簡記為 IBD 矩陣) 演算法、

秩 1 (Rank One) 演算法、對應兩個條件矩陣 A 及 B 的聯合機率矩陣 (J_AB) 演算 法，最後再根據本文之演算法，利用 MATLAB 數學軟體所設計之程式(詳如附錄)，

以實例解釋演算法。第四章結論：首先說明並比較 Song, Li, Chen, Jiang, and Kuo

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

2．理論背景

本章所引用的定義與定理依據 Arnold and Press (1989)，Arnold, Castillo, and Sarabia (2004)和 Song , Li, Chen, Jiang, and Kuo (2010)之文獻內容。

假設X,Y為兩個有限離散隨機變數，令它們的值域分別為S_X {1, 2, , I}， } J , , 2 , 1 { SY   ，並令X,Y的聯合機率分配，邊際機率分配和條件機率分配分別為 XY P 、P 、_X P 、_Y PX|Y、PY|X。 定義 2 .1 令兩個 IJ矩陣 A[a_ij] 及 B [b_ij]，對所有的 ) , ( ji 



1, 2, , I}{1, 2, , J}，若 a_ij 0 且 1 1 



 I i ij a ，j ，則A矩陣為條件機率矩陣 (conditional matrix)。 若 b_ij 0 且 1 1 



 J j ij b ，i，則B矩陣亦為條件機率矩陣。 ﹟ 定義 2 .2 假設 A[a_ij]及 B[b_ij]為兩個 IJ條件機率矩陣，若對所有的 ) , ( ji 



1, 2, , I}{1, 2, , J}，存在聯合機率分配 P_XY使得它的條件機率 分配 P_X|Y(i| j)a_ij 且 P_Y|X(j|i)b_ij，則 A及 B為相容 (compatible)。﹟

Arnold and Press (1989) 以比值矩陣的觀念提出二維離散相容的充要條件。 首先定義比值矩陣如下： 定義 2 .3 令 A[aij]及 B[bij]為兩個 IJ條件機率矩陣， 若 A B N N  且           其它。 若 , 0 N N N ) j , i ( , b a c B A ij ij ij (2 .1)

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

其中 NA {(i,j)|a_ij 0}，NB {(i,j)|b_ij 0}， 未定義值以0表示， 則矩陣 C[c_ij]為 A對 B的比值矩陣。 註：若所有 cij 0，則稱 C是完整矩陣 (complete matrix)；否則稱 C是不 完整矩陣 (incomplete matrix)。 ﹟

Arnold and Press (1989)提出相容條件的定理如下：

定理 2 . 1 定義 cij如公式 (2 .1)，則兩個條件矩陣 A[aij]及 B[bij]相容 若且唯若

(1) NA NB N

(2) (i1,i2,j1,j2)恆有 ci₁j₁ci₂j₂ ci₁j₂ci₂j₁

其中 (i₁,j₁)N，(i₁,j₂)N，(i₂,j₁)N，(i₂,j₂)N。 ﹟

Song , Li, Chen, Jiang, and Kuo (2010) 舉出上述定理之反例，並顯示當C為不完整矩 陣時，此定理可能不成立；再綜合 Arnold, Castillo, and Sarabia (2004)定理 2， Perez-Villalta (2000) 定理 3，以ROPE觀念提出解決二維離散相容性的問題。

以上論述及反例詳見 Song , Li, Chen, Jiang, and Kuo (2010)。

以下將定義秩 1 正擴張矩陣 (ROPE)：

定義 2 .4 令 C=[c_ij]為IJ的比值矩陣，若存在另一個 IJ矩陣 C [c_ij]，使得所

有 c _ij 0，及c_ijc_ij，(i, j)N，以及rank(C)1，其中N定義於定理

2 .1，則稱矩陣 C是C的秩1正擴張矩陣 (rank one positive extension

matrix)，並簡記為 ROPE矩陣。 ﹟

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

註：rank(C)1 若且唯若 存在u(u1, , uI)，v(v1, , vJ)使得Cuv。(其中 u為 u 的轉置)

以下所提定理之完整証明詳見 Song , Li, Chen, Jiang, and Kuo (2010)。

定理 2 .2 令 C為A對 B的比值矩陣，則 A及 B是相容的 若且唯若 C有秩1正擴張矩陣。 ﹟ 當I及J值都很大時，或是比值矩陣中，未定義值所代表之 0 個數較多時，常常不 易判斷比值矩陣是否有ROPE矩陣，此時可將比值矩陣轉換成 IBD 矩陣來做判斷。 首先定義 IBD 矩陣： 定義 2 .5 若比值矩陣 C經過列互換或行互換後，可轉換成以下型式：         2 1 T T 0 0 其中對角塊狀矩陣 T₁和 T₂外的元素皆為0，則稱矩陣 C為可約化 (reducible)；否則稱C為不可約化 (irreducible)。 ﹟ 定理 2 . 3 給定任何一個比值矩陣 C，則經過列互換或行互換後，可轉換成如下型式 的矩陣：                  M 2 1 T T T ) C ( T 0 0 0 0 0 0 0 0 0 0             其中塊狀對角矩陣 T1, T2, , TM (M 1)皆為不可約化矩陣，且 M T T T₁, ₂, , 外的元素皆為0。 ﹟

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

為方便起見，我們將定理 2.3 之結果以T(C)Diag(T1, , TM)表示，並稱為不可 約化塊狀對角矩陣 (IBD 矩陣)。當M1，C 本身為 IBD 矩陣。

註：Song , Li, Chen, Jiang, and Kuo (2010)指出在檢驗相容性時，只須檢驗對角線上的塊 狀矩陣即可。 定理 2 . 4 設 C為A對 B的比值矩陣，且T(C) Diag(T₁, , T_M)為C的IBD矩陣， 其中 M 1。則 A和 B是相容 若且唯若 m T 有 ROPE矩陣，m，1mM。 ﹟ 若 A 及 B 是相容的，則邊際機率分配P 、_X P 及聯合機率分配_Y P_XY可用下面兩個等 價的系求得： 系 1. 令 C為A對 B的比值矩陣，若C存在 ROPE矩陣 Cuv， ) , , (u₁  u_I   u ，v(v₁, , v_J)，則 X 的邊際機率函數，Y的邊際機率 函數及XY的聯合機率函數可分別以   u u i P i X() 、



   I i i ij Y u u b j P 1 ) ( 及   u u b j i P i ij XY( , ) 其中



   I i i u u 1 ，(i, j)N。 ﹟ 系 2. 令 C為A對 B的比值矩陣，若C存在 ROPE矩陣 Cuv， ) , , (u₁  u_I   u ，v(v₁, , v_J)，則 X 的邊際機率函數，Y的邊際機率 函數及XY的聯合機率函數可分別以 j Y v v j P ( ) 、



  J j j ij X v a v i P 1 ) ( 及 j ij XY v a v j i P (, ) 其中



  J j vj v 1 1 1 ，(i, j)N。 ﹟

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

定理 2 . 5 令 C為A對 B的比值矩陣，若 C的 ROPE矩陣可表為 C uv， ) , , (u1  uI   u ，v(v1, , vJ)，則





      I i J j j i v v u u 1 1 1 1 。 ﹟

Song , Li, Chen, Jiang, and Kuo (2010) 證明了每一個ROPE矩陣對應一個聯合機率 分配，同時也以 IBD 矩陣觀點來提供唯一性的充要條件，這兩個重要的結果敘述如下： 定理 2 . 6 令 C為A對 B的比值矩陣，集合 是C的所有 ROPE矩陣的集合，集合 是所有聯合機率分配滿足 A及B為其兩個條件機率矩陣的集合。 令H(C)PXY，其中任意C之ROPE矩陣 C uv，u(u1, , uI)， ) , , (v₁  v_J   v ， 且



  _I i i i ij XY u u b j i P 1 ) , ( ， 則映射 H：是雙射的 (bijective)。 ﹟ 定理 2 . 7 令 C為A對 B的比值矩陣，若 A及B是相容的，則下列敘述是等價的： (1) 對任何C的IBD矩陣T(C)Diag(T₁,,T_M)，恆有 M 1。 (2) C為不可約化矩陣。 (3) C恰有唯一的 ROPE矩陣。 (4) 對應A及B的聯合機率矩陣是唯一的。 ﹟ 若條件機率矩陣 A 及 B 相容且對應C之ROPE矩陣 C 非唯一時，如何求出所有可 能對應C的聯合機率函數？方法敘述在定理 2.8：

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

定理 2 . 8 若兩條件矩陣 A及 B是相容的，且C為A對 B的比值矩陣，透過 列變換基本矩陣E及行變換基本矩陣F

可得IBD矩陣 T(C) ECF Diag(T1, , TM)。

對任意 m，1mM ，令 Tm為T_m的 ROPE矩陣，且 Tm u_mv_m ； 對所有的 k_m 0， 令 u_k (k₁u₁, , k_Mu_M)，且 Eu_k (w_k1, , w_kI)； 最後對任意 (i, j)N，再令



  _I i i i ij w w b j i P 1 ) , ( k k k ， 則



P_k(i, j)|k(k₁, , k_M), k_m 0, 1mM



是以 A及B為條件機率矩 陣的所有聯合機率矩陣的集合。 ﹟ 註： 由定理 2.8，若再定義 _           M M 1 1 k , , k v v v_k  ，                      F k , , k _M M 1 1 v v tk  ，因為u_kv_k 是 ECF ) C ( T  的ROPE矩陣，所以 Eu_kv_kFw_kt_k 是C的ROPE矩陣，其中 k k u w E ，則對應C所有可能之ROPE矩陣的集合為



w_kt_k |k(k₁, , k_M) , k_m 0 , 1mM



。 註：若k k，其中 0，則由定理 2.8，對任意(i,j)N，可推得P_k(i,j)P_k(i,j) ，故所有聯合機率矩陣的集合可再表示為       _

_

_{   }  M m 1 , 0 k , 1 k , ) k , , k ( | ) j , i ( P m M 1 m m M 1  k k 。 且對應C所有可能之ROPE矩陣的集合可再表示為       _{     }



 M m 1 , 0 k , 1 k , ) k , , k ( | _m M 1 m m M 1  k t wk k 。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

3．演算法與程式設計

3.1 IB D 演算法

給定一個IJ矩陣C=[cij](i,j)



1, 2, , I

 

 1, 2, , J



，C為C矩陣的轉 置。 步驟 0 令集合S



1, 2, , I



，R 1 m 步驟 1 S min =S的最小值 令C 為_m C之第minS列的子矩陣，並令E_m 



minS



(E 將表_m C 在原矩陣_m C中之列指標集合)、Fm 



j|cmj 0



(F 將表m C 各行中有非 0 元素之行m 指標集合) 步驟 2 令Sˆ ，且R RE_m 若R之元素個數 =I，則到步驟 4，否則到步驟 3 步驟 3 For k1, 2, , I，且kR 令B=C ．_m C_(k) (其中C_(k)表C之第k列的列矩陣) 若B 0，則將C(k)併入C 以更新m C ；將m k值加入E 集合以更新m E ；m

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

更新F 集合以表更新後之_m C 各行中有非 0 元素之行指標；將m k值加入 Sˆ 集合以更新 Sˆ 。 亦即C =_m           ) k ( m C C  、E =_m E_m {k}、F =_m F_m{j|c_kj 0}及SˆSˆ{k} End For 若Sˆ ，回步驟 2 否則令 (1)SS/E_m (亦即S為原矩陣C之I列指標中，不在E₁, E₂, , E_m中之 列指標集合) (2)mm1 回步驟 1 步驟 4 令Mm For m1, 2, , M 令 o m E 為將E 集合之元素由小到大排列之有序集合、_m F 為將_mo F 集合之_m 元素由小到大排列之有序集合 End For 步驟 5 將矩陣C中之列依 o M o m o 1, , E , , E E   之列指標順序重排，再將矩陣之行

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

依 o M o m o 1 , , F , , F F   之行指標順序重排即為C之IBD矩陣，記為T(C)。 令I₁, , I_m, , I_M表對應 o M o m o 1, , E , , E E   列指標個數之列個數，故 I I M 1 m m 



 ， M m 1, , J , , J J   表對應 o M o m o 1 , , F , , F F   行指標個數之行個數，故 J J M 1 m m 



 。 因此T(C)為具下列形式塊狀對角矩陣，T(C)=                 M T T T 0 0 0 0 0 0 0 0 0 0 2 1             ， 其中子矩陣T_m [t_ij]，





     m 1 k k 1 m 1 k k 1, , I I i  ，





     m 1 k k 1 m 1 k k 1, , J J j  。 例 1：給定一個 7 × 8 矩陣C                        0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 3 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 2 0 1 0 0 1 0 2 0 0 0 0 0 0 3 0 0 2 0 C 步驟 0 } 7 , , 2 , 1 { S    R 1 m

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

步驟 1



0 2 0 0 3 0 0 0



C1  ，E ={1}、1 F ={2, 5}1 步驟 2 令Sˆ ，且R RE₁={1} (R之元素個數< 7，到步驟 3) 步驟 3 For k1, 2, , 7 且 kR，亦即 k 2, 3, 4, 5, 6, 7 檢查B=C ．₁ C_(k) 當 k=2 , 3, 4 或 5 時，皆有B0 當k 6時，B0 此時C =_{1 }      0 0 1 1 0 0 0 0 0 0 0 3 0 0 2 0 ，更新後之E ={1, 6}、₁ F ={2, 5, 6}、₁ } 6 { Sˆ 當k 7時，B0 此時            0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 3 0 0 2 0 C1 ，更新後之E ={1 ,6 ,7}、1 F ={2, 5, 6}、1 } 7 , 6 { Sˆ

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

此時Sˆ ，回步驟 2 令Sˆ ，且R RE1={1 , 6, 7}(R之元素個數< 7，到步驟 3) 步驟 3 For k1, 2, , 7 且 kR，亦即 k 2, 3, 4, 5 檢查B=C ．₁ C_(k) 當k 2時，B0 此時              0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 2 0 0 0 0 0 0 3 0 0 2 0 C1 ，更新後之E1 {1 , 6, 7, 2 }、  1 F {2, 5, 6, 4}、Sˆ {2} 當k 3, 4, 5時，皆有B0 此時Sˆ ，回步驟 2 令Sˆ ，R RE₁={1 , 6, 7, 2} (R之元素個數< 7，到步驟 3 ) 步驟 3 For k1, 2, , 7 且 kR，亦即 k 3, 4, 5 檢查B=C ．₁ C_(k)

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

當k 3, 4, 5時，皆有B0 此時Sˆ ，令SS/E₁ {3, 4, 5}， m112 回步驟 1 步驟 1  S {3 ,4, 5}，R {1, 6, 7, 2} 2 C =



1 0 2 0 0 0 0 0



，Ε₂ {3}、F₂ {1, 3} 步驟 2 令Sˆ ，R RE2={1, 6, 7, 2, 3}(R之元素個數 < 7，到步驟 3 ) 步驟 3 For k1, 2, , 7 且 kR，亦即 k 4, 5 檢查B=C ．2 C(k) 當k 4時，B0 2 C = _      1 0 0 0 0 1 0 0 0 0 0 0 0 2 0 1 ，更新之後E₂ {3 ,4}、F2 {1, 3, 8}、Sˆ{4} 當k 5時，B0 此時Sˆ ，回步驟 2

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

令Sˆ ，且R RE2={1, 6, 7, 2, 3, 4} (R之元素個數 < 7，到步驟 3 ) 步驟 3 For k1, 2, , 7 且 kR，亦即 k 5 檢查B=C ．2 C(k) 當k 5時，B0 此時Sˆ ，令SS/E₂ {5} 3 1 2 m   回步驟 1 步驟 1 } 5 { S ，R {1, 6, 7, 2, 3, 4} 3 C =



0 0 0 0 0 0 3 0



，E₃ {5}、F₃ {7} 步驟 2 令Sˆ ，且RRE₃={1, 6, 7, 2, 3, 4, 5} (R之元素個數 = 7，到步驟 4 ) 步驟 4 3 M 將E ={1, 6, 7, 2},₁ E ={3, 4},₂ E ={5}集合之元素由小到大排成有序集合₃

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

o 1 E ={1, 2, 6, 7}, E ={3, 4},o₂ E₃o={5}、並將F ={2, 5, 6, 4},₁ F ={1, 3, 8},₂ 3 F ={7}集合之元素由小到大排成有序集合F ={2, 4, 5, 6},₁o F ={1, 3, 8},₂o o 3 F ={7}。 步驟 5 先將矩陣C中之列依 o 3 o 2 o 1, E , E E 之列指標順序重排，再將矩陣之行依 o 3 o 2 o 1 , F , F F 之行指標順序重排即為C之 IBD 矩陣，記為T(C)。結果如下： <E₁o, Eo₂, Eo₃>= <1, 2, 6, 7, 3, 4, 5 > <F₁o, F₂o, F₃o> = < 2, 4, 5, 6, 1, 3, 8, 7 >                        3 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 2 0 0 0 0 0 0 3 0 2 ) C ( T

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

3.2 R ank O ne 演算法

將比值矩陣C= B A 執行IBD程式轉換成IBD矩陣T(C)=                 M 2 1 T T T 0 0 0 0 0 0 0 0 0 0             ， 其中 M(= IBD Number) 為塊狀對角矩陣T 之個數。由定理 2. 4 可知 A 及 B 是相容的_m  T 有_m ROPE矩陣m，1mM。 對每一個矩陣T =[_m t_ij]，i



1, 2, , I_m},j{1, 2, , J_m}，以下演算法將決定T 是_m 否為秩 1 矩陣，若任一個T 不是秩 1 矩陣，則表示 A 及 B 為不相容，否則 A 及 B 為相m 容，同時可得出ROPE矩陣T =m u_mv_m。 步驟 0 令u =_m                 m I i 1 u u u   為T 矩陣的第 1 行_m 步驟 1 令U[u_ij]=u_m

．

m J 1 ，其中向量 m J 1 = [1, 1, , 1]內元素 1 的個數為T 矩_m 陣之行個數J 。_m

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

再令E[e_ij] 為將T 矩陣之元素 ÷_m U矩陣對應之元素 亦即e_ij =         0 u if 0 0 u if u t ij ij ij ij



1, 2, , I },j {1, 2, , J } i  _m   _m  For j1, 2, , Jm 令M_j= max{e_ij 0|1iI_m}(M_j表 E 矩陣第j行中之非 0 元素最大值者) j m = min{e_ij 0|1iI_m}(m_j表 E 矩陣第j行中之非 0 元素最小值者) 若 Cj = 10 j j j 10 m m M _   則停止，表示T不是秩 1 矩陣。 否則令 有值 不存在 j j j j M M if if M 0 v       End For 令v_m= [v ,₁ , v_j, , m J v ]，若向量u 及_m v_m 之元素皆不為 0，亦即 0 ui  ,1iIm 且 vj 0,1 jJm， 則到步驟 3 步驟 2 令V [v_ij]= m I 1 ．v_m，其中向量 m I 1 = [1, 1, , 1]內元素 1 的個數為T 矩_m 陣之列個數I

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

再令F[f_ij]為將T 矩陣之元素 ÷_m V矩陣對應之元素 亦即令fij =         0 v if 0 0 v if v t ij ij ij ij



1, 2, , I },j {1, 2, , J } i  m   m  For i1, 2, , I_m 令M = max{_i f_ij 0|1 jJ_m}(M 表 F 矩陣第_i i列中之非 0 元素最大值者) i m = min{fij 0|1 jJm}(m 表 F 矩陣第i i列中之非 0 元素最小值者) 若 Ri  10 i i i 10 m m M  _  則停止，表示T不秩 1 是矩陣。 否則令 有值 不存在 i i i i M M if if M 0 u       End For 令u_m= [u ,₁ , u ,_i , m I u ]，若向量u 及_m v 之元素皆不為 0，亦即_m 0 u_i  ,1iI_m 且 v_j 0,1 jJ_m， 則到步驟 3 否則回步驟 1

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

步驟 3 將向量u 及_m v_m相乘可得出ROPE矩陣T ，亦即m Tm = u_m v_m。 例 2：可填值為秩 1 矩陣 給定一個 4 × 4 的矩陣T              1 0 1 0 0 8 / 3 1 0 0 3 0 4 0 0 2 1 T 步驟 0 u             0 0 4 1 為T矩陣的第 1 行 步驟 1 令U= u ．[1, 1, 1, 1] =             0 0 0 0 0 0 0 0 4 4 4 4 1 1 1 1 ，其中 1 的個數為T矩陣之 行個數 = 4 令 E = U T

=

            0 0 0 0 0 0 0 0 0 4 / 3 0 1 0 0 2 1

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

For j1, 2, , 4 1 m M₁  ₁  ，M₂ m₂ 2， 4 3 m M₃  ₃  ，M₄,m₄不存在 得出 v= [1 2 3/4 0] 4 v = 0 到步驟 2 步驟 2 令 V=             1 1 1 1 ．v =             0 4 / 3 2 1 0 4 / 3 2 1 0 4 / 3 2 1 0 4 / 3 2 1 ，其中 1 的個數為T矩陣之列個數 = 4 令F= V T =             0 0 2 / 1 0 0 2 / 1 2 / 1 0 0 4 0 4 0 0 1 1 For i1, 2, , 4 1 m M₁  ₁  ，M₂ m₂ 4， 2 1 m M₃  ₃  ， 2 1 m M₄  ₄  得出更新後 u =             2 / 1 2 / 1 4 1 4 v = 0 回步驟 1

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

步驟 1 令U=u ．[1,1,1,1] =             2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 2 / 1 4 4 4 4 1 1 1 1 ，其中 1 的個數為T矩陣之行個 數 = 4 E

=

U T =             2 0 2 0 0 4 / 3 2 0 0 4 / 3 0 1 0 0 2 1 For j1, 2, , 4 1 m M1  1  ，M2 m2 2， 4 3 m M3  3  ，M4 m4 2 得出更新後 v = [1 2 3/4 2] 0 ui  ,1i4 且 vj 0,1 j4 到步驟 3 步驟 3 將向量u 、v相乘得出T=uv =             1 8 / 3 1 2 / 1 1 8 / 3 1 2 / 1 8 3 8 4 2 4 / 3 2 1 為T之ROPE矩陣。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

例 3：無法填值秩 1 為矩陣 給定一個 4 × 4 矩陣T T=             2 1 2 0 0 1 0 3 0 0 3 2 1 2 0 1 步驟 0 令 u             0 3 2 1 表T矩陣的第 1 行 步驟 1 令 U=u．[1, 1, 1, 1] =             0 0 0 0 3 3 3 3 2 2 2 2 1 1 1 1 ，其中 1 的個數為T矩陣之行個數 = 4 令 E

=

U T =             0 0 0 0 0 3 / 1 0 1 0 0 2 / 3 1 1 2 0 1 For j1, 2, , 4 1 m M₁  ₁  ， 2 3 m M₂  ₂  ，M₃ 2， 3 1 m₃  3 j ，C3 = 3 1 3 1 2 = 5 >1010 _{則停止，表示} T不秩 1 是矩陣。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

3.3 對應矩陣 A 及 B 的聯合機率矩陣

J_AB

之演算法

由定理 2.8 知給定 A 及 B 兩IJ之條件矩陣，比值矩陣C B A  透過IBD程式轉換成

IBD矩陣T(C)後，若每個不可約化塊狀對角矩陣T ，_m 1mM，經由 Rank One 程式 執行可得出ROPE矩陣T =m u_mv_m，則表示 A 及 B 是相容的，並可求出對應矩陣 A 及 B 的聯合機率矩陣 J_AB。IBD程式提供M個有序集合 o M o m o 1, , E , , E E   ，這M個有序集 合元素個數總合為I。 步驟 1 任取 k (k1, , km, , kM)，其中km 0 且



  M 1 m m 1 k 步驟 2 令 o E = <e₁, , e_i, , e_I> 為依照E₁o, , E_mo, , Eo_M的順序排序所得 到之有序集合，亦即 o E = <Eo₁, , Eo_m, , Eo_M> 再令I×I矩陣E=[e_ij]為將基本單位矩陣之列依有序集合 o E 之列指標順序 重排，亦即 i i ij e j e j if if 1 0 e         ，也就是在第i列中，除第e 行元素為 1 外，其餘為 0。_i

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

步驟 3 M m 1, , u , , u

u   為 Rank One 程式所記錄對應ROPE矩陣

M m 1, , T , , T T   之M個行向量 令 u_k 1 I M M 2 2 1 1 k k k               u u u  步驟 4 將u 還原回比值矩陣_k C之列指標排序得w ，亦即_k Eu_k w_k=               I 2 1 w w w k k k  ， E表E矩陣的轉置。 步驟 5 將 w 正規對角化，令_k                           k k k k k k k w w 0 0 0 w w 0 0 0 w w W I 2 1        ，其中



   I 1 i i w w_{k k} 。 步驟 6 對應 A 及 B 的聯合機率矩陣JAB=WkB。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

例 4：給定兩個56條件矩陣 A 及 B A=                 0 0 1 0 0 2 / 1 0 0 0 0 3 / 1 0 0 2 / 1 0 1 0 0 1 2 / 1 0 0 3 / 2 0 0 0 0 0 0 2 / 1 B=                 0 0 2 / 1 0 0 2 / 1 0 0 0 0 1 0 0 2 / 1 0 2 / 1 0 0 2 / 1 6 / 1 0 0 3 / 1 0 0 0 0 0 0 1 比值矩陣C=                  0 0 2 0 0 1 0 0 0 0 3 / 1 0 0 1 0 2 0 0 2 3 0 0 2 0 0 0 0 0 0 2 / 1 B A ，其中 0 代表未定義值。 由IBD程式將比值矩陣C轉換成不可約化塊狀對角矩陣T(C)。 如下所示： ) C ( T =                 0 0 0 3 / 1 0 0 0 1 2 0 0 0 2 3 0 2 0 0 0 0 0 0 2 1 0 0 0 0 0 2 / 1 (M=IBD Number = 2 ) < o2 o 1,E E > = <1, 5, 2, 3, 4 > (記錄在原矩陣C中之列指標有序集合) < o 2 o 1 ,F F > = < 1, 4, 2, 3, 5, 6 > (記錄在原矩陣C中之行指標有序集合)

由IBD程式所記錄之E 及m F 將m T(C)之塊狀矩陣T 選取出來，並執行 Rank One 程m 1

T T1

2

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

式可得出ROPE矩陣T 。m  1 T _      2 1 0 2 / 1 1 u = _      1 2 / 1 (第一個不可約化塊狀對角矩陣T 跑 Rank One 填值後，所記錄之行向量)1



1 2



1  v (第一個矩陣T 跑 Rank One 填值後，所記錄之列向量)1        2 1 1 2 / 1 T1 (T 填值後結果，其中₁ T =1 u₁v₁)            0 0 0 3 / 1 0 1 2 0 2 3 0 2 T₂            3 / 1 3 / 2 2 2 u (第二個不可約化塊狀對角矩陣T 跑 Rank One 填值後，所記錄之行向量)₂ 2 v=



1 3 3/2 1



(第二個矩陣T 跑 Rank One 填值後，所記錄之列向量)₂

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

           3 / 1 2 / 1 1 3 / 1 3 / 2 1 2 3 / 2 2 3 6 2 T2 (T 填值後結果，其中₂ T =2 u₂ v₂) T(C)=                 3 / 1 2 / 1 1 3 / 1 0 0 3 / 2 1 2 3 / 2 0 0 2 3 6 2 0 0 0 0 0 0 2 1 0 0 0 0 1 2 / 1 (填值後之T(C)) 可填值後，表示 A 及 B 是相容的，將T 矩陣之行向量m u 乘以_m k 倍後填入_m u_k k u = _      2 2 1 1 k k u u =                           2 2 2 1 1 k 3 1 k 3 2 k 2 k k 2 1 (其中km 0且



  2 1 m m 1 k ) 令 Eu_k w_k，亦即                 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1                           2 2 2 1 1 k 3 1 k 3 2 k 2 k k 2 1                          1 2 2 2 1 k k 3 1 k 3 2 k 2 k 2 1

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

                                                    2 1 1 2 1 2 2 1 2 2 1 2 2 1 1 k 3 k 2 3 k k 3 k 2 3 k 3 1 k 3 k 2 3 k 3 2 k 3 k 2 3 k 2 k 3 k 2 3 k 2 1 W k                                                                  0 0 k 3 k 2 3 k 2 1 0 0 k 3 k 2 3 k 2 1 0 0 0 0 k 3 k 2 3 k 3 1 0 0 k 3 k 2 3 k 3 2 2 1 0 k 3 k 2 3 k 3 2 2 1 0 0 k 3 k 2 3 k 2 2 1 k 3 k 2 3 k 2 6 1 0 0 k 3 k 2 3 k 2 3 1 0 0 0 0 0 0 k 3 k 2 3 k 2 1 J 2 1 1 2 1 1 2 1 2 2 1 2 2 1 2 2 1 2 2 1 2 2 1 2 2 1 1 AB 將w 正規對角化得出_k 將W 乘以 B 得出_k JAB

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

設定k =₁ k =₂ 2 1                  6 / 1 3 / 1 1 2 / 1 4 / 1 k u ， 令Eu_k w_k，亦即                 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1                 6 / 1 3 / 1 1 2 / 1 4 / 1 =                 2 1 6 1 3 1 1 4 1 將w 正規對角化得出_k

W

_k=                 9 / 2 0 0 0 0 0 27 / 2 0 0 0 0 0 27 / 4 0 0 0 0 0 9 / 4 0 0 0 0 0 9 / 1 AB J = W_kB =                 0 0 9 / 1 0 0 9 / 1 0 0 0 0 27 / 2 0 0 27 / 2 0 27 / 2 0 0 9 / 2 27 / 2 0 0 27 / 4 0 0 0 0 0 0 9 / 1

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

3.4 程式設計

本程式以 MATLAB 軟體設計，程式設計包含 4 個部分：(1) 輸入端檢查 (2) IBD 程 式 (3) Rank One 程式 (4) 對應矩陣 A 及 B 的聯合機率矩陣J_AB程式。

3.4.1 輸入端檢查：

設計判斷輸入端矩陣 A 及 B 是否符合條件矩陣型式，使用者在視窗介面下呼叫執 行檔 Main，並利用 Excel 輸入或讀入 A 及 B 矩陣資料，由程式判定是否為條件矩陣， 以決定是否執行 IBD、Rank One 程式，若不符合，則顯示不合之處，供使用者調整。 如下例所示：(黑框內為 MATLAB 程式視窗之擷取畫面，紅框為作者額外標記) 例 5：給定的矩陣A及B皆不是條件矩陣

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

例 6：給定的矩陣A及B為條件矩陣但不滿足 A B N N  ，故A及B不相容

3.4.2 IBD 程式：

透過輸入端檢查後，若輸入的矩陣為滿足 A B N N  的條件矩陣 A 及 B，本程式會 詢問是否執行 IBD 程式，若同意則產生比值矩陣C並執行 IBD 程式得出T(C)矩陣(程式 以 C_IBD 表示)，IBD 程式會記錄轉換成T(C)矩陣之列指標順序 o E 集合(程式以 E_set 表示)與行指標順序 o F 集合(程式以 F_set 表示)，並列出不可約化塊狀對角矩陣個數 M (程 式以 IBD_Number 表示 )。 如下例所示：

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

例 7：給定的矩陣A及B為條件矩陣且滿足 A B N N  ，故可執行IBD程式檔，並執行之。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

3.4.3

Rank One 程式

執行 IBD 程式後，延續 IBD 程式，會再詢問是否要執行 Rank One 程式，同意後

Rank One 程式會將T(C)矩陣中的不可約化塊狀對角矩陣T ，m 1mM，其中 M 為不

可約化塊狀對角矩陣個數(程式以 IBD_Number 表示)，依T ,1 T ,2 ,T 的順序逐步檢M

驗T 是否為秩 1 矩陣；若每個m T 矩陣皆為秩 1，則顯示對應m T 矩陣的m u 及m v 向量，m

並顯示填值後的T(C)矩陣T(C)(程式以 C_IBD_bar 表示)；若其中任一個不為秩 1，如第

m 個，則顯示T 不是秩 1 矩陣，即表 A 及 B 矩陣不相容，並中斷執行。_m 承例 7 (執行 IBD 程式後) 續跑 Rank One 程式

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

顯示填值後的T(C)矩陣T(C)(程式以C_IBD_bar表示)

3.4.4 對應矩陣 A 及 B 的聯合機率矩陣

JAB

程式

延續 Rank One 程式，若執行後每個不可約化的塊狀對角矩陣T 可填值為秩 1 矩_m 陣，則表示 A 及 B 矩陣是相容的，會再詢問是否求出對應 A 及 B 的聯合機率矩陣J_AB， 同意後，輸入 M 個正值比重k1, , km, , kM後，程式首先將比重正規化，再利用 Rank One 程式所記錄每個T 填值產生的_m u 及比重正規化值得出_m u (程式將以 uk 表_k 示)，最後又再利用 IBD 程式轉換成T(C)矩陣所記錄之列指標順序 o E 集合還原u 順序_k 得出w (程式將以 wk 表示)，最後將_k w 正規對角化後與 B 矩陣相乘即可求出_k J_AB (程 式以 J_AB 表示)，詳如J_AB演算法。 如下例所示：(本例輸入k1 1 及 k2 1)

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

承 承例例 77((執執行行RRaannkkOOnnee 程程式式後後)) 續續跑跑J_AB程程式式

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

下面例子顯示A及B滿足 A B N N  的條件矩陣但不相容 例 8： 給定的矩陣 A 及 B 為條件矩陣且滿足 A B N N  ，執行 IB D 程式，再執行 Rank One 程式後，發現不相容

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

3.4.5 程式測試之設計：

本節另外設計程式，以隨機產生聯合機率矩陣 J，再由 J 矩陣造出條件矩陣 A 及 B 予本程式執行，觀察執行後產生的J_AB矩陣與 J 矩陣誤差結果是否達合理範圍，在程式 設計上我們要求兩矩陣對應元素的誤差要小於 10 10 。 J 矩陣和J_AB矩陣的可能誤差主要發生在 Rank One 程式填值時。為檢驗 J 矩陣所造 的 C 矩陣是否為秩 1 矩陣，程式設計上，因為做矩陣相除的關係，可能產生數值誤差， 取值過程中，若兩個值的誤差除以兩者較小值不超過 10 10 ，則設定此兩個值相等，並取 最大值為共同值。最後所求之J_AB矩陣與原造 J 矩陣內相對應元素誤差之絕對值皆小於 10 10 ，本程式設計已符合需求。 如下例所示：(黑框內為 MATLAB 程式視窗之擷取畫面，箭頭與藍字為作者額外標記) 矩 矩陣陣大大小小以以11000000xx11000000為為例例 由 由JJ矩矩陣陣造造出出條條件件矩矩陣陣AA及及BB

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

圖為 Matlab 程式視窗下所顯示之條件矩陣 A 及 B

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

註： 本文末附錄為本程式所寫之使用說明與程式碼，供讀者參考。 隨機產生之聯合機率矩陣 J，IBD Number =1 執行程式後所得出之聯合機率矩陣 JAB 兩矩陣相對應元素誤差之絕對值 兩矩陣相對應元素誤差之絕對值加總後的結果

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

4．結論

4. 1 比較兩種演算法之差異

本文與 Song , Li, Chen, Jiang, and Kuo (2010) , Compatibility of finite discrete conditional

distributions.所發表之演算法比較差異如下：

在 IB D 部分：

兩種演算法皆針對非 0 元素，做列與行的互換重排，使得非 0 元素往對角線靠近，0

元素往非對角線移動，以形成塊狀矩陣。Song , Li, Chen, Jiang, and Kuo (2010)之演算法針

對矩陣元素，優先做列的移動，必要時再做行的移動，因此若第一列第一行之元素為 0， 將在第一行中非 0 之元素所對應的最小列經列互換成為新的第一列後，再以鎖定熱區建構 方式，將非 0 元素往對角線上移動，而本文對非 0 元素的移動方式則在鎖定第一列後，針 對第一列中的非 0 元素，是否在同一行其它位置有非 0 元素值，如果有，則將該非 0 元素 所在的列與其它列做列互換，並與第一列合併為新矩陣的前幾列，最後才將前幾列中有非 0 元素的行經行互換成為新矩陣的前幾行，兩者在移動的步驟上有所區別。

在 R ank O ne 部分：

兩種演算法皆鎖定各行與第一行的比值是否一致來填值，若不一致，則表示不為 Rank One 矩陣，若一致，則再進而鎖定各列與第一列的比值是否一致來填值，若不一致，則表

示不為 Rank One 矩陣；若一致，則確定為 Rank One 矩陣，其填值在觀念上是一致的，差

異的部分是在操作的方法上，本文演算法是以原矩陣第一行為參考行，令其它行皆與第一

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

依此比例所得之值當作參考列，再取其它列使與參考列相同，而得出新建構矩陣，並將原 矩陣除以新建構矩陣，觀察各列是否成比例，若成比例，則依此比例所得之值當作新參考 行，重複行之，一旦所求之參考行與參考列中元素皆完成填值動作後，將參考行與參考列

相乘即為 ROPE 矩陣，Song , Li, Chen, Jiang, and Kuo (2010)演算法則鎖定熱區，計算對應

熱區內各行、列中非 0 元素的比值是否一致，若比值不一致，則停止，且不為 Rank One 矩陣；若比值一致，則找出對應在同一行、列中的 0 元素，再以比值填入，以形成 ROPE 矩陣。

在聯合機率矩陣

JAB

部分：

兩演算法皆利用比值矩陣轉換成 IBD 矩陣所記錄的列指標順序集合 o E (即比值矩陣 C B A  轉換成 IBD 矩陣之列指標順序)，與 IBD 矩陣中，將每個塊狀對角矩陣T 可填值為_m 秩 1 所得到的u (每個塊狀對角矩陣_m T 可填值為 ROPE 矩陣_m T 之行向量)，來求出對應矩m 陣 A 及 B 的聯合機率矩陣J_AB。首先令k_m 0乘入對應的行向量u_m

(

其中



  1 m m 1 k )

，

得出 k u ，再利用 o E 還原u 順序得出_k w ，將_k w 正規對角化後得出_k W ，最後將_k W 與_k B矩陣相 乘而求出對應矩陣 A 及 B 的聯合機率矩陣JAB，在步驟上是大同小異。

4.2 程式設計之測試模擬結果

本文在 3.4.5 以模擬方式對程式做測試，由程式隨機造出聯合機率矩陣J，再由J矩陣 產生條件矩陣A及B，輸入此A及B矩陣予程式執行，得出對應 A 及 B 的聯合機率矩陣 AB J ，為了比較兩個矩陣是否一致，故檢驗兩個矩陣J與JAB元素的誤差是否極小，亦即是 否能在 10 10 以下，程式模擬的結果詳見如下之兩表： (E 表以 10 為底的次方，E-15 表 15 10 )

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

表 1：1000×1000 矩陣之模擬誤差 模擬矩陣大小：1000×1000 模擬次數 1 2 3 4 5 6 7 8 9 10 兩矩陣每個 元素值誤差 之絕對值加 總結果

8.33E-16 1.35E-15 2.78E-16 5.11E-16 1.03E-15 8.47E-16 2.04E-15 7.84E-16 1.05E-15 3.58E-16

1.06E-21 2.33E-21 1.06E-22 8.47E-22 1.06E-21 7.41E-22 5.82E-22 4.24E-22 2.12E-21 2.12E-22

1.69E-21 2.54E-21 4.24E-22 4.24E-22 2.65E-22 1.27E-21 1.06E-21 4.24E-22 7.41E-22 1.06E-22

1.59E-22 2.12E-21 2.12E-22 1.59E-22 8.47E-22 1.48E-21 6.35E-22 8.47E-22 7.41E-22 2.12E-22

1.27E-21 6.35E-22 2.12E-22 6.35E-22 7.41E-22 2.12E-21 2.33E-21 1.27E-21 1.69E-21 4.24E-22

1.06E-21 2.54E-21 4.24E-22 6.35E-22 1.91E-21 1.59E-22 3.81E-21 1.06E-21 6.35E-22 8.47E-22

1.32E-22 8.47E-22 1.06E-22 4.24E-22 8.47E-22 1.27E-21 8.47E-22 1.06E-21 2.54E-21 4.24E-22

4.24E-22 2.33E-21 4.24E-22 6.35E-22 9.53E-22 3.71E-22 2.33E-21 1.32E-22 2.54E-21 4.24E-22

1.06E-21 1.91E-21 4.24E-22 3.18E-22 7.41E-22 7.41E-22 3.18E-21 1.48E-21 8.47E-22 6.35E-22

1.27E-21 8.47E-22 2.12E-22 2.12E-22 1.27E-21 7.94E-23 5.82E-22 6.35E-22 6.35E-22 2.12E-22

比較兩矩陣 第一行前 10 個元素的誤 差之絕對值

1.69E-21 1.91E-21 4.24E-22 5.29E-23 1.06E-21 5.29E-22 2.96E-21 1.06E-21 2.38E-22 1.59E-22

表 2：700×800 矩陣之模擬誤差 模擬矩陣大小：700×800 模擬次數 1 2 3 4 5 6 7 8 9 10 兩矩陣每個 元素值誤差 之絕對值加 總結果

1.43E-15 4.66E-16 5.41E-16 1.34E-15 4.49E-16 5.73E-16 5.05E-16 3.15E-16 5.98E-16 6.03E-16

3.39E-21 1.69E-21 2.12E-21 3.39E-21 4.24E-22 4.24E-22 2.12E-21 2.12E-22 6.35E-22 1.27E-21 1.91E-21 8.47E-22 2.12E-21 1.91E-21 1.27E-21 1.27E-21 8.47E-22 4.24E-22 1.27E-21 8.47E-22 1.48E-21 1.69E-21 2.65E-22 2.54E-21 6.35E-22 1.06E-21 4.24E-22 2.65E-23 8.47E-22 7.94E-23 2.12E-21 6.35E-22 2.12E-21 4.24E-22 8.47E-22 2.54E-21 8.47E-22 2.12E-22 1.69E-21 8.47E-22 1.27E-21 8.47E-22 1.27E-21 4.24E-21 1.06E-22 1.69E-21 4.24E-22 3.18E-22 1.06E-21 3.18E-22 2.75E-21 1.27E-21 2.12E-22 3.71E-22 4.24E-22 1.59E-22 1.27E-21 1.27E-21 1.69E-21 8.47E-22 2.75E-21 8.47E-22 6.35E-22 4.24E-21 1.99E-23 1.06E-21 5.29E-22 8.47E-22 6.35E-22 1.06E-21 1.06E-21 6.35E-22 1.27E-21 1.46E-22 6.35E-22 1.06E-22 1.69E-21 1.27E-21 5.29E-22 1.69E-21 2.54E-21 3.18E-22 2.12E-21 1.69E-21 3.18E-22 2.12E-22 1.06E-21 1.27E-21 2.12E-21 6.35E-22

比較兩矩陣 第一行前 10 個元素的誤 差之絕對值

3.39E-21 2.12E-22 2.12E-21 1.59E-22 3.18E-22 1.69E-21 3.18E-22 1.59E-22 2.54E-21 1.69E-21 由上兩表可發現，針對方陣或不是方陣的矩陣做模擬的結果，誤差值皆在合理範圍內。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

4.3 程式應用之建議

本程式是利用 MATLAB 軟體所編寫，後續研究者或使用者可朝有關程式介面的相容 性、使用者輸入操作…等再加以修改設計，使程式能符合更多研究者與使用者的需求，增 進程式的應用效果，亦可利用其它適當的軟體，編寫新的程式，以適應不同的特殊環境。 本文研究範疇只在 2 維，亦可進一步推廣至 3 維、或更高維度，使其研究能有更多貢獻。

4.4 IB D 之第二種演算法的建議

有關IBD演算法，給定任意矩陣C經過基本列、行互換後得出IBD矩陣T(C)，亦可以 圖示連通(connect)的原理來做構想，步驟為針對任意IJ矩陣C，令SI {1, 2, , I}表矩 陣C之列指標集合，SJ {1, 2, , J}表矩陣C之行指標集合，找出矩陣C中之非 0 元素， 亦即令集合E{(r,c)0|rSI,cSJ}，觀察E中每個元素列行指標對應連通情形，將有 連 通 之 列 、 行 指 標 集 合 分 別 由 小 到 大 排 序 ， 可 得 E , F , E , F2o, o 2 o 1 o 1 等 ， 再 依  , F , F , , E , E₁o o_{2 1}o ₂o 重新排序，則可得出 IBD 矩陣T(C)。 令為T ，m 1mM，則IBD矩陣T(C)=                 M 2 1 T T T 0 0 0 0 0 0 0 0 0 0            

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

列指標 1 2 3 4 5 6 7 行指標 1 2 3 4 5 6 7 8 以本文第三章IBD演算法所舉之例 1 來說明： 矩陣                        0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 3 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 2 0 1 0 0 1 0 2 0 0 0 0 0 0 3 0 0 2 0 C } 7 , 6 , 5 , 4 , 3 , 2 , 1 { S_I  表矩陣C之列指標集合，S_J {1, 2, 3, 4, 5, 6, 7, 8}表矩陣C 之行指標集合，將矩陣C中之大於 0 元素找出，得出E={(1,2)、(1,5)、(2,4)、(2,6)、(3,1)、 (3,3)、(4,3)、(4,8)、(5,7)、(6,5)、(6,6)、(7,6)}。其中(1,2)表將列指標 1 與行指標 2 以直線 相連，其它坐標以此類推，觀察列行指標連通情形發現： 列指標集合{1,2,6,7}令為 o 1 E 與行指標集合{2,4,5,6}令為F 有連通，以1o 紅色線 表之。 列指標集合{3,4}令為 o 2 E 與行指標集合{1,3,8}令為F 有連通，以 黑色線 表之。2o 列指標集合{5}令為 o 3 E 與行指標集合{7}令為 o 3 F 有連通，以 藍色線 表之。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

依E₁o, Eo₂,, F₁o, F₂o,之順序，將列、行指標由小到大重新排序，可得出 IBD 矩陣T(C)。                        3 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 2 0 0 0 0 0 0 3 0 2 ) C ( T

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

參考文獻

[1] Arnold, B. C., and Press, S. J. (1989). Compatible Conditional Distributions. J. Amer. Statist. Assoc. 84, 152-156.

[2] Arnold, B. C., Castillo, E., and Sarabia, J. M. (2004), Compatibility of Partial or Complete Conditional Probability Specifications. J. Statist. Plann. Inference. 123, 133-159.

[3] Perez-Villalta, R. (2000), Variables finitas condicionalmente esecificadas. Questioo 24, 425-448.

[4] Song, C. C., Li, A., Chen, C. H., Jiang, T. J., and Kuo, K. L. (2010), Compatibility of finite discrete conditional distributions. Statistica Sinica. 20, 423-440.

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

附

錄

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

附錄 1： 程式使用說明 ( p50–p56)

將 4 個 MATLAB 程式檔 Main、CheckAB、IBD、RankOne 放入預設路徑 MATLAB 資料 夾下，開啟 MATLAB 軟體後，在命令視窗(Command Window)，將該視窗最大化。

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

命令視窗會顯示：請至 Excel 輸入 AB 矩陣資料? Y/N [Y]:

此時，離開命令視窗界面，並開啟 Excel 軟體在 Excel 界面下，輸入 A、B 矩陣的資料， 輸入的格式如下：

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

結果顯示如下： 在 A 這個頁面下，輸入的格式為數字，每個格子代表 A 矩陣中的每個元素，故 5x6 矩陣 則輸入 5 列與 6 行的格子，在 B 頁面也是如此。 這是 A 頁面

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

輸入完畢後，將檔案存檔，存檔位置或檔名可自行決定，本例是將檔案存在桌面並將檔 名改為 Set1。 這是 B 頁面

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

回到命令視窗界面，按下 Y，彈出視窗找到 Set1 檔案，點開啟。 程式會執行，並在命令視窗顯示結果。 下例為本文第 3 章的例 7

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

是否執行 IBD 程式檔? Y/N [Y]：按下 Y

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

是否求出對應 A、B 的聯合機率質量函數 J_AB? Y/N [Y]：按下 Y，並輸入比重

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

附錄 2： Main 主程式碼 ( p57–p60 )

clearall closeal clc formatrat format('compact')

reply_Excel = input('請至Excel輸入AB矩陣資料? Y/N [Y]: ','s');

if isempty(reply_Excel), reply_Excel ='Y';end if reply_Excel =='Y'

%========================== Data Import================================

[FileName,PathName,FilterIndex]=uigetfile('*.xls','Choose a dir'); A=xlsread([PathName FileName],'A');

B=xlsread([PathName FileName],'B'); %=================================================================== %=========================檢查Matrix是否符合型式======================= [C Index_Check]=CheckAB(A,B); %=================================================================== if Index_Check==1

reply_IBD = input('是否執行IBD程式檔? Y/N [Y]: ','s');

if isempty(reply_IBD), reply_IBD ='Y';end if reply_IBD =='Y'

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

%=============================== IBD=================================

[C_IBD E_set F_set T_Set IBD_Number]=IBD(C); C E_set F_set IBD_Number C_IBD %====================================================================

reply_RankOne = input('是否執行RankOne程式檔? Y/N [Y]: ','s');

if isempty(reply_RankOne), reply_RankOne ='Y';end if reply_RankOne =='Y' %=============================== RankOne ============================== C_IBD_bar =zeros(size(C)); col=1; row=1; reply_JAB='N'; for q=1:length(T_Set) sprintf('\b\b\b\b\b\b\n--- %s %s %s

---',['T',num2str(q)],['u',num2str(q)],['v',num2str(q)])

T=T_Set{q}

[T_ba{q} u_record{q} v_record{q}]=RankOne(T_Set{q});

C_IBD_bar (row+[1:size(T_ba{q},1)]-1,col+[1:size(T_ba{q},2)]-1)=T_ba{q}; row=row+size(T_ba{q},2);

col=col+size(T_ba{q},2);

if sum(sum(isnan(T_ba{q})))~=0

‧

國

立

政 治 大

學

‧

N

a

tio

na

l C

h engchi U

ni

ve

rs

it

y

else u=u_record{q}(:,end) v=v_record{q}(end,:) end sprintf('\b\b\b\b\b\b---') end Index_Check_RankOne=sum(sum(C_IBD_bar)); %================================ JAB ================================= if ~isnan(Index_Check_RankOne) C_IBD_bar

reply_JAB = input('是否求出對應A、B的聯合機率矩陣J_AB? Y/N [Y]: ','s');

if isempty(reply_JAB), reply_JAB ='Y';end end if reply_JAB =='Y' if IBD_Number==1 k_set=1; else for q=1:IBD_Number str=['請輸入比重k',num2str(q),':']; k_set(q) = str2num(input(str,'s')); end end k_set=k_set/sum(k_set); E=eye(length(E_set));E=E(E_set,:); IE=inv(E);