粗差分布及其隨機數產生

李國重,夏云青,焦慧平

(1.信息工程大學理學院,鄭州 450001;2.中州大學信息工程學院,鄭州 450044)

0.引言

在對加工工件的測量數據處理中,經常需要對新的測量平差理論或方法進行評估,其中一項重要工作是檢驗新理論或新方法抵御粗差的能力 (粗差是指離群的大誤差[1,2],由失誤等引起,它實際不可避免)。如果新理論或新方法能夠有效地抵御粗差的不良影響,那么認為新理論或新方法有效,從而在測量數據處理中可以大膽使用;否則,我們懷疑新理論或新方法的效能,需謹慎使用或者棄之不用[3]。理論上講,一個性能優良的平差理論能夠抵御觀測數據的任何粗差,無論這些粗差分布在哪些觀測數據中以及粗差的量級多大。然而,檢驗平差理論抵御粗差的能力并非易事,尤其在觀測數據量很大時 (不妨設有 1000個觀測數據)。如果這些觀測數據中僅有一個粗差,粗差在觀測數據中的位置有1000種可能性 (每個觀測數據都可能包含此粗差)。為此我們需要進行 1000次檢驗,若成功探測出粗差 1000次,則新理論探測粗差的效能達到 100%;若成功 900次,則探測粗差的效能就是 90%。如果測量數據中含有 2個粗差,這兩個粗差在觀測數據中的位置就有 C21000=499500種可能性,為此要進行 499500次檢驗,如果有 3個粗差呢?要進行=166167000次檢驗,要進行如此多次數的檢驗,很明顯這是不現實的[4]。如何設計檢驗方案,達到既能科學評估新的平差理論和方法效能的目的,又能使檢驗次數有限,在實際工作中非常必要。這就需要探討粗差的分布,依據其分布設計檢驗方案。

值得注意的是,以往討論粗差分布的文章很少[5],且粗差在觀測數據中的位置和大小具有隨意性,基于隨意性基礎上的平差理論評估是不可靠的[6]。

1.粗差個數的分布

大量觀測數據的統計分析表明,粗差在觀測數據中占少數,一般情況,含粗差的觀測數據占數據總量的 1%～10%[2,4],大部分觀測數據是正常的。

假定有 n個觀測數據,且觀測數據相互獨立。每個觀測數據含有粗差的概率相等,都等于 p(p一般為 1%～10%)。設這 n個觀測數據中含有 X個粗差,顯然 X是隨機變量,可能取到 0,1,2,…,n這些值。由概率論知識知:粗差個數 X服從二項分布,即 X～b(n,p)。

在測量數據處理中可以取 X的平均值 (數學期望)E(X)=np作為粗差個數的估值。因粗差個數是整數,所以一般要對平均值舍入取整數,假定取整后粗差個數是 m=[np]。

2.粗差位置和大小的聯合分布

用上述方法確定了 n個觀測數據中含有 m個粗差,那么這 m個粗差到底分布在哪些觀測數據上呢?粗差的大小又該如何確定呢?下面將詳細分析這些問題。

2.1 單個粗差位置和大小的聯合分布

首先討論單個粗差Δ位置 Y和大小 Z的聯合分布。單個粗差Δ以等可能分布在每個觀測數據上,所以粗差的位置 Y分布是離散均勻分布,即

單個粗差Δ的大小分布可以按如下方法確定。由于測量誤差通常服從正態分布,一般認為 3σ(σ是標準差)之外的誤差為粗差。基于上述認識,可以構造一個粗差分布,粗差Δ在 (-3σ,+3σ)內概率密度為 0,在 3σ之外服從正態分布,即粗差Δ大小的概率密度函數為:

由于此分布是在正態分布基礎上提出的,而且只在兩端有密度,故本文稱之為截尾正態分布。

綜上所述,可得粗差Δ所在位置 Y和大小 Z的聯合分布函數 (稱為離散均勻 -截尾正態分布)為:

很明顯,粗差Δ所在的位置和大小 (Y,Z)的聯合分布既不是離散型隨機變量,也不是連續型隨機變量,它屬于非離散非連續型隨機變量。

如果 n個觀測數據的權重不同,則觀測誤差不同,因此粗差的界定不同。不妨設第 k個觀測數據的權重為 pk,此時(3)式可改寫為:

根據 (2)式和 (5)式,可以推出粗差Δ所在的位置和大小 (Y,Z)的聯合概率分布為:

2.2 多個粗差位置和大小的聯合分布

因為每個觀測數據含有粗差的概率相同,自然認為這 m個粗差Δ(粗差向量)離散均勻分布在 n個觀測數據上。這相當于 m個學生離散均勻分布在 n個不同的座位上,共有可能性,即

由此,知粗差向量所在的位置和大小 (Y,Z)也服從離散均勻——截尾正態分布,其聯合分布為:

如果 n個觀測數據的權重不同,不妨設第 k個觀測數據的權重為 pk,經過和(6)式類似地推導可得粗差向量Δ所在的位置和大小 (Y,Z)的聯合概率分布為:

3.粗差分布的隨機產生方法

假定觀測數據相互獨立而且等權,下面將討論利用Matlab[8]軟件產生粗差隨機數的方法。

第一步,用函數 round()對粗差個數的數學期望舍入取整,得到粗差個數 。

第二步,用 unidrnd()產生離散均勻分布函數的隨機數。

第三步,可通過編程來產生截尾正態分布的隨機數。首先產生正態隨機數,保存正態隨機數的正負號,然后對正態隨機數取絕對值,再向右平移 3σ,最后再把保存的正負號再添加到相應的經過平移后的隨機數上。

具體的程序如下:

clc;clear

disp(’設定觀測數據n=40,粗差概率p=0.1,取整后粗差數 m=4,標準差 sigma=1,重復試驗次數 k=10′)

n=40;p=0.10;m=round(n＊p);sigma=1;k=10;

disp(’產生k組符合要求的m個離散均勻分布′)

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(’產生k組符合要求的m個截尾正態隨機數(每組一行,每一行 m個粗差)

t=nor mrnd(0,1,k,m);

%disp(’保留每個隨機數的正負號’)

fuhao=sign(t);

%disp(’每個隨機數取絕對值’)

tt=abs(t)+3＊sigma;

%disp(’截尾正態分布隨機數’)

for j=1:m

for i=1:k

x(i,j)=fuhao(i,j)＊tt(i,j);

end

end

x

整理程序運行結果后,得到:

設定觀測數據 n=40,粗差概率 p=0.1,取整后粗差數m=4,標準差 sigma=1.0,重復試驗次數 k=10

表1

4.結論

(1)粗差個數服從二項分布,粗差的位置和大小的聯合分布為離散均勻——截尾正態分布,并可用Matlab軟件產生了聯合分布的隨機數。這是科學評價新的平差理論和方法的基礎條件和先決條件。

(2)在產生離散均勻——截尾正態分布隨機數的過程中有可能在一個觀測數據上出現兩個或多個粗差的情況,這在測量數據處理中是不允許出現的,此時在程序中應舍棄本次試驗,重新產生符合要求的隨機數的方法來處理。

[1]周江文.經典誤差理論與抗差估計[J].測繪學報,1989,18(2):115-120.

[2]周江文,歐吉坤,楊元喜.測量誤差理論新探 [M].北京:地震出版社,1999:6-8.

[3]哈爾濱工業大學,上海工業大學.機床夾具設計 [M].上海:上海科學技術出版社,1985.

[4]劉友才,肖繼德.機床夾具設計[M].北京:機械工業出版社,1992.

[5]華中工學院標準化與計量測試教研室.互換性與技術測量[M].武漢:華中工學院出版社,1983.

[6]劉登平.機械制造工藝及機床夾具設計 [M].北京:北京理工大學出版社,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.