李春紅，黃紹軍，廖娟芬

（廣西大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院，廣西 南寧 030051）

IPF算法與N-R算法的對(duì)比研究及應(yīng)用

李春紅，黃紹軍，廖娟芬

（廣西大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院，廣西 南寧 030051）

IPF算法和N-R算法在列聯(lián)表的相關(guān)模型研究中應(yīng)用廣泛.本文通過MATLAB實(shí)現(xiàn)了IPF算法和N-R算法，并通過數(shù)據(jù)模擬研究了兩種算法的優(yōu)劣.模擬結(jié)果顯示，IPF算法具有穩(wěn)健、簡(jiǎn)潔等特點(diǎn)，而N-R算法在迭代次數(shù)、精度方面有優(yōu)勢(shì)，算法的實(shí)現(xiàn)和對(duì)比研究對(duì)實(shí)際應(yīng)用具有重要地參考意義.

列聯(lián)表；IPF算法；N-R算法

聯(lián)表在醫(yī)學(xué)、生物學(xué)、工農(nóng)業(yè)和社會(huì)科學(xué)等有著廣泛的應(yīng)用，引起人們不斷地去研究，提出了對(duì)數(shù)線性模型、Logistic回歸模型和多元分析模型等多種模型[1].這些分析模型應(yīng)用到高維列聯(lián)表時(shí)，可以把高維列聯(lián)表劃分為獨(dú)立模型和相關(guān)模型[2]，對(duì)于獨(dú)立模型可以直接用觀測(cè)數(shù)據(jù)進(jìn)行估計(jì)，而相關(guān)模型則不能直接用觀測(cè)數(shù)據(jù)來估計(jì)，于是人們提出了各種迭代算法.常用的算法有，迭代比例擬合算法（Iterative Proportional Fitting Algorithm，IPF）[2-3]和牛頓—勒夫遜算法（Newton-Raphson，N-R）[1].文獻(xiàn)[2]用分別IPF算法和N-R算法，文獻(xiàn)[3]提出了IPF算法對(duì)三維列聯(lián)表進(jìn)行擬合，文獻(xiàn)[4]給出了在對(duì)高維列聯(lián)表進(jìn)行探索性分析時(shí)IPF算法比較簡(jiǎn)單、穩(wěn)健，但并沒有深入的對(duì)這兩種算法進(jìn)行比較分析.本文利用Matlab對(duì)兩種算法進(jìn)行模擬，通過比較分析，闡明了IPF算法與N-R算法的優(yōu)劣性．

1 算法簡(jiǎn)介

假設(shè)有n個(gè)個(gè)體根據(jù)三個(gè)屬性X、Y和Z進(jìn)行分類.屬性X有r類，記為X1，…，Xr，屬性Y有c類，記為Y1，…，Yc，屬性Z有t類，記為Z1，…，Zt.n個(gè)個(gè)體中屬于Xi、Yj和Zk類的有nijk個(gè)，比如表1中，n111=22，根據(jù)這三個(gè)屬性繪制一個(gè)交叉表，就得到一張三維的r×c×t列聯(lián)表.例如在研究手術(shù)后并發(fā)癥是否與手術(shù)類型有關(guān)的課題中，得到見表1.

表1 并發(fā)癥與手術(shù)類型案例數(shù)據(jù)Fig.1 Cases of complications with operation types

這就是一張2×2×2三維列聯(lián)表.實(shí)際上，屬于Xi、Yj和Zk類的個(gè)體是以一定值的概率發(fā)生，其概率值記為Pijk.那么于全部個(gè)體n，屬于Xi、Yj和Zk類的理論個(gè)數(shù)為npijk，該理論值稱為期望頻數(shù)，記mijk=npijk.對(duì)應(yīng)的，實(shí)際屬于Xi、Yj和Zk類的個(gè)數(shù)nijk稱為實(shí)際頻數(shù)，并且有Enijk=mijk.

如果兩個(gè)屬性X和Y之間毫無關(guān)系，就稱X和Y相互獨(dú)立，反之就稱X和Y相關(guān).在列聯(lián)表中，如果兩個(gè)屬性X和Y之間存在相關(guān)關(guān)系，就稱X和Y之間有二次交互作用，如果兩個(gè)屬性X、Y和Z之間存在相關(guān)關(guān)系，就稱X、Y和Z之間有三次交互作用．

三維列聯(lián)表的相關(guān)模型有兩種.一種是飽和模型，其期望頻數(shù)不能分解，此模型三個(gè)屬性X、Y和Z不僅兩兩之間有交互作用，而且三個(gè)之間也有交互作用，其期望頻數(shù)的估計(jì)值就是實(shí)際頻數(shù)；另一種是期望頻數(shù)可以分解，此模型三個(gè)屬性X、Y和Z僅兩兩之間有交互作用，但三者之間沒有交互作用.可分解模型可記為mijk=αijβikγjk，其中αij、βik和γjk分別表示（X，Y）、（X，Z）和（Y，Z）之間交互作用，可分解相關(guān)模型的期望頻數(shù)的估計(jì)值要求用迭代算法來估計(jì)，常用的有IPF算法和N-R算法.下面就這兩個(gè)算法進(jìn)行簡(jiǎn)單介紹．

1.1 IPF算法簡(jiǎn)介

對(duì)于可分解相關(guān)模型mijk=αijβikγjk，求其極大似然估計(jì)有代估計(jì)就是期望頻數(shù)mijk的極大似然估計(jì)mijk．

1.2 N-R算法簡(jiǎn)介

N-R算法是求解非線性方程組的常用方法，對(duì)于相關(guān)模型mijk=αijβikγjk同樣可以通過N-R算法來估計(jì).相關(guān)模型 mijk=αijβikγjk中含有的獨(dú)立參數(shù)個(gè)數(shù)可以通過對(duì)數(shù)線性模型變換來體現(xiàn)，其對(duì)數(shù)線性模型的表達(dá)如下：

2 數(shù)據(jù)模擬

本文通過Matlab數(shù)學(xué)軟件進(jìn)行數(shù)值模擬.模擬分為兩部分，第一部分求N-R算法的有效區(qū)間，第二部分是模擬比較兩種算法.第二部分模擬了100次，模擬過程分為三步驟進(jìn)行，第一步由Poison分布隨機(jī)生成一個(gè)三維的列聯(lián)表，第二步分別IPF算法和N-R算法對(duì)列聯(lián)表進(jìn)行數(shù)值計(jì)算，第三步對(duì)計(jì)算結(jié)果進(jìn)行組合整理后輸出數(shù)值結(jié)果.模擬具體過程如下.

2.1 求解有效區(qū)間的算法代碼及輸出結(jié)果

說明，Meandx為平均有效區(qū)間的輸出結(jié)果，即N-R算法的平均有效區(qū)間為(β0-0.4089,β0+3.8890),β0為初始值.

2.2 IPF算法與N-R算法模擬比較的代碼及輸出結(jié)果

說明，Meanlter、SumOutChi、SumOutLh 都是結(jié)果輸出.Meanlter中6.3000是N-R算法的平均迭代次數(shù)，8.9200是IPF算法的平均迭代次數(shù)；SumOut?Chi是兩種算法通過卡方檢驗(yàn)值比較得出的計(jì)算結(jié)果，輸出結(jié)果有6個(gè)元素，分別表示為，①模擬無效個(gè)數(shù)；②IPF卡方檢驗(yàn)值大于N-R卡方檢驗(yàn)值的個(gè)數(shù)；③卡方檢驗(yàn)值相等的個(gè)數(shù)；④IPF卡方檢驗(yàn)值小于N-R卡方檢驗(yàn)值的個(gè)數(shù)；⑤卡方檢驗(yàn)值的差值總和；⑥卡方檢驗(yàn)值差值比率的總和（以N-R值為基數(shù)），即對(duì)應(yīng)的結(jié)果為[0 94 0-6 139.5953 0.2380]；同樣，SumOutLh是兩種算法通過似然比值比較得出的計(jì)算結(jié)果，其各個(gè)元素的代表的意義跟SumOutChi一樣.

3 結(jié)果分析

即只要初始值不為零且相等，就進(jìn)行迭代運(yùn)算．

2）從迭代次數(shù)來看，IPF算法平均迭代次數(shù)為8.9200，N-R算法在β0處的平均迭代次數(shù)為6.3000.由此我們可以看出，N-R算法的迭代次數(shù)在初始值β0出具有較強(qiáng)的收斂速度．

3）從檢驗(yàn)值的比較來看，IPF算法的似然比值大于N-R算法似然比值所占的比率為99%，IPF算法的似然比值小于N-R算法所占的比率為1%，而且IPF算法的似然值比N-R算法似然比值差值比率的平均值為0.9471%，可見N-R算法比IPF算法在似然比值的標(biāo)準(zhǔn)下有優(yōu)勢(shì).同樣的，通過卡方檢驗(yàn)值來比較也得出相同的效果．

4）如果要應(yīng)用到實(shí)際的數(shù)據(jù)，去掉生成矩陣A的函數(shù)，把A設(shè)定為主函數(shù)輸入?yún)?shù)即可.對(duì)于N-R算法，當(dāng)有的觀測(cè)值為0時(shí)，有可能不收斂，所以建議把觀測(cè)值都加上0.5再進(jìn)行運(yùn)算．

4 小結(jié)

本文通過Matlab來對(duì)比分析IPF算法與N-R算法，得出了IPF算法具有穩(wěn)健、簡(jiǎn)潔等特點(diǎn)，而N-R算法在迭代次數(shù)、精度方面有優(yōu)勢(shì).就擬合結(jié)果來看，IPF與N-R算法分析結(jié)果基本上一致，但N-R的能夠詳細(xì)的給出參數(shù)的估計(jì)值，因此本文提倡用N-R算法.另外，對(duì)于SAS、SPSS等統(tǒng)計(jì)軟件不熟悉的科研工作者，可以參考利用本文代碼在Matlab上進(jìn)行數(shù)值求解．

[1]張堯庭,夏立顯,安希忠,等.定性資料的統(tǒng)計(jì)分析[M].廣西:廣西師范大學(xué)出版社,1991:78-85.

[2]王靜龍,梁小筠.定性數(shù)據(jù)統(tǒng)計(jì)分析[M].北京:中國統(tǒng)計(jì)出版社,2008:147-155.

[3]Santner T J,Duffy D E.The Statistical Analysis of Discrete Data[M].New York，Spring-Verlag,1989:113-199．

[4]張巖波,何大衛(wèi).對(duì)數(shù)線性模型的IPF算法及其軟件實(shí)現(xiàn)[J].中國衛(wèi)生統(tǒng)計(jì),1999,16(5):258-259．

Comparative Analysis and Application of IPF Algorithm and N-R Algorithm

LI Chunhong,HUANG Shaojun,LIAO Juanfen
（College of Mathematics and Information Science，Guangxi University，Nanning530004，China）

IPF algorithm and N-R algorithm are widely used in the dependence modeling of contingency table.This ar?ticle wrote IPF algorithm and N-R algorithm codes in the MATLAB Software and compared the two algorithms by means of digital simulation.The simulation showed that IPF algorithm hold the robustness and briefness features,N-R algo?rithm had advantages in iterations and precision.The codes and conclusion offered important reference value in practical application.

contingency table；IPF algorithm；N-R algorithm

O 212.1

A

1674-4942（2011）02-0128-03

2011-01-25

國家自然基金資助項(xiàng)目（11061002）

畢和平