李國(guó)重,夏云青,焦慧平
(1.信息工程大學(xué)理學(xué)院,鄭州 450001;2.中州大學(xué)信息工程學(xué)院,鄭州 450044)
在對(duì)加工工件的測(cè)量數(shù)據(jù)處理中,經(jīng)常需要對(duì)新的測(cè)量平差理論或方法進(jìn)行評(píng)估,其中一項(xiàng)重要工作是檢驗(yàn)新理論或新方法抵御粗差的能力 (粗差是指離群的大誤差[1,2],由失誤等引起,它實(shí)際不可避免)。如果新理論或新方法能夠有效地抵御粗差的不良影響,那么認(rèn)為新理論或新方法有效,從而在測(cè)量數(shù)據(jù)處理中可以大膽使用;否則,我們懷疑新理論或新方法的效能,需謹(jǐn)慎使用或者棄之不用[3]。理論上講,一個(gè)性能優(yōu)良的平差理論能夠抵御觀測(cè)數(shù)據(jù)的任何粗差,無(wú)論這些粗差分布在哪些觀測(cè)數(shù)據(jù)中以及粗差的量級(jí)多大。然而,檢驗(yàn)平差理論抵御粗差的能力并非易事,尤其在觀測(cè)數(shù)據(jù)量很大時(shí) (不妨設(shè)有 1000個(gè)觀測(cè)數(shù)據(jù))。如果這些觀測(cè)數(shù)據(jù)中僅有一個(gè)粗差,粗差在觀測(cè)數(shù)據(jù)中的位置有1000種可能性 (每個(gè)觀測(cè)數(shù)據(jù)都可能包含此粗差)。為此我們需要進(jìn)行 1000次檢驗(yàn),若成功探測(cè)出粗差 1000次,則新理論探測(cè)粗差的效能達(dá)到 100%;若成功 900次,則探測(cè)粗差的效能就是 90%。如果測(cè)量數(shù)據(jù)中含有 2個(gè)粗差,這兩個(gè)粗差在觀測(cè)數(shù)據(jù)中的位置就有 C21000=499500種可能性,為此要進(jìn)行 499500次檢驗(yàn),如果有 3個(gè)粗差呢?要進(jìn)行=166167000次檢驗(yàn),要進(jìn)行如此多次數(shù)的檢驗(yàn),很明顯這是不現(xiàn)實(shí)的[4]。如何設(shè)計(jì)檢驗(yàn)方案,達(dá)到既能科學(xué)評(píng)估新的平差理論和方法效能的目的,又能使檢驗(yàn)次數(shù)有限,在實(shí)際工作中非常必要。這就需要探討粗差的分布,依據(jù)其分布設(shè)計(jì)檢驗(yàn)方案。
值得注意的是,以往討論粗差分布的文章很少[5],且粗差在觀測(cè)數(shù)據(jù)中的位置和大小具有隨意性,基于隨意性基礎(chǔ)上的平差理論評(píng)估是不可靠的[6]。
大量觀測(cè)數(shù)據(jù)的統(tǒng)計(jì)分析表明,粗差在觀測(cè)數(shù)據(jù)中占少數(shù),一般情況,含粗差的觀測(cè)數(shù)據(jù)占數(shù)據(jù)總量的 1%~10%[2,4],大部分觀測(cè)數(shù)據(jù)是正常的。
假定有 n個(gè)觀測(cè)數(shù)據(jù),且觀測(cè)數(shù)據(jù)相互獨(dú)立。每個(gè)觀測(cè)數(shù)據(jù)含有粗差的概率相等,都等于 p(p一般為 1%~10%)。設(shè)這 n個(gè)觀測(cè)數(shù)據(jù)中含有 X個(gè)粗差,顯然 X是隨機(jī)變量,可能取到 0,1,2,…,n這些值。由概率論知識(shí)知:粗差個(gè)數(shù) X服從二項(xiàng)分布,即 X~b(n,p)。
在測(cè)量數(shù)據(jù)處理中可以取 X的平均值 (數(shù)學(xué)期望)E(X)=np作為粗差個(gè)數(shù)的估值。因粗差個(gè)數(shù)是整數(shù),所以一般要對(duì)平均值舍入取整數(shù),假定取整后粗差個(gè)數(shù)是 m=[np]。
用上述方法確定了 n個(gè)觀測(cè)數(shù)據(jù)中含有 m個(gè)粗差,那么這 m個(gè)粗差到底分布在哪些觀測(cè)數(shù)據(jù)上呢?粗差的大小又該如何確定呢?下面將詳細(xì)分析這些問(wèn)題。
首先討論單個(gè)粗差Δ位置 Y和大小 Z的聯(lián)合分布。單個(gè)粗差Δ以等可能分布在每個(gè)觀測(cè)數(shù)據(jù)上,所以粗差的位置 Y分布是離散均勻分布,即
單個(gè)粗差Δ的大小分布可以按如下方法確定。由于測(cè)量誤差通常服從正態(tài)分布,一般認(rèn)為 3σ(σ是標(biāo)準(zhǔn)差)之外的誤差為粗差?;谏鲜稣J(rèn)識(shí),可以構(gòu)造一個(gè)粗差分布,粗差Δ在 (-3σ,+3σ)內(nèi)概率密度為 0,在 3σ之外服從正態(tài)分布,即粗差Δ大小的概率密度函數(shù)為:
由于此分布是在正態(tài)分布基礎(chǔ)上提出的,而且只在兩端有密度,故本文稱之為截尾正態(tài)分布。
綜上所述,可得粗差Δ所在位置 Y和大小 Z的聯(lián)合分布函數(shù) (稱為離散均勻 -截尾正態(tài)分布)為:
很明顯,粗差Δ所在的位置和大小 (Y,Z)的聯(lián)合分布既不是離散型隨機(jī)變量,也不是連續(xù)型隨機(jī)變量,它屬于非離散非連續(xù)型隨機(jī)變量。
如果 n個(gè)觀測(cè)數(shù)據(jù)的權(quán)重不同,則觀測(cè)誤差不同,因此粗差的界定不同。不妨設(shè)第 k個(gè)觀測(cè)數(shù)據(jù)的權(quán)重為 pk,此時(shí)(3)式可改寫為:
根據(jù) (2)式和 (5)式,可以推出粗差Δ所在的位置和大小 (Y,Z)的聯(lián)合概率分布為:
因?yàn)槊總€(gè)觀測(cè)數(shù)據(jù)含有粗差的概率相同,自然認(rèn)為這 m個(gè)粗差Δ(粗差向量)離散均勻分布在 n個(gè)觀測(cè)數(shù)據(jù)上。這相當(dāng)于 m個(gè)學(xué)生離散均勻分布在 n個(gè)不同的座位上,共有可能性,即
由此,知粗差向量所在的位置和大小 (Y,Z)也服從離散均勻——截尾正態(tài)分布,其聯(lián)合分布為:
如果 n個(gè)觀測(cè)數(shù)據(jù)的權(quán)重不同,不妨設(shè)第 k個(gè)觀測(cè)數(shù)據(jù)的權(quán)重為 pk,經(jīng)過(guò)和(6)式類似地推導(dǎo)可得粗差向量Δ所在的位置和大小 (Y,Z)的聯(lián)合概率分布為:
假定觀測(cè)數(shù)據(jù)相互獨(dú)立而且等權(quán),下面將討論利用Matlab[8]軟件產(chǎn)生粗差隨機(jī)數(shù)的方法。
第一步,用函數(shù) round()對(duì)粗差個(gè)數(shù)的數(shù)學(xué)期望舍入取整,得到粗差個(gè)數(shù) 。
第二步,用 unidrnd()產(chǎn)生離散均勻分布函數(shù)的隨機(jī)數(shù)。
第三步,可通過(guò)編程來(lái)產(chǎn)生截尾正態(tài)分布的隨機(jī)數(shù)。首先產(chǎn)生正態(tài)隨機(jī)數(shù),保存正態(tài)隨機(jī)數(shù)的正負(fù)號(hào),然后對(duì)正態(tài)隨機(jī)數(shù)取絕對(duì)值,再向右平移 3σ,最后再把保存的正負(fù)號(hào)再添加到相應(yīng)的經(jīng)過(guò)平移后的隨機(jī)數(shù)上。
具體的程序如下:
clc;clear
disp(’設(shè)定觀測(cè)數(shù)據(jù)n=40,粗差概率p=0.1,取整后粗差數(shù) m=4,標(biāo)準(zhǔn)差 sigma=1,重復(fù)試驗(yàn)次數(shù) k=10′)
n=40;p=0.10;m=round(n*p);sigma=1;k=10;
disp(’產(chǎn)生k組符合要求的m個(gè)離散均勻分布′)
l=unidrnd(n,k,m);
for i=1:k
I=0;
for j=1:m-1
for jj=j+1:m
if l(i,j)= =l(i,jj)
I=I+1;
end
end
end
while I>0
l(i,:)=unidrnd(n,1,m);I=0;
for j=1:m-1
for jj=j+1:m
if l(i,j)= =l(i,jj)
I=I+1;
end
end
end
end
end
l
disp(’產(chǎn)生k組符合要求的m個(gè)截尾正態(tài)隨機(jī)數(shù)(每組一行,每一行 m個(gè)粗差)
t=nor mrnd(0,1,k,m);
%disp(’保留每個(gè)隨機(jī)數(shù)的正負(fù)號(hào)’)
fuhao=sign(t);
%disp(’每個(gè)隨機(jī)數(shù)取絕對(duì)值’)
tt=abs(t)+3*sigma;
%disp(’截尾正態(tài)分布隨機(jī)數(shù)’)
for j=1:m
for i=1:k
x(i,j)=fuhao(i,j)*tt(i,j);
end
end
x
整理程序運(yùn)行結(jié)果后,得到:
設(shè)定觀測(cè)數(shù)據(jù) n=40,粗差概率 p=0.1,取整后粗差數(shù)m=4,標(biāo)準(zhǔn)差 sigma=1.0,重復(fù)試驗(yàn)次數(shù) k=10
表1
(1)粗差個(gè)數(shù)服從二項(xiàng)分布,粗差的位置和大小的聯(lián)合分布為離散均勻——截尾正態(tài)分布,并可用Matlab軟件產(chǎn)生了聯(lián)合分布的隨機(jī)數(shù)。這是科學(xué)評(píng)價(jià)新的平差理論和方法的基礎(chǔ)條件和先決條件。
(2)在產(chǎn)生離散均勻——截尾正態(tài)分布隨機(jī)數(shù)的過(guò)程中有可能在一個(gè)觀測(cè)數(shù)據(jù)上出現(xiàn)兩個(gè)或多個(gè)粗差的情況,這在測(cè)量數(shù)據(jù)處理中是不允許出現(xiàn)的,此時(shí)在程序中應(yīng)舍棄本次試驗(yàn),重新產(chǎn)生符合要求的隨機(jī)數(shù)的方法來(lái)處理。
[1]周江文.經(jīng)典誤差理論與抗差估計(jì)[J].測(cè)繪學(xué)報(bào),1989,18(2):115-120.
[2]周江文,歐吉坤,楊元喜.測(cè)量誤差理論新探 [M].北京:地震出版社,1999:6-8.
[3]哈爾濱工業(yè)大學(xué),上海工業(yè)大學(xué).機(jī)床夾具設(shè)計(jì) [M].上海:上??茖W(xué)技術(shù)出版社,1985.
[4]劉友才,肖繼德.機(jī)床夾具設(shè)計(jì)[M].北京:機(jī)械工業(yè)出版社,1992.
[5]華中工學(xué)院標(biāo)準(zhǔn)化與計(jì)量測(cè)試教研室.互換性與技術(shù)測(cè)量[M].武漢:華中工學(xué)院出版社,1983.
[6]劉登平.機(jī)械制造工藝及機(jī)床夾具設(shè)計(jì) [M].北京:北京理工大學(xué)出版社,2008.
[7]Peiliang Xu.Sign-constrained robust least squares,subjective breakdown point and the effectofweightsof observations on robustness[J].J.of Geod.2005(79):146-159.
[8]Huber P J.Robust Statistics[M].New York:John W iley&Sons,1981.