10

我想从平均值 = 0 和 sd -1 的正态分布中抽取 50,000 个值。但我想将值限制为 [-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]
4

3 回答 3

15

像你的代码这样的东西肯定会起作用，但你大大高估了你需要多少值。鉴于它是一个已知分布和相当多的样本，您知道有多少样本会出现多于或少于 3 个。

(1-pnorm(3))*2 * 50000
[1] 134.9898

因此，鉴于您在 50,000 的平局中可能只得到约 135 超出范围，因此很容易再抽出几个，但仍然不是一个非常大的数字并对其进行修剪。只需取 50,500 个中小于或大于 3 的前 50,000 个。

x <- rnorm(50500)
x <- x[x < 3 & x > -3]
x <- x[1:50000]

我运行了前 2 行 40,000 次，它每次返回的长度都大于 50000。一个小的布尔检查可以保证它总是这样。

x <- 1
while (length(x) < 50000){
    x <- rnorm(50500)
    x <- x[x < 3 & x > -3]}
x <- x[1:50000]

对我来说，这在 6 毫秒内执行了几乎 100% 的时间。这是在 R 中执行此操作的一种简单方法，执行速度非常快，易于阅读，并且不需要附加组件。

于 2012-12-25T22:57:55.097 回答
11

如果您真的关心效率，那么这段简短的 Rcpp 代码将很难被击败。将以下内容存储在文件中，例如/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>
于 2012-12-25T21:51:56.857 回答
10

John 和 Dirk 给出了拒绝抽样的好例子，这对于给定的问题应该没问题。但要给出另一种方法，当您拥有累积分布函数及其逆函数（或其合理近似值）时，您可以从均匀分布生成数据并进行变换：

x <- qnorm( runif(50000, pnorm(-3), pnorm(3)) )
range(x)
hist(x)

对于给定的问题，我不希望这比拒绝采样方法好得多（如果有的话），但是如果您想从截断的正常 0,1 生成 2 到 3 之间的数据，那么这种方法可能会很多更高效。它确实取决于累积及其倒数（在这种情况下为 pnorm 和 qnorm），因此不会像拒绝抽样那样简单，因为没有那些容易获得的分布。

于 2012-12-26T16:55:01.443 回答