我稍微整理了一下，不晓得楼主可否提供数据以方便测试呢？

原文在附件中

原文中的程序，是模拟用的，应该不需要程序。
请懂的朋友运行一下，看一下怎么做出原作者的
结果。

> cnpksim=function(n,iter,lsl,usl){
+
+    #internal function 1 =========================
+
+    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) }
+
+
+
+    #internal function 2 =========================
+
+    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) }
+
+
+
+    #internal function 3 =========================         
+
+    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=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=(cnpkstar-cnpk)/std
+
+          H1=cnpkstar-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) }
+
+
+
+    # ===================================================
+
+    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) }
>
没有结果

希望大家看看和帮帮忙！