cnpkf=function(x,lsl,usl){
x50=quantile(x,0.5)
x99.5=quantile(x,0.995)
x0.5=quantile(x,0.005)
cnpl=(x50-lsl)/(x50-x0.5)
cnpu=(usl-x50)/(x99.5-x50)
cnpk=min(cnpl,cnpu)
return(cnpk)
}
cnpkt=function(lsl,usl){
z50=qnorm(0.5)
z99.5=qnorm(0.995)
z0.5=qnorm(0.005)
cnpl=(z50-lsl)/(z50-z0.5)
cnpu=(usl-z50)/(z99.5-z50)
true.cnpk=min(cnpl,cnpu)
return(true.cnpk)
}
cnpkbootpm=function(x,lsl,usl,b)
{
coverage=matrix(0,1,4)
n=length(x)
cnpk=cnpkt(lsl,usl)
T1=matrix(0,b,1)
H1=matrix(0,b,1)
T1S=matrix(0,b,1)
H1S=matrix(0,b,1)
cnpkstar=matrix(0,b,1)
cnpkstars=matrix(0,b,1)
cnpks=matrix(0,b,1)
sampleSD=sd(x)
sigma=min(sampleSD,IQR(x)/1.349)
h=1.587*sigma*n^(-1/3)
cnpkhat=cnpkf(x,lsl,usl)
for(i in 1:b)
{
xstar=sample(x,n,replace=T)
cnpkstar[i]=cnpkf(xstar,lsl,usl)
cnpkstar1=matrix(0,b,1)
for(j in 1:100)
{
xstar1=sample(xstar,n,replace=T)
cnpkstar1[j]=cnpkf(xstar1,lsl,usl)
}
std=sd(cnpkstar1)
T1[i]=(cnpkstar[i]-cnpk)/std
H1[i]=cnpkstar[i]-cnpk
}
cnpks=sort(cnpkstar)
se=sd(cnpkstar)
# standard bootstrap------------------------
BS.CL=cnpkhat-1.96*(se)
BS.CU=cnpkhat+1.96*(se)
if ((cnpk>=BS.CL)&&(cnpk<=BS.CU))
coverage[1]=1
width1=BS.CU-BS.CL
# percentile bootstrap----------------------
BSP.CL=cnpks[as.integer(b*.05)]
BSP.CU=cnpks[as.integer(b*.95)]
if ((cnpk>=BSP.CL)&&(cnpk<=BSP.CU))
coverage[2]=1
width2=BSP.CU-BSP.CL
# calculating CI cnpk- SE(T)--------------------
tsort=sort(T1)
d1=b*.95
d2=b*.05
l1=tsort[as.integer(d1)]
u1=tsort[as.integer(d2)]
BST.CL=cnpkhat-se*l1
BST.CU=cnpkhat-se*u1
if ((cnpk>=BST.CL)&&(cnpk<=BST.CU))
coverage[3]=1
width3=BST.CU-BST.CL
# calculating CI cnpk-H--------------------------
hsort=sort(H1)
e1=b*.95
e2=b*.05
l2=hsort[as.integer(e1)]
u2=hsort[as.integer(e2)]
BSH.CL=cnpkhat-l2
BSH.CU=cnpkhat-u2
if ((cnpk>=BSH.CL)&&(cnpk<=BSH.CU))
coverage[4]=1
width4=BSH.CU-BSH.CL
return(coverage)
}
cnpksim=function(n,iter,lsl,usl)
{
covmat=matrix(0,iter,4)
b=1000
for (i in 1:iter)
{
x=rnorm(n,0,1)
covmat[i, ]=cnpkbootpm(x,lsl,usl,b)
}
return (covmat)
}
以上程序如何显示结果？