Ugo - 1 year ago 87
R Question

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.

`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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download