2013年蘆山MS7．0地震序列參數(shù)的早期特征：傳染型余震序列模型計算結(jié)果*

蔣長勝 莊建倉 龍 鋒 韓立波 郭路杰

1）中國北京 100081 中國地震局地球物理研究所

2）日本東京 106-8569 數(shù)理統(tǒng)計研究所

3）中國成都 610041 四川省地震局

引言

據(jù)中國地震臺網(wǎng)測定，2013年4月20日8時2分46．0秒，四川省蘆山縣發(fā)生MS7．0強烈地震．至4月23日6時，該地震已造成193人死亡，25人失蹤，1萬2 211人受傷（中華人民共和國中央人民政府，2013）．蘆山MS7．0地震是繼2008年汶川MS8．0地震后在龍門山斷裂帶上發(fā)生的又一次災(zāi)害性地震，兩次強烈地震的余震區(qū)相距僅約45km（劉杰等，2013）．因涉及到與汶川MS8．0地震的關(guān)系（陳運泰等，2013）及龍門山斷裂帶南段的地震危險性等問題，蘆山地震的發(fā)生引起了社會、政府和地震學(xué)界的高度關(guān)注．

余震序列的衰減特征和激發(fā)余震的能力可通過序列的統(tǒng)計參數(shù)來表征，不同構(gòu)造區(qū)域、序列類型、主震的震源機制類型、震源區(qū)局部構(gòu)造應(yīng)力水平及大地?zé)崃鞯鹊厍蛭锢硖卣鞫伎赡鼙憩F(xiàn)為序列參數(shù)的顯著差異（Kaganetal，2010）．對余震序列參數(shù)尤其是早期參數(shù)特征（蔣海昆等，2007）的研究一直是地震學(xué)和地球物理學(xué)研究的熱點問題．近年來，作為“大森－宇津”公式（Omori，1894；Utsu，1961）的推廣，“傳染型余震序列”（epidemic-type aftershock sequence，簡寫為ETAS）模型（Ogata，1988，1989，2001）得到長足發(fā)展．ETAS模型除可對地震序列參數(shù)進行接近真實的估計外，還被用于地震活動平靜識別（Ogata，1992，2001；Zhuang，2000）與地震活動的觸發(fā)（Pengetal，2012）等研究．

與蘆山MS7．0地震相關(guān)的一個重要科學(xué)問題，就是地震序列參數(shù)及其表征的震源區(qū)地震學(xué)特征．本文使用ETAS模型和最大似然法對蘆山MS7．0地震序列進行參數(shù)估計，并討論參數(shù)估計的可靠性以及參數(shù)的地震學(xué)含義．

1 時間序列ETAS模型

余震序列隨時間的衰減常用修正的“大森－宇津”公式（Omori，1894；Utsu，1961）來表示：

式中，f（t）為t時刻對應(yīng)的余震發(fā)生率，t為主震發(fā)生后的離逝時間，μ表示背景地震的發(fā)生率，p表示余震序列衰減的快慢，c為主震后余震頻次達到峰值時對應(yīng)的時間，K表示余震的活躍程度．式中的μ，K，c和p均為非負常數(shù)．

實際上，真實的余震序列遠比簡單的修正的“大森－宇津”公式復(fù)雜．針對每次余震都可激發(fā)自己的余震的情況，Ogata（1988，1989，2001）在“大森－宇津”公式基礎(chǔ)上，提出了基于隨機點過程的ETAS模型．ETAS模型假設(shè)研究區(qū)的背景地震發(fā)生率μ為常數(shù)，所有余震遵從“大森－宇津”公式激發(fā)自己的余震，被激發(fā)的余震數(shù)僅與其“母震”的震級相關(guān)，且震級的分布是獨立的．ETAS模型假定主震的發(fā)生為初始零時刻，在其后的一個觀測時間段［0，T］內(nèi)的地震序列｛（ti，Mi）；i＝1，2，…，N｝的強度函數(shù)可表示為（Ogata，1988）

其中，Mi和ti分別表示第i個事件的震級和發(fā)生時間；M0為參考震級，一般可取截止震級．p與修正的“大森－宇津”公式中的p有相同的物理含義，均表示序列衰減的快慢，p越大衰減越快，反之亦然．α表示觸發(fā)次級余震的能力（Ogata，1989，1992）．對于震群型序列，一般情況下α＜1，而序列中無明顯的次級余震時，α則較大（Ogata，2001）．

ETAS模型參數(shù)可使用最大似然法進行估計．由于主震發(fā)生后短時間內(nèi)主震的波形振幅較大，使得后續(xù)余震檢測的信噪比顯著降低，余震頻次難以達到峰值，因此擬合常在［S，T］范圍內(nèi)進行，其中0＜S＜T．似然函數(shù)L的形式為

其中，待估參數(shù)為μ，K，c，α，p．盡管上式中的參數(shù)在［S，T］內(nèi)進行擬合，但實際上［0，S］范圍內(nèi)的事件也參與了計算，具體體現(xiàn)在上式右側(cè)的λ（t）中（Ogata，2001）．對于模型的適用性檢測，常使用赤池準(zhǔn)則（Akaike information criteria，簡寫為AIC）（Akaike，1974）進行判斷：

AIC值越小表示適用效果越好．值得注意的是，上式中k為待定擬合參數(shù)的數(shù)量．因此，AIC準(zhǔn)則既考慮了模型對數(shù)據(jù)的適用程度，也考慮了對增加參數(shù)個數(shù)提高擬合效果的“懲罰”．ETAS模型研究中常使用AIC準(zhǔn)則判斷不同模型的擬合效果（Guo，Ogata，1997），以及對擬合變點問題（change-point problem）的討論（Ogata，1992；Zhuang，2000）等．

2 研究所用資料和序列目錄的完整性分析

本研究所用的蘆山MS7．0地震序列目錄來自中國地震臺網(wǎng)中心“全國地震編目系統(tǒng)”提供的《四川地震臺網(wǎng)速報目錄》①全國地震編目系統(tǒng)．http：∥10．5．202．22／bianmu／index．jsp．查閱日期：2013年5月15日．．選用2013年4月1日—5月14日0級以上地震，限定地震序列的空間范圍為29．8°—30．6°N，102．6°—103．2°E．

研究中采用“震級－序號”法（Huang，2006；蔣長勝，吳忠良，2011；Zhuangetal，2012）定性討論蘆山地震序列目錄的完整性震級Mc．震級－序號法按地震發(fā)生時間的先后順序排序，地震密度較大的區(qū)域的連線大致為Mc的時序變化．這種分析一方面避免了主震發(fā)生后短期內(nèi)余震較為密集，在時間軸上難以分析余震監(jiān)測能力的快速變化；另一方面，由于監(jiān)測分析工作中的人為因素，Mc在震后短期內(nèi)常出現(xiàn)分段和不連續(xù)性的變化，震級－序號法則易于分辨這種人為因素（圖1）．由圖1a給出的震級－序號點圖和放大的子圖可見，蘆山地震發(fā)生后的短時間內(nèi)，Mc僅為ML3．0左右，其后余震區(qū)的監(jiān)測能力逐漸提升，Mc逐漸減??；至第190個事件，即震后0．31天左右，Mc可達ML2．0．

圖1 蘆山MS7．0地震序列目錄的完整性分析（a）蘆山MS7．0地震序列的震級－序號圖．子圖為序列第1—240個事件的放大圖，垂直虛線標(biāo)出了ML2．0以上地震完整的位置，約為第190個事件，即主震后0．31天；（b）震級－序號法給出的地震密度分布Fig．1 Catalogue completeness of the Lushan MS7．0earthquake sequence（a）Magnitude-rank distribution of the Lushan MS7．0earthquake sequence．Grey subplot shows the partial enlarged view of the first 240events，and the vertical dashed line indicates the completeness of ML≥2．0aftershocks，which is the position of the 190th event or 0．31days after the main shock；（b）Seismic rate distribution given by magnitude-rank analysis

3 蘆山MS7．0地震序列的ETAS模型參數(shù)的最大似然估計

在對蘆山MS7．0地震序列進行ETAS模型參數(shù)的最大似然估計時，選定2013年4月1日—5月14日的地震序列，設(shè)定Mc＝ML2．0，參數(shù)擬合所用的時段為震后0．31—24．12天（5月14日24時）．由此獲得的蘆山MS7．0地震序列的ETAS模型參數(shù)為μ＝0．0000，K＝0．009 5，c＝0．045 7，α＝1．886 4和p＝1．220 5，結(jié)果如圖2所示．圖2a為利用最大似然估計獲得的ETAS模型條件強度曲線，即單位時間內(nèi)地震發(fā)生率隨時間的變化．由該圖可見，在震后10天內(nèi)，序列衰減較為平穩(wěn)；10天后由于序列地震事件的顯著減少，較大余震事件引起的強度曲線變化較為明顯．

為檢測ETAS模型對實際地震序列的擬合效果，除AIC準(zhǔn)則外，一般還常使用“轉(zhuǎn)換時間”域的“殘差分析”（Ogata，1988；Daley，Vere-Jones，2003）．對條件強度函數(shù)λ（t）可采用下式對時間序列｛ti｝進行轉(zhuǎn)換：

根據(jù)上式，｛ti｝被逐個地轉(zhuǎn)化為｛τi｝，且｛τi｝服從單位速率的穩(wěn)態(tài)泊松分布（Zhuanget al，2012），經(jīng)過轉(zhuǎn)化的數(shù)據(jù)稱為“殘差點過程”（residual point process，簡寫為RPP）．如果模型與數(shù)據(jù)擬合得較好，RPP在轉(zhuǎn)換時間－累積圖上就表現(xiàn)為線性，接近標(biāo)準(zhǔn)的穩(wěn)態(tài)泊松過程的理論直線．圖3a給出了蘆山MS7．0地震序列累積地震數(shù)與ETAS擬合曲線在“轉(zhuǎn)圖2 利用ETAS模型擬合給出的蘆山MS7．0地震序列ML2．0以上地震的條件強度曲線（a）和M-t圖（b）．縱坐標(biāo)上的頻度指每天發(fā)生的地震次數(shù)

Fig．2 Temporal variation of the conditional intensity from fitting the ETAS model to the LushanMS7．0earthquake sequence with cutoff magnitude ML2．0（a）andM-tplot（b）．The frequency on the ordinate axis means the earthquake number per day換時間”域τ的比較情況，即RPP分析結(jié)果．由圖3b給出的M-τ圖可見，蘆山MS7．0地震序列在“轉(zhuǎn)換時間”域τ已被轉(zhuǎn)換為穩(wěn)態(tài)泊松分布．實際上就一般強余震預(yù)測而言，對“轉(zhuǎn)換時間”域τ的預(yù)測和預(yù)測效能進行檢驗可能更有意義（Jiang，Wu，2012）．

圖3 蘆山MS7．0地震序列ML2．0以上地震的累積地震數(shù)與ETAS模型擬合曲線的比較（a）累積地震數(shù)與ETAS擬合曲線在“轉(zhuǎn)換時間”域τ的比較；（b）M-τ圖，τ為將時間轉(zhuǎn)換為穩(wěn)態(tài)泊松分布的“轉(zhuǎn)換時間”Fig．3 Comparison of the observed and modeled cumulative numbers of the Lushan MS7．0earthquake sequence with magnitude larger than ML2．0．（a）Observation（black curve）and ETAS formula（thick gray dashed line）in the transformed time（τ）domain；（b）M-τplot whereτis the transformed time which is according to the stationary Poisson process

4 序列參數(shù)早期特征的影響因素分析

由于真實地震序列的復(fù)雜性，參數(shù)估計的不確定性一直是ETAS模型研究關(guān)注的重要問題．例如截止震級對序列分支比（Sornette，Werner，2005）、模型參數(shù)的估計偏差（Schoenbergetal，2010）的影響，以及地震序列時空選取對ETAS模型參數(shù)估計的影響（Wangetal，2010a）等．為考察蘆山MS7．0地震序列ETAS模型參數(shù)早期特征的穩(wěn)定性，本研究從截止震級Mc的選取及擬合截止時間兩個角度分別討論．

4．1 截止震級對ETAS模型參數(shù)估計的影響

由于ETAS模型假定每一次地震都可觸發(fā)其后發(fā)生的任何震級的地震，因此，截止震級Mc在確保地震目錄完整性的同時，也可能切斷了Mc以下地震與以上地震的觸發(fā)關(guān)系，從而造成“鏈接缺失”現(xiàn)象（Wangetal，2010a）．為考察Mc的選取對ETAS模型參數(shù)估計的影響，研究中設(shè)定Mc＝ML2．0，2．1，…，3．0，分別考察模型參數(shù)的最大似然估計結(jié)果．此外，還采用基于對數(shù)似然函數(shù)的Hessian矩陣（Ogata，1978）估計擬合參數(shù)的標(biāo)準(zhǔn)差，結(jié)果見表1．

由表1可見，當(dāng)Mc自ML2．0—3．0逐漸升高時，參數(shù)K總體有減小的趨勢，α有明顯的線性增加，表明序列觸發(fā)次級余震的能力隨著截止震級的增加逐漸變?nèi)酰畃值則隨Mc的變化不大，約為1．22—1．26，表明在蘆山MS7．0地震序列早期階段，序列衰減快的特征較為穩(wěn)定．

表1 不同截止震級下ETAS模型擬合參數(shù)Table 1 Parameters of ETAS model with different cutoff magnitude

4．2 擬合截止時間對ETAS模型參數(shù)估計的影響

由于板內(nèi)地震所處構(gòu)造運動速率不同，地震序列持續(xù)時間存在明顯差別（Stein，Liu，2009），震源的破裂特征、斷層愈合速度、周邊構(gòu)造場的調(diào)整等諸多方面也會使地震序列的衰減特征、激發(fā)余震能力等出現(xiàn)明顯差異．因此，所用序列的長度或擬合截止時間不同，獲得的ETAS模型參數(shù)可能不同，擬合參數(shù)的標(biāo)準(zhǔn)差也將存在差異（Wangetal，2010b）．為考察震后短期內(nèi)序列參數(shù)早期特征的穩(wěn)定性，設(shè)定序列擬合的截止時間tend＝［2．0，2．5，3．0，3．5，4．0，5．0，…，23，24．12］，利用最大似然法分別進行參數(shù)估計．

除考察ETAS模型主要參數(shù)α和p外，研究中還考慮了相應(yīng)的地震序列b值的穩(wěn)定性．在b值計算中采用最大似然法（Aki，1965；Utsu，1965）：

式中，〈M〉為平均震級；ΔMbin為地震目錄的震級滑動寬度，一般ΔMbin＝0．1．對b值的不確定度δb估計采用Shi和Bolt（1982）給出的方法：

計算獲得的α值、p值和b值結(jié)果如圖4所示．由該圖可見，在震后10天內(nèi)上述參數(shù)均存在顯著變化，標(biāo)準(zhǔn)差幅值較大．其中α值變化范圍為2．47—2．10，p值變化范圍為0．93—1．24，b值則逐漸由0．64增加至0．71．在震后10—24．12天（5月14日），α值逐漸減小至1．89，p值逐漸增加至1．22，b值相對平穩(wěn)地保持在0．72左右．由此可見，在估算蘆山MS7．0地震序列的ETAS模型參數(shù)時，震后10天內(nèi)獲得的結(jié)果隨時間變化較為明顯，而其后各參數(shù)盡管有所變化但總體較為平穩(wěn)．

圖4 ETAS模型參數(shù)α，p和b隨擬合截止時間的變化（a）和M-tend圖（b）Fig．4 Parametersα，pand b against the end time tendin ETAS fitting（a）and the M-tendplot（b）

4 討論與結(jié)論

為考察2013年4月20日蘆山MS7．0地震序列參數(shù)的早期特征，本文利用ETAS模型和最大似然法進行了參數(shù)估計．選用截止震級Mc＝ML2．0，擬合所用的時段為震后0．31—24．12天，計算得到α＝1．89和p＝1．22．序列總體表現(xiàn)為觸發(fā)次級余震的能力較弱，序列衰減速率較快；利用最大似然法估計的b值也較低，約為0．72，表明余震區(qū)應(yīng)力水平相對較高．

為檢測上述結(jié)果穩(wěn)定性，設(shè)定Mc＝ML2．0，2．1，…，3．0，以及選用不同的擬合截止時間分別進行參數(shù)擬合和參數(shù)標(biāo)準(zhǔn)差估計．結(jié)果表明，截止震級Mc的選取對α和K值影響明顯，但對p值影響則較?。送猓鸷?0天內(nèi)得到的參數(shù)擬合結(jié)果隨時間變化較為明顯，而其后各參數(shù)變化總體較為平穩(wěn)．

根據(jù)蔣海昆等（2007）給出的中國大陸294次M＞5．0中強震余震序列參數(shù)ETAS模型擬合結(jié)果，西南地區(qū)的平均序列參數(shù)為b＝0．871，α＝0．933，p＝0．760；中國大陸地區(qū)斜滑型主震對應(yīng)的平均序列參數(shù)為b＝0．865，α＝0．979，p＝0．876；中國大陸M＞7．0地震的平均序列參數(shù)為b＝0．956，α＝0．509，p＝0．684．與上述結(jié)果比較并考慮結(jié)果的離散范圍，本研究獲得的蘆山MS7．0地震序列b＝0．72，低于上述相同區(qū)域（西南地區(qū)）、主震類型（逆沖型）和主震震級（M＞7．0）的平均結(jié)果；α＝1．89顯著高于上述結(jié)果，并高于日本地區(qū)板內(nèi)地震的平均結(jié)果（Guo，Ogata，1997），表明該地震對次級地震的激發(fā)能力較弱；此外p＝1．22也明顯高于上述各結(jié)果，表明序列衰減速率較快．

本文采用時間序列ETAS模型并重點研究了地震序列參數(shù)的早期特征．計算結(jié)果顯示蘆山MS7．0地震序列參數(shù)與以往研究的中國大陸震例存在顯著差異，可能有兩方面原因：一方面是龍門山斷裂帶南段的構(gòu)造特點和介質(zhì)屬性，但與2008年汶川MS8．0地震衰減速率和斷層愈合過程相比差別較大來看，這種影響可能性很小；另一方面是受區(qū)域構(gòu)造應(yīng)力場的影響，此次地震的序列參數(shù)顯示了區(qū)域較強的構(gòu)造應(yīng)力背景，而這種影響可能造成斷層愈合過程較快，或發(fā)生強余震、后續(xù)強震的可能性較大．但對斷層愈合過程的討論仍需更多直接的觀測證據(jù)證明．實際上，目前已有多種形式的ETAS模型，例如，基于隨機除叢法的時－空ETAS模型（Ogata，1998；Zhuangetal，2002；Zhuang，Ogata，2006），考慮變化的背景地震活動的時間序列ETAS模型（Pengetal，2012）等；本文中的ETAS模型選取主要出于使研究問題盡量簡化的考慮．此外，在參數(shù)估計中未考慮地震目錄中震級等的測定誤差，這可能是本研究的不足之處．

本文使用了中國地震臺網(wǎng)中心“全國地震編目系統(tǒng)”提供的四川省區(qū)域地震臺網(wǎng)速報目錄，中國地震局“四川省蘆山‘4·20’7．0級強烈地震科學(xué)考察”總指揮吳忠良研究員對本文給予了指導(dǎo)，周仕勇教授和蔣海昆研究員提出了有益建議．作者在此一并表示感謝．

陳運泰，楊智嫻，張勇，劉超．2013．從汶川地震到蘆山地震［J］．中國科學(xué)：地球科學(xué)，43（6）：1064-1072．

蔣長勝，吳忠良．2011．玉樹MS7．1地震前的中長期加速矩釋放（AMR）［J］．地球物理學(xué)報，54（6）：1501-1510．

蔣海昆，鄭建常，吳瓊，曲延軍，李永莉．2007．傳染型余震序列模型震后早期參數(shù)特征及其地震學(xué)意義［J］．地球物理學(xué)報，50（6）：1778-1786．

劉杰，易桂喜，張致偉，官致軍，阮祥，龍鋒，杜方．2013．2013年4月20日四川蘆山MS7．0級地震介紹［J］．地球物理學(xué)報，56（4）：1404-1407．

中華人民共和國中央人民政府．2013．民政部：蘆山地震已造成193人死亡，25人失蹤［EB／OL］．［2013-04-23］．http：∥www．gov．cn／jrzg／2013-04／23／content＿2386791．htm．

Akaike H．1974．A new look at the statistical model identification［J］．IEEETransAutomatControl，19（6）：716-723．Aki K．1965．Maximum likelihood estimate ofbin the formula logn＝a－bmand its confidence limits［J］．BullEarthquake ReInst，43（2）：237-239．

Daley D D，Vere-Jones D．2003．AnIntroductiontoTheoryofPointProcesses——Volume1：ElementaryTheoryand Methods［M］．Second Edition．Springer：New York：17-33．

Guo Z，Ogata Y．1997．Statistical relations between the parameters of aftershocks in time，space and magnitude［J］．J GeophysRes，102（B2）：2857-2873．

Huang Q．2006．Search for reliable precursors：A case study of the seismic quiescence of the 2000Western Tottori prefecture earthquake［J］．JGeophysRes，111（B4）：B04301，doi：10．1029／2005JB003982．

Jiang C S，Wu Z L．2012．Testing the forecast of aftershocks：A simple method with an example of application［J］．ResearchinGeophysics，2（1）：29-33．

Kagan Y Y，Bird P，Jackson D D．2010．Earthquake patterns in diverse tectonic zones of the globe［J］．PureApplGeophys，167（6／7）：721-741．

Ogata Y．1978．The asymptotic behavior of maximum likelihood estimators for stationary point processes［J］．AnnInst StatistMath，30（2）：243-261．

Ogata Y．1988．Statistical models for earthquake occurrences and residual analysis for point processes［J］．JAmerStatist Assoc，83（401）：9-27．

Ogata Y．1989．Statistical model for standard seismicity and detection of anomalies by residual analysis［J］．Tectonophysics，169（1／3）：159-174．

Ogata Y．1992．Detection of precursory relative quiescence before great earthquakes through a statistical model［J］．J GeophysRes，97（B13）：19845-19871．

Ogata Y．1998．Space-time point-process models for earthquake occurrences［J］．AnnInstStatistMath，50（2）：379-402．

Ogata Y．2001．Increased probability of large earthquakes near aftershock regions with relative quiescence［J］．JGeophys Res，106（B5）：8729-8744．

Omori F．1894．On aftershocks of earthquakes［J］．JCollSciImpUnivTokyo，7：111-200．

Peng Y J，Zhou S Y，Zhuang J C．2012．An approach to detect the abnormal seismicity increase in Southwestern China triggered co-seismically by 2004SumatraMW9．2earthquake［J］．GeophysJInt，189（3）：1734-1740．

Schoenberg F P，Chu A，Veen A．2010．On the relationship between lower magnitude thresholds and bias in ETAS parameter estimates［J］．JGeophysRes，115（B4）：B04309，doi：10．1029／2009JB006387．

Shi Y，Bolt B A．1982．The standard error of the magnitude frequencyb-value［J］．BullSeismolSocAm，72（5）：1677-1687．

Sornette D，Werner M J．2005．Apparent clustering and apparent background earthquakes biased by undetected seismicity［J］．JGeophysRes，110（B9）：B09303，doi：10．1029／2005JB003621．

Stein S，Liu M．2009．Long aftershock sequences within continents and implications for earthquake hazard assessment［J］．Nature，462：87-89．

Utsu T．1961．A statistical study on the occurrence of aftershocks［J］．GeophysMag，30：521-605．

Utsu T．1965．A method for determining the value ofbin a formula logn＝a－bmshowing the magnitude frequency for earthquakes［J］．GeophysBullHokkaidoUniv，13：99-103．

Wang Q，Jackson D D，Zhuang J C．2010a．Missing links in earthquake clustering models［J］．GeophysResLett，37（21）：L21307，doi：10．1029／2010GL044858．

Wang Q，Schoenberg F P，Jackson D D．2010b．Standard errors of parameter estimates in the ETAS model［J］．Bull SeismolSocAm，100（5A）：1989-2001．

Zhuang J C．2000．Statistical modeling of seismicity patterns before and after the 1990Oct 5Cape Palliser earthquake，New Zealand［J］．NZJGeolGeophys，43（3）：447-460．

Zhuang J，Harte D，Werner M J，Hainzl S，Zhou S．2012．BasicModelsofSeismicity：TemporalModels［M／OL］．Community Online Resource for Statistical Seismicity Analysis．doi：10．5078／corssa-79905851．http：∥www．corssa．org．

Zhuang J，Ogata Y．2006．Properties of the probability distribution associated with the largest event in an earthquake cluster and their implications to foreshocks［J］．PhysicalReviewE，73（4）：046134，doi：10．1103／PhysRevE．73．046134．

Zhuang J，Ogata Y，Vere-Jones D．2002．Stochastic declustering of space-time earthquake occurrences［J］．JAmerStat Assoc，97（458）：369-380．