if语句——在R中生成两个具有一定相关性和特定条件的级数

Question

我想在R中生成两个大小为100的数据序列，其中一个是来自Exp(mean=1)分布的缓解时间tr，另一个是来自Exp(mean=2.5)分布的生存时间t。我希望它们呈负相关(例如，相关性为-0.5)。但同时我希望R避免t[I]对于数据点 I 小于tr[I]的值，因为生存时间应该大于缓解时间。我已经能够使用以下代码在两个变量之间产生一些相关性(尽管相关性没有完全重现):

rho <- -0.5
mu <- rep(0,2)
Sigma <- matrix(rho, nrow=2, ncol=2) + diag(2)*(1 - rho)
library(MASS)
rawvars <- mvrnorm(100, mu=mu, Sigma=Sigma)
pvars <- pnorm(rawvars)
    
tr<-rep(0,100)
for(i in 1:100){
tr[i] <- qexp(pvars[,1][i], 1/1)
    }
    
t<-rep(0,100)
for(i in 1:100){
repeat { 
    t[i] <- qexp(pvars[,2][i], 1/2) 
    if (t[i]>tr[i]) break
}
}
    
cor(tr,t)
sum(tr>t) # shows number of invalid cases

但是如何有效地诱导条件，使R只生成t大于相应的tr的值?

此外，是否有更好的方法(更快的方法)在R中完成整个事情?

Answer 1

这里的问题是qexp是分位数函数，并且对于相同的概率pvars[,2][i]将返回相同的值。因此，当pvars[i,]中的任何一个与t[i]<=tr[i]相同时，您的代码很容易进入无限循环。为了避免这种情况，您必须为每个不符合条件的t[i], tr[i]对重新生成rawvars。此外，不需要对pvars进行循环，因为qexp和运算符>都是矢量化的。下面的代码可以满足您的要求:

rho <- -0.5
mu <- rep(0,2)
Sigma <- matrix(rho, nrow=2, ncol=2) + diag(2)*(1 - rho)
library(MASS)
set.seed(1)  ## so that results are repeatable
compute.tr.t <- function(n, paccept) {
  n <- round(n / paccept)
  rawvars <- mvrnorm(n, mu=mu, Sigma=Sigma)
  pvars <- pnorm(rawvars)
  tr <- qexp(pvars[,1], 1/1)
  t <- qexp(pvars[,2], 1/2)
  keep <- which(t > tr)
  return(data.frame(t=t[keep],tr=tr[keep]))
}
n <- 10000  ## generating 10000 instead of 100, this can now be large
paccept <- 1
res <- data.frame()
while (n > 0) {
  new.res <- compute.tr.t(n, paccept)
  res <- rbind(res, new.res)
  paccept <- nrow(new.res) / n
  n <- n - nrow(res)
}

指出:

1. 函数compute.tr.t借用了这里的拒绝采样技术。它的输入参数是我们想要的请求样本数量和期望接受的概率。用这个:

   • 它为tr和t生成n = n / paccept指数变量，因为你这样做是为了说明接受的概率
   • 只保留满足t > tr条件的。

   compute.tr.t返回的可能小于请求的n样本。然后，我们可以使用这些信息来计算我们还需要多少个样本，以及更新后的期望接受概率是多少。

2. 我们在while循环中生成满足我们条件的样本。在这个循环中:

   • 我们调用compute.tr.t，要求生成的样品数量和期望的接受率。最初，这些将分别设置为我们想要的总样本数和1。
compute.tr.t的结果被附加到res的结果数据帧中。
3. 更新接受概率就是返回的样本数与请求的样本数之比。
4. 更新请求的样本数量只是我们需要从我们想要的总数中增加多少。
5. 当下一个请求的样本数量小于或等于0(即我们有足够的样本)时，我们停止。
结果数据帧可能包含比我们想要的样本总数更多的样本。

运行这段代码，我们得到:

print(cor(res$tr,res$t))
[1] -0.09128498
print(sum(res$tr>res$t)) # shows number of invalid cases
##[1] 0

我们注意到反相关性明显弱于预期。这是由于你的身体状况。如果我们通过将compute.tr.t修改为

来删除此条件:

compute.tr.t <- function(n, paccept) {
  n <- round(n / paccept)
  rawvars <- mvrnorm(n, mu=mu, Sigma=Sigma)
  pvars <- pnorm(rawvars)
  tr <- qexp(pvars[,1], 1/1)
  t <- qexp(pvars[,2], 1/2)
  return(data.frame(t=t,tr=tr))
}

则得到:

print(cor(res$tr,res$t))
##[1] -0.3814602
print(sum(res$tr>res$t)) # shows number of invalid cases
##[1] 3676

这种相关性现在更加合理了，但是无效病例的数量还是很大。