实验室测量不确定度和统计分析新进展
通风柜
中国计量科学研究院
刘智敏
刘智敏
VIM3-ISO.IEC Guide99
.测量结果
量值和不确定度
.测量不确定度
基于所用信息,表征赋予被测量量值之分散性的非负参数
定义不确定度
不确定度
2. GUM
Y=f(X1, …,XN)
y=f(x1, …,xN)
u2(y)= Σ(∂f/ ∂xi) 2 u2(xi) + 相关项
= Σ ui 2 + 相关项
U=k u(y) , 包含区间[y-U, y+U]
k=tp (ν)
ν=u4(y)/ [∑ui4/νi]
简化计算
p=0.95, k=2+ 2.4/ν
MCS
概述
蒙特卡洛(Monte Carlo)方法是随机统计模拟解问题的方法
随机数可调用 MATLAB
3. 应用
(1) 求π
[-1,1]*[-1,1]作均匀随机数 x,y.当x2+y2≤1,
(x,y)在圆内.取M点,圆内a点
圆内面积:正方形面积
= π:4=a:M
π=4*a:M
M=M0; rand('seed',0); a=0;
for i=1:M
x=-1+2*rand;
y=-1+2*rand;
r2=x*x+y*y;
if r2<=1
a=a+1;
else
end
end
PPI=4*a/M
M=104, π=3.15
M=106, π=3.145
M=108, π=3.1416
(2)求N(0,1)在[-1,1]概率P
P=
f(x)=(2 π) - 1/2exp(-x2/2)
[-1,1]*[0,1]作均匀随机数 x,y.
.取M点,当y≤f计a点
P:2=a:M
P=2*a:M
M=M0 ; rand('seed',0)
a=0;
for i=1:M
x=-1+2*rand;
y=rand;
f=exp(-x*x/2)/sqrt(2*pi);
if y<=f
a=a+1;
else
end
end
P=2*a/M
M=104, P=0.6856
M=106, P=0.6824
M=108, P=0.6827
二.不确定度评定原理
1.基础
(1)内容
建模 Y=f(X1, …,XN)
分布传播
给出
y, u(y), Y的包含区间
(2)途径: A:解析法; G: GUM; M: MCM
2. GUM
Y=f(X1, …,XN)
y=f(x1, …,xN)+Σ(∂f/ ∂xi)(Xi-xi)+…
一阶
u2(y)= Σ(∂f/ ∂xi) 2 u2(xi) + 相关项
U=tp(ν) u(y) , 包含区间[y-U, y+U]
三.MCM作不确定度分布传播
1.特点
非线性
非对称
2.步骤
Y=f(X1, …,XN)=f(X)
(1)对Xi 抽样,得 Y=f(X)
(2)作M次试行
形成 Y的分布G
(3)用 G 得 y, u(y),
最短包含区间[ylow, yhigh]
最短包含区间 [IL ,IR]
N(μ,σ2)概率对称包含区间
[μ-U,μ+U]=[μ-kσ,μ+kσ]
非对称分布概率对称包含区间
[μ-U L,μ+UR], U L ≠ UR
包含概率为p的包含区间中最短者[IL ,IR]
f(IL )=f(IR)
3.输入量的分布
(1)矩形分布 R(a,b)
g (ξ)=1/(b-a) ,a≤ ξ ≤b
=0 ,其它
E=(a+b)/2, V=(b-a)2/12
ξ=a+(b-a)r
(2)限不确的矩形分布 CTrap(a,b,d)
矩形分布的限为a ±d, b ±d
g(ξ)=ln[(w+d)/(x-ξ)]/[4d] ,a-d≤ξ≤a+d
=ln[(w+d)/(w-d)]/[4d] ,a+d≤ξ≤b-d
=ln[(w+d)/(ξ-x)]/[4d] ,b-d≤ξ≤b+d
=0 ,其它
x=(a+b)/2,w=(b-a)/2
E=(a+b)/2, V=(b-a)2/12+d2/9
as=(a-d)+2dr1,,bs=(a+b)-as,ξ=as+(b s–as )r2
梯形分布 Trap(a,b,β)
独立R(a1 ,b1), R(a2 ,b2)的和
g(ξ)=(ξ-x+λ2)/(λ22 -λ1 2 ), x- λ2 ≤ξ≤ x- λ1
=1/(λ1 + λ2) , x- λ1 ≤ξ≤ x+ λ1
=(x+λ2 -ξ)/(λ22 -λ1 2 ), x+ λ1 ≤ξ≤ x+ λ2
=0 ,其它
x=(a+b)/2, a=a1+a2,b=b1+b2;
λ1= (b1+a2 -a1- b2)/2, λ2= (b -a)/2, β= λ1 /λ2
E=(a+b)/2, V=(b-a)2 (1+ β2) /24
ξ=a+(b–a )[(1+ β)r1 +(1 - β)r2]/2
(4)三角分布T(a,b)= Trap(a,b,0)
g(ξ)=(ξ-a)/w2, a ≤ξ≤ x
=(b -ξ)/w2, x ≤ξ≤b
=0 ,其它
x=(a+b)/2, w= λ2= (b -a)/2,
E=(a+b)/2, V=(b-a)2 /24
ξ=a+(b–a )[r1 +r2]/2
(5)反正弦分布 As(a,b)
当Φ为R(0,2π)则
[a+b+(b-a)sinΦ]/2为As(a,b)
g(ξ)=(2/π)[(b-a)2-(2ξ-a-b)2]-1/2, a ≤ξ≤ b
=0 ,其它
E=(a+b)/2, V=(b-a)2 /8
ξ=[b+a+(b–a )sin(2 πr)]/2
当 (X-μ)/σ为t(v), 称X为tν(μ,σ2)
g(ξ)=(.5v+.5)*
[1+(ξ- μ)2 /v σ2)-.5v-.5/[(πv).5 Γ(.5v) σ]
E= μ, V= vσ2/(v-2)
ξ=μ+σt
1) 平均
xk ,k=1,…,n 独立高斯,平均 值为tν(μ,σ2): v=n-1, μ=Σ xk /n,
σ2= [Σ xk-( Σ xk/n)]2/ [n(n-1)]
2)证书:X估计x, Up(…kp,…v)
X为tν(x, (Up/kp )2)
(9)投影分布
g(ξ)=[A 2(2 ξ- ξ2)]- 0.5
E=(1-cos A)/3, σ=3 (1-cos A)/10
ξ=1- cos(Ar)
4. MCM执行
由 yr , r=1,2,…, M 得Y 分布
样本平均 ˜y =Σ yr / M
样本方差 u2(˜y )= Σ (yr - ˜y)2 /( M-1)
M=104 --- 106
两结果差应小于 (有效末位的 ½)
四.例
1.正态分布
Y=X1+X2+X3+X4, Xi 独立N(0,1)
A: Y~N(0,22)
G: y=0.00, u(y)=2.00, U=1.96u(y)=3.92(p=0.95)
[y-U,y+U]=[-3.92,3.92]
M: M=105, y=0.00, u(y)=2.00,
[y low,y high]=[-3.92,3.92]
M=M0
Randn(‘seed’,0)
A=randn(1,M); B=randn(1,M);
C=randn(1,M);D=randn(1,M);
Y=A+B+C+D;
[mean(Y),std(Y)]
[prctile(Y,2.5),prctile(Y,97.5)]
2.端度规校准
d=L(1+αθ)- L s(1+α sθ s)
L=L s +d - L s(θδα + α sδθ)
δL=L s +D+d1+d2 - L s[(θ0+Δ)δα
+αsδθ)]-Lnom
Ls~tν(μ,σ2);μ=50000623nm,σ=25nm, v=8
D~tν(μ,σ2);μ=215nm,σ=6nm, v=24
d1 ~tν(μ,σ2); μ=0nm,σ=4nm, v=5
d2 ~tν(μ,σ2); μ=0nm,σ=7nm, v=8
αs ~R(a,b); a=9.5*10-6 ºC-1,, b=13.5*10-6 ºC-1,,
θ0 ~ N(μ,σ2);μ=-0.1 ºC ,σ= 0.2 ºC
Δ ~As(a,b); a= -0.5 ºC ,b= 0.5 ºC
δα ~CTrap(a,b,d);a= -1.0 *10-6 ºC-1,b=1.0 *10-6 ºC-1
d=0.1 *10-6 ºC-1
δθ~CTrap(a,b,d);a= -0.050 ºC,b=0.050 ºC
d=0.025 ºC
M=M0; n=n0
rand('seed',0);randn('seed',0)
for i=1:M
Ls=50000623+25*trnd(18); D=215+6*trnd(24);
d1=4*trnd(5); d2=7*trnd(8);
as1=-0.0000012+0.0000004*rand;dA=as1-2*as1*rand;
c0=normrnd(-0.1,0.2); De=-0.5+rand;
As=0.0000095+0.000004*rand;
as2=-0.1+0.1*rand; dc=as2-2*as2*rand;
dL(i)=Ls+D+d1+d2-Ls*(dA*(c0+De)+As*dc)-50000000;
end
[mean(dL),std(dL)]
for i=1:n %以下求最短0.95 CI
CR(i)=prctile(dL,99+i/n); CL(i)=prctile(dL,i/n);
C(i)=CR(i)-CL(i);
end
a=min(C);
for i=1:n
if C(i) == a
i;
C(i)
[CR(i),CL(i)]
else
end
end
G:δL=838nm, u(δL)=32nm,
最短0.99CI=[745nm,931nm]
M: δL=838nm, u(δL)=36nm,
最短0.99CI=[745nm,932nm]
当相容
Y=KCRV
=xref
xi的等效度:DoE(di , U(di ))
di = xi –xref
U(di )=2{u2 (xi)-u2 (xref)}1/2
xi ,xJ的等效度:DoE(diJ , U(diJ ))
diJ = xi –xJ
U(di J)=2{u2 (xi)+u2 (xJ)}1/2
稳健统计
对X等精度独立测量得 x1, x2,..., xn,
次序量 x（1）≤ x（2）≤... ≤ x（n)
a)分位数
中位数
med x =x ( (n+1)/2 ) , n为奇
={x ( n/2 ) +x ( n/2 +1 )}/2 , n为偶
,下四分位数 {x}1/4=x(n/4+3/4)
上四分位数 {x} 3/4=x(3n/4+1/4)
{x}: 0.8, 2.1, 2.2 , 2.4, 2.6
med x =2.2
{x}1/4 =x(5/4+3/4) =x(2)=2.1
{x}3/4= x(3*5/4+1/4)=x(4)=2.4
a 位置 中位数
切尾平均
{x}; 得分 9.1, 9.1 , 9.2, 9.3, 9.8
xT= (9.1+9.2+9.3)/3 =9.2
b 分散
s=NIQR
=.7413(x3/4-x1/4)
b)中位
s= med| x – med x|/0.6745
=1.483 med| x – med x|
)
C 并尾法
a)初值
x*=med x
s*=1.483 med| x – x*|
b)迭代新值
Φ=1.5 s*
X(i)*= x*-Φ, X(i)< x*-Φ
= x*+Φ, X(i)> x*+Φ
= X(i),
x*=Σ X(i)* /n
s*=1.13{Σ (X(i)* - x*)2 /(n-1)}1/2
|
迭代
|
0次
|
1次
|
2次
|
|
Φ
|
|
1.424
|
1.478
|
|
x*-Φ
|
|
18.876
|
18.909
|
|
x*+Φ
|
|
21.724
|
21.865
|
|
X(1) *
|
17.570
|
18.876
|
18.909
|
|
X(2) *
|
19.500
|
19.500
|
19.500
|
|
X(3) *
|
20.100
|
10.100
|
10.100
|
|
X(4) *
|
20.155
|
20.155
|
20.155
|
|
X(5) *
|
20.300
|
20.300
|
20.300
|
|
X(6) *
|
20.705
|
20.705
|
20.705
|
|
X(7) *
|
20.040
|
20.940
|
20.940
|
|
X(8) *
|
21.185
|
21.185
|
21.185
|
|
X(9) *
|
24.140
|
20.724
|
21.865
|
|
X *
|
20.300
|
20.387
|
20.407
|
|
s *
|
0.949
|
0.985
|
1.009
|
能力验证
[5]ISO13528
1)均匀性稳定性
ss ≤0.3
│ x..-y.. │≤0.3
2)性能统计
Z=(x-X)/σ^
En=(x-X)/(Ulab-Uref)0.5
3) Youden 椭圆
ISO/TC69-2008
上一篇:实验室知识产权管理
下一篇:化学实验室 标配仪器设备