在单位圆上使用 Monte Carlo 模拟估算 pi 时，norm 函数出错

Question

## simulate `N` uniformly distributed points on unit square
N <- 1000
x <- matrix(runif(2 * N), ncol = 2)

## count number of points inside unit circle
n <- 0; for(i in 1:N) {if (norm(x[i,]) < 1) {n <- n + 1} } 
n <- n / N 

## estimate of pi
4 * n

但我得到：

"Error in norm(x[i,]): 'A' must be a numeric matrix"

不确定哪里出了问题。

Answer 1

norm 给你错误，因为它要求一个矩阵。但是，x[i, ]不是矩阵，而是向量。换句话说，当您从矩阵中提取单个行/列时，它的维度将被删除。您可以使用 x[i, , drop = FALSE] 来维护矩阵 class.

第二个问题是，你要在这里L2-norm。所以设置 type = "2" inside norm。总而言之，使用

norm(x[i, , drop = FALSE], type = "2") < 1

norm 不是唯一的解决方案。您还可以使用以下任一方法：

sqrt(c(crossprod(x[i,])))
sqrt(sum(x[i,] ^ 2))

事实上，它们更有效率。它们还支持在下面的矢量化方法中使用 rowSums 的想法。

矢量化

我们可以通过以下方式避免循环：

n <- mean(sqrt(rowSums(x ^ 2)) < 1)  ## or simply `mean(rowSums(x ^ 2) < 1)`

sqrt(rowSums(x ^ 2)) 为所有行提供 L2 范数。与 1（半径）比较后，我们得到一个逻辑向量，TRUE 表示 "inside the circle"。现在，您想要的值 n 就是 TRUE 的数量。您可以对这个逻辑向量求和然后除以 N，或者简单地对这个向量取平均值。

在单位圆上使用 Monte Carlo 模拟估算 pi 时，norm 函数出错

Error with `norm` function when estimating `pi` using Monte Carlo simulation on a unit circle

r

simulation

montecarlo

pi