Kolmogorov-Smirnov检验置信带在Matlab中的实现.docx

Kolmogorov-Smirnov检验置信带在Matlab中的实现张三(北京交通大学土木建筑工程学院，北京100044)摘要：K-S检验是非参数检验的一个十分精确的方法，构造高精度的基于K-S统计量的置信带可以更加直观的判断某随机变量是否服从某一分布。本文首先介绍了K-S检验的统计量，构造了K-S检验置信带，进而采用Matlab软件分析所构造的置信带的准确性，最终得到了较好的结果。关键词：K-S检验；置信带；MatlabImplementation of Kolmogorov-Smirnov Test Confidence Belt in MatlabZHANG San(School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China)Abstract: K-S test is a very accurate method of non-parametric test. High-precision K-S confidence belt based on K-S statistics can be used to tell whether some random variables are submitting to a certain distribution. This paper introduced the K-S test statistics, constructed the K-S test confidence belts, and then used Matlab to analyze the accuracy of confidence belts which were constructed on the basis of different distributions, and ultimately get satisfied results.Key words: K-S test; confident belts; Matlab在可靠性分析过程中，常需要估计一些随机变量服从的分布，再对假设的分布函数进行非参数检验，常用的非参数检验方法有很多种，比如卡方检验、Kolmogorov-Smirnov检验（以下简称K-S检验）等。卡方检验的结果受分组的影响较大，因而容易受到人为因素的影响，故相比而言K-S检验结果更为精确1。利用K-S法检验某种分布类型无非有两种形式，一是已知参数估计分布类型，二是未知参数，先估计参数再检验分布类型2，由于缺乏实际数据，故本文进行已知参数的K-S检验。K-S检验较卡方检验的先进之处在于，K-S检验认为当样本容量n很大时，经验分布Fn(x)与总体分布函数F(x)二者的偏差不应太大，即K-S检验是对样本每一点的检验。因此K-S检验的统计量可记为Dn=max-x0有，P(max|Fnx-F(x)| )2e-2n2令Lx=maxFnx-n,0，Hx=minFnx+n,1，其中n=12nln(2) ，则对所有的Fx和所有的n，都有，PLxFxHx1- (x R)记K-S统计量的1-分位点为Dn1-=n，则由上述结论可得，Lx=maxFnx-Dn1-,0，Hx=minFnx+Dn1-,1，x R，记Li1-=maxin-Dn1-,0，Hi1-=minin+Dn1-,1，可见事件maxFnx-FxDn1- 等价于Li1-FX(i)Hi1-，i=1 , 2 , , n，其中X(i)为次序统计量。假设处在X(i-1)和X(i)之间的x成线性变化，则可得到下面的式子3，Lx=0 (xX1) Li1-+Li+11-Li1-x-XiXi+1-Xi(XixXi+1,i=1,2,n-1)Ln1- (xXn) Hx=H01- (xX1) Hi-11-+Hi1-Hi-11-x-XiXi+1-Xi(XixXi+1,i=1,2,n-1)1 (xXn) 可见得到的置信带是关于x的分段函数。2 K-S检验置信带在Matlab中的实现可靠性分析中常见的分布主要有正态分布、对数正态分布、指数分布、Weibull分布及Gamma分布等，表1给出了这5种分布的密度函数和参数。表1 常用分布密度函数及参数Table 1 Common distribution density function and parameters分布密度函数参数正态分布fx=12e-(x-)222(-x+)，对数正态分布fx=12xe-(lnx-)222(0x0)0 (x0)EX=DX=2Weibull分布fx=ba(xa)b-1e-(xa)b (x0)0 (x0)EX=a(1+1b)DX=a21+2b-(1+1b)2Gamma分布fx=1ba(a)xa-1e-xb(0x+)EX=abDX=ab2由表1可见，上述五种分布的参数都可以用样本的均值和方差表示，故另样本均值mu=10.9，标准差sigma=0.5，在Matlab中生成上述5种分布的随机数各1000个。相关代码如下，function X=fun_rand(n,type,sigma,mu,a,b) switch type %根据type的值进行判断，产生不同分布的随机数 case 1 t=randn(1,n); X=sigma*t+mu; %产生正态分布随机数 case 2 mlg=log(mu/sqrt(1+(sigma/mu)2);siglg=sqrt(log(1+sigma/mu)2); X=lognrnd(mlg,siglg,1,n); %产生对数正态分布随机数 case 3 X=exprnd(mu,1,n); %产生指数分布随机数 case 4 X=wblrnd(a,b,1,n); %产生Weibull分布随机数 case 5 X=gamrnd(a,b,1,n); %产生Gamma分布随机数 endend经验分布函数Fnx很容易得到，将产生的随机数组X进行从小到大的排序，计算Li1-和Hi1-，进而计算置信带的分段函数Lx和Hx。相关代码如下，X=sort(X); %将随机数组X从小到大排序Fn=fun_Fn(n); %由函数fun_Fn生成经验分布函数Dlim=1.36/sqrt(n); %置信水平为5%的分位点Dn1-for i=1:n+1 LX(i)=max(Fn(i)-Dlim,0); HX(i)=min(Fn(i)+Dlim,1); %计算Li1-和Hi1-endLXX,HXX=fun_LH(LX,HX,x,X,n); %由函数fun_LH计算Lx和Hx下面是函数fun_LH的代码，function LXX,HXX=fun_LH(LX,HX,x,X,n)N=length(x);LXX=zeros(1,N);HXX=zeros(1,N);for i=1:N LXX(i)=0*(x(i)X(n); HXX(i)=HX(1)*(x(i)X(n); for j=2:nLXX(i)=LXX(i)+(LX(j)+(LX(j+1)-LX(j)*(x(i)-X(j-1)/(X(j)-X(j-1)*(x(i)X(j-1)&x(i)X(j-1)&x(i)=X(j); endendend3 计算结果利用上述程序和初始数据进行计算，计算结果见图1，2，3，4，5，其中红线均为置信上限，绿线均为置信下限，蓝线为概率累积分布函数图像。图1 正态分布K-S检验置信带(=5%)Fig.1 The K-S test confident belt of normal distribution (=5%)图2 对数正态分布K-S检验置信带(=5%)Fig.2 The K-S test confident belt of lognormal distribution (=5%)图3 指数分布K-S检验置信带(=5%)Fig.3 The K-S test confident belt of exponential distribution (=5%)图4 Weibull分布K-S检验置信带(=5%)Fig.4 The K-S test confident belt of Weibull distribution (=5%)图5 Gamma分布K-S检验置信带(=5%)Fig.5 The K-S test confident belt of gamma distribution (=5%)4 结论和展望由上述计算结果可以看出，K-S检验置信带具有良好的精度。对于实际工程来说，可以将上述的随机数组X由实际的随机变量代替进行计算，如果由实际的随机变量得出的概率累积分布函数值均位于置信上下限之内，说明可以有95%的把握认为随机变量服从某假设分布。参考文献1 葛素娟. 中原地区住宅建筑结构荷载调查与统计模型研究D. 郑州:郑