Search code examples
rsimulationmontecarlopi

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

## 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

But I get:

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

Not sure what is wrong.

Solution

  • norm gives you error, because it asks for a matrix. However, x[i, ] is not a matrix, but a vector. In other words, when you extract a single row / column from a matrix, its dimension is dropped. You can use x[i, , drop = FALSE] to maintain matrix class.

    The second issue is, you want L2-norm here. So set type = "2" inside norm. Altogether, use

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

    norm is not the only solution. You can also use either of the following:

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

    and in fact, they are more efficient. They also underpin the idea of using rowSums in the vectorized approach below.

    Vectorization

    We can avoid the loop via:

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

    sqrt(rowSums(x ^ 2)) gives L2-norm for all rows. After comparison with 1 (the radius) we get a logical vector, with TRUE indicating "inside the circle". Now, the value n you want is just the number of TRUE. You can sum over this logical vector then divide N, or simply take mean over this vector.