一个R程序求助！如何知道R模拟出的结果

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

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=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)
}
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)
}
以上程序如何显示结果？如何知道模拟出的结果？命令是什么？谢谢！
本文来自: 人大经济论坛 文献求助专区 版，详细出处参考：https://bbs.pinggu.org/viewthread.php?tid=1112611&page=1&from^^uid=11232

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：R程序 Bootstrap function quantile Bootstra 求助 程序 结果 模拟

1. cnpksim=function(n,iter,lsl,usl){
   
2.    #internal function 1 =========================
   
3.    cnpkf=function(x,lsl,usl){
   
4.       x50=quantile(x,0.5)
   
5.       x99.5=quantile(x,0.995)
   
6.       x0.5=quantile(x,0.005)
   
7.       cnpl=(x50-lsl)/(x50-x0.5)
   
8.       cnpu=(usl-x50)/(x99.5-x50)
   
9.       cnpk=min(cnpl,cnpu); return(cnpk) }
   
10. 
    
11.    #internal function 2 =========================
    
12.    cnpkt=function(lsl,usl){
    
13.       z50=qnorm(0.5); z99.5=qnorm(0.995)
    
14.       z0.5=qnorm(0.005)
    
15.       cnpl=(z50-lsl)/(z50-z0.5)
    
16.       cnpu=(usl-z50)/(z99.5-z50)
    
17.       true.cnpk=min(cnpl,cnpu)
    
18.       return(true.cnpk) }
    
19. 
    
20.    #internal function 3 =========================         
    
21.    cnpkbootpm=function(x,lsl,usl,b){
    
22.       coverage=matrix(0,1,4); n=length(x)
    
23.       cnpk=cnpkt(lsl,usl)
    
24.       T1=matrix(0,b,1); H1=matrix(0,b,1)
    
25.       T1S=matrix(0,b,1); H1S=matrix(0,b,1)
    
26.       cnpkstar=matrix(0,b,1)
    
27.       cnpkstars=matrix(0,b,1)
    
28.       cnpks=matrix(0,b,1); sampleSD=sd(x)
    
29.       sigma=min(sampleSD,IQR(x)/1.349)
    
30.       h=1.587*sigma*n^(-1/3)
    
31.       cnpkhat=cnpkf(x,lsl,usl)
    
32. 
    
33.       for(i in 1:b){
    
34.          xstar=sample(x,n,replace=T)
    
35.          cnpkstar=cnpkf(xstar,lsl,usl)
    
36.          cnpkstar1=matrix(0,b,1)
    
37.          for(j in 1:100){
    
38.             xstar1=sample(xstar,n,replace=T)
    
39.             cnpkstar1[j]=cnpkf(xstar1,lsl,usl)}
    
40.          std=sd(cnpkstar1)
    
41.          T1=(cnpkstar-cnpk)/std
    
42.          H1=cnpkstar-cnpk }
    
43.          cnpks=sort(cnpkstar); se=sd(cnpkstar)
    
44. 
    
45.       # standard bootstrap------------------------
    
46.       BS.CL=cnpkhat-1.96*(se)
    
47.       BS.CU=cnpkhat+1.96*(se)
    
48.       if ((cnpk>=BS.CL)&&(cnpk<=BS.CU))
    
49.          coverage[1]=1; width1=BS.CU-BS.CL
    
50. 
    
51.       # percentile bootstrap----------------------
    
52.       BSP.CL=cnpks[as.integer(b*.05)]
    
53.       BSP.CU=cnpks[as.integer(b*.95)]
    
54.       if ((cnpk>=BSP.CL)&&(cnpk<=BSP.CU))
    
55.          coverage[2]=1; width2=BSP.CU-BSP.CL
    
56. 
    
57.       # calculating CI cnpk- SE(T)--------------------
    
58.       tsort=sort(T1); d1=b*.95; d2=b*.05
    
59.       l1=tsort[as.integer(d1)]
    
60.       u1=tsort[as.integer(d2)]
    
61.       BST.CL=cnpkhat-se*l1; BST.CU=cnpkhat-se*u1
    
62.       if ((cnpk>=BST.CL)&&(cnpk<=BST.CU))
    
63.          coverage[3]=1; width3=BST.CU-BST.CL
    
64. 
    
65.       # calculating CI cnpk-H--------------------------
    
66.       hsort=sort(H1); e1=b*.95; e2=b*.05
    
67.       l2=hsort[as.integer(e1)]
    
68.       u2=hsort[as.integer(e2)]
    
69.       BSH.CL=cnpkhat-l2; BSH.CU=cnpkhat-u2
    
70.       if ((cnpk>=BSH.CL)&&(cnpk<=BSH.CU))
    
71.          coverage[4]=1; width4=BSH.CU-BSH.CL
    
72.       return(coverage) }
    
73. 
    
74.    # ===================================================
    
75.    covmat=matrix(0,iter,4); b=1000
    
76.    for (i in 1:iter) {
    
77.       x=rnorm(n,0,1)
    
78.       covmat[i, ]=cnpkbootpm(x,lsl,usl,b) }
    
79.    return (covmat) }

复制代码

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

原文在附件中

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

> 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) }
>
没有结果

希望大家看看和帮帮忙！

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)
}

#############
n=500
lsl=-2.7822
usl= 2.7822
b=1000
set.seed(123)  
x=rnorm(n,0,1)
cnpk=cnpkt(lsl,usl)
cnpkstar=matrix(0,b,1)
cnpks=matrix(0,b,1)
cnpkhat=cnpkf(x,lsl,usl)    # 1.15713
cnpkhat
for(i in 1:b)
{
xstar=sample(x,n,replace=T)
cnpkstar=cnpkf(xstar,lsl,usl)
}
cnpks=sort(cnpkstar)
se=sd(cnpkstar)     
# standard bootstrap------------------------
BS.CL=cnpkhat-1.96*(se)    # 1.010139
BS.CL
BS.CU=cnpkhat+1.96*(se)    # 1.304122
BS.CU
width1=BS.CU-BS.CL         # 0.2939837   文献 0.2935
width1  
# percentile bootstrap----------------------
#BSP.CL=cnpks[as.integer(b*.05)]
BSP.CL=cnpks[as.integer(b*.025)]  #0.9550955
BSP.CL
#BSP.CU=cnpks[as.integer(b*.95)]
BSP.CU=cnpks[as.integer(b*.975)]  #1.254852
BSP.CU
width2=BSP.CU-BSP.CL              #0.2997566   文献 0.297
width2

###### coverage probabilities
#####Standard
mean(cnpkstar>=BS.CL & cnpkstar<=BS.CU)  #0.928  文献 89.6%

#####Percentile
mean(cnpkstar>=BSP.CL & cnpkstar<=BSP.CU) #0.951  文献 93.3%

#####
MSVAR
   OX_MSVAR.rar (277.23 KB)

2011-10-13 08:40:10 上传

非常感谢epoh大师！

谢谢epoh老师

感谢epoh老师