如何在给定Ranked Set Sampling相关R代码中获取秩j对应的观测值

解决方法 #

你要获取的秩j对应的Xj，是秩集抽样中每组排序后第j位的观测值，你编写的rss函数返回的样本矩阵列名m=1到m=m正好对应秩j=1到j=m的观测值，直接在模拟循环中按列提取即可，完整代码如下：

# 复用你已经编写的两个基础函数
rankedsets<-function(X,m,s=m){
  if(s==m){
    x=sample(X,(m^2),replace=F)
    n=matrix(x,ncol=m,nrow=m,byrow=T)
    ms=matrix(0,ncol=m,nrow=m)
    for (i in 1:m){
      ms[i,]=sort(n[i,])
    }
  }else {
    x=sample(X,(m*s),replace=F)
    n=matrix(x,ncol=m,nrow=s,byrow=T)
    ms=matrix(0,ncol=m,nrow=s)
    for (i in 1:s){
      ms[i,]=sort(n[i,])
    }
  }
  return(ms)
}

rss<-function(X,m,r=1,sets=FALSE){
  rss=numeric()
  set=matrix(0,ncol=m,nrow=(m*r))
  if (is.vector(X)){
    a=0
    for (j in 1:r){
      ms=rankedsets(X,m)
      for (i in 1:(m)){
        set[i+a,]=ms[i,]
        rss[i+a]=ms[i,i]
      }
      a=a+m
    }
    rss=matrix(rss,ncol=m,nrow=r,byrow=T)
    cn=rn=numeric()
    for (i in 1:r){
      rn[i]=paste("r","=",i)
    }
    for (i in 1:m){
      cn[i]=paste("m","=",i)
    }
    rownames(rss)=rn
    colnames(rss)=cn
    if (sets){
      s=list(sets=set,sample=rss)
      return(s)
    } else {
      return(rss)}
  }else stop(" X must be a vector!",call.=F)
}

# 模拟参数设定
sims = 1000
target_j = 2 # 替换为你需要获取的秩j的数值
m = 5 # 秩集大小和你之前的示例保持一致，可按需修改
Xj = rep(NA, sims)
pop_data = rnorm(10000, 1, 3) # 抽样总体，可按需修改分布参数

# 模拟循环直接写在参数定义之后即可
for (k in 1:sims) {
  rss_result = rss(X = pop_data, m = m, r = 1)
  Xj[k] = rss_result[, target_j] # 提取对应秩j的观测值
}

# 运行完成后Xj就存储了1000次模拟得到的对应秩j的观测值，可通过以下代码查看验证
head(Xj)

补充说明 #

• 如果你的r参数设置大于1，需要取多次循环的所有对应j的结果，把Xj改为矩阵存储即可，rss_result[, target_j]会返回长度为r的向量，直接按行/列赋值到Xj矩阵对应位置即可
• 如果需要同时存储多个j的结果，把Xj改为多列矩阵，每次循环提取多列分别赋值即可

内容的提问来源于stack exchange，提问作者Felix