2016-11-26 1 views
1
## 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 

をpi`を推定するしかし、私は得る:`norm`機能付きエラー単位円上のモンテカルロ・シミュレーションを使用して`

規範で「エラー(X [I、]): 'A'数値行列」

間違っているのかわからないでなければなりません。それが行列を要求するため

答えて

2

normは、あなたにエラーが発生します。しかし、x[i, ]は行列ではなくベクトルです。つまり、行列から単一の行/列を抽出すると、その次元は削除されます。 x[i, , drop = FALSE]を使用して、行列クラスを維持することができます。

第二の問題は、あなたがここにL2ノルムをしたい、です。したがって、標準の内部にtype = "2"を設定してください。全体を使用する

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が「円内」を示すと、論理ベクトルを得ます。さて、あなたが望むn値はTRUEの数だけあります。この論理ベクトルを合計してNを除算するか、単にこのベクトルを平均します。

関連する問題