Ugo - 1 year ago 72

R 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

But I get:

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

Not sure what is wrong.

Answer Source

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