2

私は統計学のクラスを教えており、学生に R を使用したシミュレーションを通じて確率と統計学の問題を探させています。最近、5 つのサイコロを振ったときにちょうど 2 つの 6 が出る確率について混乱がありました。答えは choose(5,2)*5^3/6^5 ですが、一部の学生は「順序は問題にならない」と確信していました。つまり、答えは choose(5,2)*choose(25,3)/choose(30,5) になります。5つのサイコロを何千回も振ってシミュレートし、各実験の経験的確率を追跡し、実験を何度も繰り返すと面白いと思いました。問題は、上記の 2 つの数値が十分に近いため、統計的に有意な方法で違いを引き出すためのシミュレーションを取得するのが非常に難しいことです (もちろん、私のやり方が間違っている可能性もあります)。5個のサイコロを100000回振ってから、実験を10000回繰り返してみました。これを i7 Linux マシンで実行するのに 1 時間ほどかかりましたが、25% の確率で正しい答えが choose(5,2)*choose(25,3)/choose(30,5) になる可能性がありました。そこで、1 回の実験でサイコロを振る回数を 10^6 に増やしました。現在、コードは 2 日以上実行されており、終了の兆候は見られません。操作の数を 1 桁だけ増やしただけなので、実行時間が 10 時間近くになるはずであることを意味しているため、これには混乱しています。現在、コードは 2 日以上実行されており、終了の兆候は見られません。操作の数を 1 桁だけ増やしただけなので、実行時間が 10 時間近くになるはずであることを意味しているため、これには混乱しています。現在、コードは 2 日以上実行されており、終了の兆候は見られません。操作の数を 1 桁だけ増やしただけなので、実行時間が 10 時間近くになるはずであることを意味しているため、これには混乱しています。

2 番目の質問: これを行うためのより良い方法はありますか? 以下に投稿されたコードを参照してください。

probdist = rep(0,10000)

for (j in 1:length(probdist))
{
   outcome = rep(0,1000000)
   for (k in 1:1000000)
   {
      rolls = sample(1:6, 5, replace=T)
      if (length(rolls[rolls == 6]) == 2) outcome[k] = 1
   }

   probdist[j] = sum(outcome)/length(outcome)
}
4

3 に答える 3

3

経験則としては、ループを絶対に記述しないことforRです。別の解決策は次のとおりです。

doSample <- function()
{
   sum(sample(1:6,size=5,replace=TRUE)==6)==2
}

> system.time(samples <- replicate(n=10000,expr=doSample()))
user  system elapsed 
0.06    0.00    0.06 
> mean(samples)
[1] 0.1588
> choose(5,2)*5^3/6^5
[1] 0.160751

$10,000$ のサンプルではあまり正確ではないようです。100,000ドルの方が良い：

> system.time(samples <- replicate(n=100000,expr=doSample()))
user  system elapsed 
0.61    0.02    0.61 
> mean(samples)
[1] 0.16135
于 2013-10-30T16:24:27.670 に答える
2

ほとんどの場合、for ループよりもベクトル化が優先されます。この場合、最初にすべてのサイコロを投げてから、5 個の各グループのいくつが 6 に等しいかを確認することで、大幅なスピードアップが見られるはずです。

set.seed(5)
N <- 1e6
foo <- matrix(sample(1:6, 5*N, replace=TRUE), ncol=5)
p <- mean(rowSums(foo==6)==2)
se <- sqrt(p*(1-p)/N)
p
## [1] 0.160382

95% 信頼区間は次のとおりです。

p + se*qnorm(0.975)*c(-1,1)
## [1] 0.1596628 0.1611012

真の答えans1（ans2真の答えをテストするときの p 値は 0.31 ですが、偽の答えの場合は 0.0057 です。

(ans1 <- choose(5,2)*5^3/6^5)
## [1] 0.160751
pnorm(abs((ans1-p)/se), lower=FALSE)*2
## [1] 0.3145898

ans2 <- choose(5,2)*choose(25,3)/choose(30,5)
## [1] 0.1613967
pnorm(abs((ans2-p)/se), lower=FALSE)*2
## [1] 0.005689008

一度にすべてのサイコロを投げていることに注意してください。メモリが問題になる場合は、元の投稿で行ったように、これを分割して組み合わせることができます。これがおそらく、予想外のスピードアップを引き起こした原因です。スワップ メモリを使用する必要がある場合は、大幅に遅くなります。その場合は、ループ内のロールの数ではなく、ループを実行する回数を増やすことをお勧めします。

于 2013-12-09T18:46:58.130 に答える
2

私はもともと、M. Berk が R の duplicate() 関数を使用するよう提案したことに対して、正解チェックを授与していました。さらなる調査により、以前の支持を取り消すことを余儀なくされました。replica() は sapply() の単なるラッパーであり、実際には for ループよりもパフォーマンス上の利点はありません (これはよくある誤解のようです)。いずれにしても、実行時間を比較するために、3 つのバージョンのシミュレーションを用意しました。 :

# dice26dist1.r: For () loop version with unnecessary array allocation
probdist = rep(0,100)

for (j in 1:length(probdist))
{
  outcome = rep(0,1000000)
  for (k in 1:1000000)
  {
    rolls = sample(1:6, 5, replace=T)
    if (length(rolls[rolls == 6]) == 2) outcome[k] = 1
  }
  probdist[j] = sum(outcome)/length(outcome)
}

system.time(source('dice26dist1.r'))
ユーザー システム経過時間
596.365 0.240 598.614

# dice26dist2.r: For () loop version
probdist = rep(0,100)

for (j in 1:length(probdist))
{
  outcomes = 0
  for (k in 1:1000000)
  {
    rolls = sample(1:6, 5, replace=T)
    if (length(rolls[rolls == 6]) == 2) outcomes = outcomes + 1
  }
  probdist[j] = outcomes/1000000
}

system.time(source('dice26dist2.r'))
ユーザー システム経過時間
506.331 0.076 508.104

# dice26dist3.r:  replicate() version
doSample <- function()
{
   sum(sample(1:6,size=5,replace=TRUE)==6)==2
}

probdist = rep(0,100)

for (j in 1:length(probdist))
{
  samples = replicate(n=1000000,expr=doSample())
  probdist[j] = mean(samples)
}

system.time(source('dice26dist3.r'))
ユーザー システム経過時間
804.042 0.472 807.250

このことから、replicate() バージョンは、いずれの for ループ バージョンよりも、どの system.time メトリックでもかなり遅いことがわかります。私の問題は、100 万文字の output[] 配列を割り当てることによるキャッシュ ミスが主な原因であると当初は考えていましたが、dice26dist1.r と dice26dist2.r の時間を比較すると、これはパフォーマンスにわずかな影響しか与えていないことがわかります (ただし、システムへの影響はかなりの時間がかかります: >300% の差。

3 つのシミュレーションすべてでまだ for ループを使用していると主張する人もいるかもしれませんが、私が知る限り、ランダム プロセスをシミュレートする場合、これは完全に避けられません。毎回、実際にランダムなプロセス (この場合は 5 つのサイコロを振る) をシミュレートする必要があります。for ループの使用を回避できるようにする手法について知りたいです (もちろん、パフォーマンスを向上させる方法で)。この問題が並列化に非常に効果的に役立つことは理解していますが、単一の R セッションの使用について話しているのですが、これを高速化する方法はありますか?

于 2013-10-30T16:43:28.747 に答える