• Tweet
• 

• Tweet
• 

今回の話題はこれです。

ジョーカーを除いたトランプ５２枚の中から１枚のカードを抜き出し、
表を見ないで箱の中にしまった。
そして、残りのカードをよく切ってから３枚抜き出したところ、
３枚ともダイヤであった。
このとき、箱の中のカードがダイヤである確率はいくらか。
答えが１／４ってのは納得出来ない！
１０／４９だろ！！

定期的にあちこちで話題になる確率の問題です。

調べると2chまとめサイトなどでも話題にされていますし、考察しているブログなども多数見つかります。
Yahoo!知恵袋では早稲田大学の過去問だという方がいました。

論調もいろいろで、1/4だよ派と10/49だよ派だけに限らず、本文の解釈の仕方で違うだの、この問題は単純な確率の問題ではないだのと様々な意見が出されています。

実際には数学的にきちんと議論できて正解も出せるのですが、それでは納得できない人は大勢います。数式で考えるというのは抽象化のレベルが相当高いので、計算上はともかく問題文を本当に反映できているのか?という疑問が起きるのは不思議なことではありません。

あれこれ議論する前に、実際にやってみればいいんですね。

今回はその結果の話と、結果の「正しさ」の扱いについて少し書いておきます。

手作業ではなくコンピュータでやります。コンピュータのほうがシャッフルもむしろ信頼できますし、何より回数が稼げます。早い話がモンテカルロ法です。

数式はあえて最小限にしてあります。

---

実装

「AをしたところBであった」という事象を用意する最もシンプルな実装は「Aをしてみて、Bであった場合のみ考え、あとは捨てる」という方法でしょう。今回もこれでいきます。

実際には3枚引いてダイヤになることなど滅多にないので、ほとんどの試行が無視されることになります。

C++で書いたソースコードは以下のようになりました。

#include <iostream>
#include <ctime>
#include <cmath>
#include <random>
#include <functional>
#include <vector>
#include <algorithm>
#include <iomanip>

bool is_dia(const int card){
  if(card < 13) return true;
  return false;
}

int main(int argc, char** argv){
  std::vector deck(52);
  int count = 0;
  for(auto& it : deck){
    it = count;
    count++;
  }
  time_t timer;
  time(&timer);
  auto rand = std::bind(std::uniform_int_distribution(0,51), std::mt19937(static_cast(timer)));

  long long int diaCount = 0;
  long long int notDiaCount = 0;
  long long int totalCount = 0;

  while(true){
    std::random_shuffle(deck.begin(), deck.end(), rand);
    int first_card = deck[0];
    if(is_dia(deck[1]) && is_dia(deck[2]) && is_dia(deck[3])){
      if(is_dia(first_card)) diaCount++;
      else notDiaCount++;
    }
    totalCount++;
    if(totalCount >= 1000000000) break;
  }

  long long int together = diaCount+notDiaCount;
  double ratio = static_cast(diaCount) / static_cast(diaCount+notDiaCount);
  double variance = ratio*(1.0-ratio)/together;
  double stddev = std::sqrt(variance);
  std::cout << "total: " << totalCount << std::endl
    << "dia: " << diaCount << ", other: " << notDiaCount << ", together: " << together << std::endl
    << "ratio: " << std::setprecision(10) << ratio << std::endl
    << "variance: " << variance << ", stddev: " << stddev << std::endl;
  return 0;
}

中央下あたりのwhile文が大量に試行している部分です。例によってC++11の機能を使っていますので、コンパイルには -std=c++11 などのオプションが必要になります。

52枚のカードのセットは0~51の数字の書かれた配列として表現して、これをシャッフルして使います。0~12の数字(のカード)をダイヤとしています。

捨てた回も合わせて全部で10億回ほど試行してみました。

---

結果

結果を見てみましょう。

total: 1000000000
dia: 2641328, other: 10299836, together: 12941164
ratio: 0.2041028149
variance: 1.255256914e-08, stddev: 0.0001120382485

トータル1,000,000,000回のうち有効な試行(3枚ともダイヤであった回)は12,941,164回、その中でさらに箱の中の1枚もダイヤであった回数は2,641,328回でした。

ここから、求める確率は0.20410(±0.00011)と出ました。
標準偏差(±の部分)については後述します。

この結果から1/4(=0.25)か10/49(=0.20408)のうちどちらかを選ぶなら、10/49が正解であると強く主張できることになります。

---

モンテカルロ法あれこれ

モンテカルロ法を行うときはいくつかの重要な論点があります。

1つは、本当に全体の事象を一様に再現できているかの確認が必要です。

まずは乱数に気を使う必要があります。特に古い言語(FortranとかCとか)ではデフォルトの乱数が線形合同法というクラシックな方法を使っていることが多く、モンテカルロ法には適さないことがよくあります。
今回ではカードをシャッフルするときに乱数が使われています。今回の実装ではC++11から標準で実装された乱数生成器mt19937(メルセンヌ・ツイスタ)を指定してあります。

多少余談ですが、もっと複雑な問題、例えば物理現象などを扱うときは、この状態を一様に用意するというのがかなり難しい話になってきます。

例えば「温度が300Kの理想気体の状態を100万セットほど用意したい」という要求があったとき(こういう問題はよくあります)、取り方に恣意性がないように気をつけなければいけません。

条件に合うものだけ取得するといってもあり得る状態は無限にありますし、300Kの状態の中には、確率こそ少ないもののほとんどの原子が止まっている瞬間や逆に全て激しく動いている瞬間、あるいは原子が隅に寄っているパターンもありえます。
ではどう揃えると真の「まんべんなく」になるか、というのは結構微妙な問題です。実際にこの辺りはしっかりと議論されています。

結果の正しさ

出てきた結果について「20回やったら該当5回だったからちょうど1/4だったよ」と主張しても、それは回数が少なくて偶然だったのではないかという主張が通ってしまいます。

こういう時に統計の出番です。何回やればどれだけ正しいかを表現する道具が用意されています。

今回は結果がダイヤであるかないかを集計するので、試行の結果は二項分布に従います。
確率pで起きる現象を全部でN回試行をしたとき、トータルで観測される回数は、

N回の平均:Np
N回の分散:Np(1-p)

で表現されます。ここで分散が試行回数Nに比例するというのがミソで、これの平方根をとった標準偏差σは当然Nの平方根に比例します。

N回の標準偏差:

ということは、1回あたりで見ると、両方Nで割ればいいので

確率:p
確率の標準偏差:

となります。回数を増やすほど値のブレが少なくなるんですね。(この辺りを正しく検討する場合はベータ分布の議論が必要になります。割愛しますが、結果はほぼ同じです。)

正規分布(に非常に近い二項分布も)においては、標準偏差はそのまま確からしさにつながります。工学の分野では3σという言葉があり、平均から±標準偏差の3倍の範囲にある確率は99.7%(この範囲であることはほぼ確実)であることがわかっています。

今回は結果に期待値だけでなく標準偏差も載せておきました。

私の研究室の先生は、工学上の問題においてこの統計的な考察がとても重要だとよく話しています。統計というのは正確な議論をしようとすると一筋縄ではいかない分野ですが、統計がどんな問題を扱えて何を保証してくれるのかは知っておくと便利だと思います。

また上記の議論から、モンテカルロ法は精度を10倍にするためには100倍の計算が必要になることがわかります。精度を詰めるほど計算量は勢い良く増えていくわけですね。

---

答えの補足

モンテカルロ法でやりっぱなしもあれなので、問題について一応の補足を。

巷では条件付き確率を加味した場合だけ答えが変わるだの、たまたま1回の場合と何度もやった場合で違うだのと様々に言われていますが、確率というのはそこまで主観的なものではありません。

今回の問題では条件として起こりうる全てのパターンが列挙できますので、それら全パターンをもれなくすべて書き出し、その中で条件に合うものを数え上げ、総数で割ったものがここでの「確率」の定義となります。ここでの確率の意味はそれ以上でもそれ以下でもありません。

1/4ではないかとの議論が巻き起こる理由はおそらく「はじめに選んだものの確率が、後からの行為で変わるわけがない」という信念に基づくものだと思います。
しかし実際には変わります。「他から3枚選んだ」までなら変わりませんが、そのカードの中身を知ること、つまり残りのカードの山に対する新しい知識を得たことは確率を変えるのです。

3枚でなく「残り全部のカードを見た」なら結果はさらに明快になります。(このとき箱の中の1枚はすぐ判明します。)見たのがたった3枚でも結果には少し影響しているのです。

ここから答えを10/49とぽんと出すには、さらに「箱に入れる1枚を先に引いても3枚引いた後で引いても答えは同じ」というアイデアに至ればわかります。といってもこのアイデアは、(正しいけど)証明や細かい議論なしに即座に正しいと認められるものでもないと思います。
ソースコードを読むとこの理屈が比較的すんなりと受け入れられると思います。

これ、出題もよくできているなあと思います。「正解が1/4になっているのは納得できない!」という書き方になっているところが憎いですね。赤本か何かで1/4と書いてあったのでしょうか。

---

そんなわけで、正しく使うと強力で頼もしいモンテカルロ法のお話でした。