我希望从均值为0和标准差为-1的正态分布中采样50000个值。但是我想将这些值限制在[-3,3]之间。我已经编写了相应的代码，但不确定它是否最有效？希望能够得到一些建议。
lower <- -3
upper <- 3
x_norm<-rnorm(75000,0,1)
x_norm<-x_norm[which(x_norm >=lower & x_norm<=upper)]
repeat{
x_norm<-c(x_norm, rnorm(10000,0,1))
x_norm<-x_norm[which(x_norm >=lower & x_norm<=upper)]
if(length(x_norm) >= 50000){break}
}
x_norm<-x_norm[1:50000]
像你的代码一样肯定可以工作，但你高估了需要多少个值。考虑到这是一个已知分布和相当大量的样本，你知道会出现多少少于或多于3个值。
(1-pnorm(3))*2 * 50000
[1] 134.9898
因此，在抽取50000个数的情况下，你可能只有约135个数超出范围。因此，很容易再多抽几个数，但仍不会显著增加数量并对其进行修剪。只需取小于或大于3的前50,500个数中的前50,000个。
x <- rnorm(50500)
x <- x[x < 3 & x > -3]
x <- x[1:50000]
我将前两行代码运行了40,000次，每次都返回长度大于50,000。一个小的布尔检查可以保证它始终如此。
x <- 1
while (length(x) < 50000){
x <- rnorm(50500)
x <- x[x < 3 & x > -3]}
x <- x[1:50000]
对我来说，这个方法几乎在6毫秒内执行了100%的时间。这是一种在R中实现非常快速、易读且不需要附加组件的简单方法。
/tmp/rnormClamp.cpp：#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector rnormClamp(int N, int mi, int ma) {
NumericVector X = rnorm(N, 0, 1);
return clamp(mi, X, ma);
}
/*** R
system.time(X <- rnormClamp(50000, -3, 3))
summary(X)
*/
使用sourceCpp()（同时也可以使用Rcpp）进行构建和运行。在我的电脑上，实际的绘制和夹紧大约需要4毫秒:
R> sourceCpp("/tmp/rnormClamp.cpp")
R> system.time(X <- rnormClamp(50000, -3, 3))
user system elapsed
0.004 0.000 0.004
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.00000 -0.67300 -0.00528 0.00122 0.68500 3.00000
R>
clamp() 糖函数在 Romain 的 SO 回答中被介绍，并且指出您需要 Rcpp 的版本为 0.10.2。
编辑：根据 Ben 的提示，我似乎误解了。这是 C++ 和 R 的混合内容：
// [[Rcpp::export]]
List rnormSelect(int N, int mi, int ma) {
RNGScope scope;
int N2 = N * 1.25;
NumericVector X = rnorm(N2, 0, 1);
LogicalVector ind = (X < mi) | (X > ma);
return List::create(X, ind);
}
R> system.time({ Z <- rnormSelect(50000, -3, 3);
+ X <- Z[[1]][ ! Z[[2]] ]; X <- X[1:50000]})
user system elapsed
0.008 0.000 0.009
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.00000 -0.68200 -0.00066 -0.00276 0.66800 3.00000
R>
我将重新审视逻辑索引和行子集，这部分内容我需要查阅。也许明天吧。但9毫秒还算不错 :)
编辑2: 看起来我们确实没有逻辑索引。我们需要添加它。此版本是“手动”执行的，但比从R索引略快。
// [[Rcpp::export]]
NumericVector rnormSelect2(int N, int mi, int ma) {
RNGScope scope;
int N2 = N * 1.25;
NumericVector X = rnorm(N2, 0, 1);
LogicalVector ind = (X >= mi) & (X <= ma);
NumericVector Y(N);
int k=0;
for (int i=0; i<N2 & k<N; i++) {
if (ind[i]) Y(k++) = X(i);
}
return Y;
}
输出结果：
R> system.time(X <- rnormSelect2(50000, -3, 3))
user system elapsed
0.004 0.000 0.007
R> summary(X)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.99000 -0.66900 -0.00258 0.00223 0.66700 2.99000
R> length(X)
[1] 50000
R>
John和Dirk给出了拒绝抽样的好例子，对于给定的问题应该很合适。但是为了提供另一种方法，当您拥有累积分布函数及其反函数（或合理的近似函数）时，您可以从均匀分布中生成数据并进行转换：
x <- qnorm( runif(50000, pnorm(-3), pnorm(3)) )
range(x)
hist(x)
对于所提出的问题，我不指望这种方法比拒绝抽样方法更好（如果有更好的话），但如果您想要从一个截断的标准正态分布中生成介于2和3之间的数据，则这种方法可能会更加高效。它依赖于累积分布函数及其反函数（在本例中为pnorm和qnorm），因此对于没有这些易于获得的分布的拒绝抽样来说，并不像那么简单。
truncnorm::rtruncnorm()感兴趣。 - hplieninger