熊 雄,庹 恒
(1.湘潭大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院,湖南 湘潭 411105;2.中南大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,湖南 長沙 411000)
Cauchy分布[1]是概率統(tǒng)計(jì)問題中常用的分布,其密度為f(x)=λ/π(λ2+(x-μ)2),其中μ,λ為參數(shù),μ∈(-∞,+∞),λ>0.Cauchy分布在許多領(lǐng)域有廣泛的應(yīng)用.如在背景圖像中,像素點(diǎn)觀測的時(shí)序波動(dòng)來源于白噪聲,2幀圖像中對應(yīng)像素間的灰度比值的分布服從Cauchy分布.[2]但是因?yàn)镃auchy的各階原點(diǎn)矩均不存在,所以難以深入研究與數(shù)字特征有關(guān)的問題.雙截尾Cauchy分布的密度為
引入如下記號:記隨機(jī)變量X有分布函數(shù)F(x),并記
α(F)=inf{x:F(x)>0},ω(F)=sup{x:F(x)<1}.
引理1[5]設(shè)X1,n,X2,n,…,Xn,n為任意在[A,B]取值的連續(xù)型隨機(jī)變量的順序統(tǒng)計(jì)量,則對于?ε>0,有
引理2[6]若向量函數(shù)Ψ(X)存在m(m>2)階連續(xù)偏導(dǎo)數(shù),X是方程Ψ(X)=0的單根,當(dāng)初值X0充分接近精確解X時(shí),Newton迭代法收斂且至少為二階收斂.
序列{bn}取為bn=inf{x:1-F(x)≤1/n}.
引理7[10]設(shè)h(x;θ)為分布族,X1,X2,…,Xn為idd樣本,參數(shù)θ是一維的且參數(shù)空間Θ是開集,若lnh(x;θ)關(guān)于參數(shù)θ是可微的且分布族是可識(shí)別的(即對于?θ1,θ2∈Θ,{x:h(x;θ1)=h(x;θ2)}不是零測集),則以概率1似然方程在n→∞時(shí)有解,且此解關(guān)于θ是相合的.
假設(shè)參數(shù)μ,λ已知,考慮A,B的矩估計(jì),就是考慮如下方程組的解:
(1)
其中:
設(shè)由方程組(1)得到A,B的矩估計(jì)分別為AM,BM.首先,對于矩估計(jì)的可靠程度,由大數(shù)定律容易得到如下結(jié)論:
定理1矩法估計(jì)量AM,BM分別是參數(shù)A,B的弱相合估計(jì),即
(2)
要得到方程組(1)的公式解,其難度異乎尋常,因此筆者研究其數(shù)值解.下面以標(biāo)準(zhǔn)雙截尾Cauchy分布為例,說明其矩估計(jì)的算法步驟.
(ⅰ)待解方程組為
(ⅱ)尋找初值點(diǎn).由引理2和引理2,可以取(X1,n,Xn,n)作為計(jì)算(A,B)的初值點(diǎn).
(ⅲ)利用Newton迭代法計(jì)算.將fj在初值點(diǎn)(X1,n,Xn,n)進(jìn)行一階Taylor展開,得
(3)
忽略方程組(3)中的余項(xiàng),得
(4)
易知方程組(4)的解即為方程組(1)的近似解.
現(xiàn)考慮線性方程組(4)的解,其系數(shù)矩陣就是函數(shù)對于未知參數(shù)的Jacobi矩陣:
令f=(f1,f2)T,X=(A*,B*)T,X0=(X1,n,Xn,n)T,那么(2)式可以寫成f(X0)+J(X0)(X-X0)=0,其解為X1=X0-J-1(X0)f(X0).重復(fù)這個(gè)步驟,便得到如下Newton迭代法的公式:
Xk+1=Xk-J-1(Xk)f(Xk)k=0,1,2,….
Newton迭代法的優(yōu)勢在于平方范數(shù)收斂性,即‖Xk+1-X‖≤C‖Xk-X‖ .鑒于Newton迭代法的特性,剩余量范數(shù)不一定遞減.為了保證剩余量范數(shù)遞減,即‖f(Xk+1)‖≤‖f(Xk)‖(k=0,1,2,…),考慮加入一個(gè)修正因子ωk,ωk>0,將解改寫為
Xk+1=Xk-J-1ωk(Xk)f(Xk)k=0,1,2,….
選擇合適的ωk使得范數(shù)遞減,這個(gè)辦法被稱為阻尼Newton迭代法.
(ⅳ)確定算法的終止條件.由于阻尼Newton迭代法保證了剩余量范數(shù)遞減,因此只要找到某一步的運(yùn)算結(jié)果Xk并使得其滿足一定的條件即可.給定允許誤差值δ>0,若存在k∈N,s.t.‖Xk-X0‖≤δ,則有‖Xn-X0‖≤δ(?n≥k).
這里的矩估計(jì)算法只針對參數(shù)μ,λ已知的情形作出討論.事實(shí)上,若A,B已知,要求解參數(shù)μ,λ的矩估計(jì)時(shí),算法完全一致,只是初值點(diǎn)需要重新選擇.
定理2設(shè)Xk~C(μ,λ;A,B),X1,n,X2,n,…,Xn,n是其順序統(tǒng)計(jì)量,若λ,μ已知,則參數(shù)A,B分別有極大似然估計(jì)X1,n,Xn,n.
證明考慮對數(shù)似然函數(shù)
由X1,n≥A,Xn,n≤B,可得lnL(A,B)≤lnL(X1,n,Xn,n).再由lnx的單調(diào)性即得參數(shù)A,B的極大似然估計(jì)分別為X1,n,Xn,n.
證明由引理6,只要證明
現(xiàn)只要證明當(dāng)n→∞時(shí),以概率1似然方程有解.顯然雙截尾柯西分布的密度函數(shù)滿足引理7的條件,因其密度函數(shù)f(x)在定義域上連續(xù),故以概率1似然方程有解.
最后,先給出非標(biāo)準(zhǔn)雙截尾Cauchy分布極端順序統(tǒng)計(jì)量的漸近分布,再導(dǎo)出參數(shù)A,B的近似區(qū)間估計(jì).
bn=inf{x:1-F(x)≤1/n}.
證明對于Xk的分布函數(shù)F(x),有ω(F)<∞,則對于?x>0,由L'Hospital法則有
故取γ=2.由引理4,取an=0,bn=t,則
下面令ε→0,對于?x>1,有
再由引理4可得定理5成立.
cn,dn分別取為cn=μ,dn=sup{x:F(x)≤1/n}-μ.
證明對于雙截尾Cauchy分布的分布函數(shù)F(x),有α(F)=A,則分布函數(shù)F*(x)(F*(x)=F(A-1/x),x<0)滿足:對于 ?x>0,由L'Hospital法則有
故取γ=1.根據(jù)引理3,結(jié)合定理5的證明方法可得定理6成立.
定理7設(shè){Xk,1≤k≤n}~C(μ,λ;A,B),參數(shù)λ,μ已知,X1,n,X2,n,…,Xn,n為其順序統(tǒng)計(jì)量,則當(dāng)樣本量n充分大時(shí),參數(shù)A,B分別有近似區(qū)間估計(jì)
其中0<α<1為顯著性水平.
P{X1,n≤Ax+μ(1-x)}≈L2,1(x)n→∞.
(5)
特別地,取L2,1(x)=1-α,得x=-lnα.經(jīng)過整理,(5)式變成
[1] 鄧集賢,楊維權(quán),司徒榮,等.概率論及數(shù)理統(tǒng)計(jì)[M].4版.北京:高等教育出版社,2009:295-300.
[2] 明 英,蔣晶玨.視覺監(jiān)視中基于柯西分布的統(tǒng)計(jì)變化檢測[J].中國圖象圖形學(xué)報(bào),2008,13(2):328-334.
[3] 何永濟(jì).統(tǒng)籌方法中一任務(wù)的完成時(shí)間X應(yīng)服從截尾Cauchy分布[J].武漢水利電力學(xué)院學(xué)報(bào),1984(1):85-92.
[4] 匡能暉,陳 勇.雙截尾的Cauchy分布順序統(tǒng)計(jì)量的漸進(jìn)分布[J].北京大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,47(3):385-388.
[5] 林元烈,梁宗霞.隨機(jī)數(shù)學(xué)引論[M].北京:清華大學(xué)出版社,2003:123-136.
[6] 黃云清,舒 適.數(shù)值計(jì)算方法[M].北京:科學(xué)出版社,2004:158-168.
[7] NADARAJAH S.Explicit Expressions for Momentsof Order Statistic[J].Statistics & Probability Letters,2008,78(2):196-205.
[8] THOMAS P YAGEEN,SAMUEL PHILIP.Recurrence Relations for the Moments of Order Statistics from a BetaDdistribution[J].Statistical Papers,2008,49(1):139-146.
[9] 茆詩松.高等數(shù)理統(tǒng)計(jì)[M].北京:北京大學(xué)出版社,2003:123-136.
[10] 熊 雄,匡能暉.三參數(shù)的Pareto分布順序統(tǒng)計(jì)量的漸進(jìn)分布[J].四川大學(xué)學(xué)報(bào)(自然科學(xué)版),2012,49(5):975-978.