0x19f (Shinya Kato) の日報

主にプログラミング関連の話をします

山本一成さんの「エレファントな解法」に対抗して動的計画法による「エレガントな解法」を考えてみた

エレファントな解法(?)

先日Ponanza開発者の山本一成さんがブログでこのような記事を書いておられました.

note.mu

この記事の中では以下の問いが投げかけられています.

問. コインを100回投げて表か裏が10回連続で出る確率は?

山本さんのブログでは「エレファントな解法」としてモンテカルロ法を用いた力任せの解法が紹介されています. つまりは, 乱数を使って100回コインを投げるというシミュレーションを何度も繰り返して条件を満たす確率を近似的に求めるというわけです. この方法ならきちんと確率を求めることができなくても, それなりに確率を正しく求められます. 困った時の最終手段としては非常に便利だと思います.

ちなみに, こういった方法が「エレファントな解法」と呼ばれているのは初めて知りました.

エレガントな解法

「エレファントな解法」のことは一旦置いておいて, 山本さんがエレガントな解法について「私は正直解ける気がしません」と言っておられたのが少し気になりました. 競プロerならぱっと見で「動的計画法で解けそう?」と感じた人も多いんじゃないでしょうか, 僕もそう感じました.

というわけで, この問題を動的計画法(DP)で解いてみたいと思います.

DPで解いてみる

先に断っておくと, 解いてみてから気づいたのですがちょっと冗長なDPをしています. あと, 答えはあったものの, 若干怪しいところがあります(おい).

基本的な考え方

DPそのものの説明はいろいろと調べてもらうのが早いと思うのですが, 僕は「コインの投げ方って2N通りあるけれど, 実はそんなに場合分けしなくても"何回コインを投げて何回連続で同じ面が出てるか"で場合分けすれば十分じゃない?」というものだと思ってます.

例えば, 10回コインを投げて最後に表が3回連続で出るような投げ方って26通りありますが, 今回の問題では最初の6回までに表裏どっちが出たかって知らなくても答えは出せます.

というわけで, これらをまとめて計算して, さらに上手に前の計算結果を利用していくことで計算量を落とすことができます.

解法

まず, コインを投げる回数をN, 求めたい確率の連続回数をMとします. そして, 以下のようなDPを考えます.

dp[i][j][k][l]:
i回コインを投げた時に, 最後にj回連続で(k = 0: 表, 1: 裏)が出て, 条件を(l = 0: 満たしていない, 1: 満たしている)確率

遷移に関しては, (k, l)の4パターンについてそれぞれ表が出た場合と裏が出た場合を考えます. 表が連続している場合に限って説明します.

まだ条件を満たしていいない場合は, 表が出たら連続回数を1増やし裏が出た場合は連続回数を1に戻してあげます. このとき, 表がM回以上連続していたらl = 1, つまり条件をすでに満たした状態に遷移してあげます.

すでに条件を満たしている場合も, 表が出たら連続回数を1増やし裏が出た場合は連続回数を1に戻してあげます.

裏が連続している場合もほぼ同様です.

最終的な答えは, N回コインを投げて最後の連続数とコインの裏表に関係なく条件を満たしている確率です. つまり, dp[N][j][k][1] (0 <= j < N, k = 0, 1)の和が答えです.

というわけで, O(N2)の時間計算量で答えが出せます.

実装

これをC++で書いてみると以下のようになります.

#include <cstdio>

int main(void) {
  const int N = 100, M = 10; // コインを投げる回数, 連続する回数
  const double P = 0.5;      // 表が出る確率

  // dp[i][j][k][l]:
  //   i回コインを投げて,
  //   最後のj回連続して(k = 0: 表, 1: 裏)が出て,
  //   (l = 0 : まだ条件を満たしていない, 1: すでに条件を満たしている)確率
  double dp[N + 1][N + 1][2][2];
  for(int i = 0; i < N + 1; i++) {
    for(int j = 0; j < N + 1; j++) {
      dp[i][j][0][0] = 0;
      dp[i][j][0][1] = 0;
      dp[i][j][1][0] = 0;
      dp[i][j][1][1] = 0;
    }
  }
  dp[0][0][0][0] = P;
  dp[0][0][1][0] = (1 - P);

  for(int i = 0; i < N; i++) {
    for(int j = 0; j < i + 1; j++) {

      // (k, l) = (0, 0): 表が連続, 条件を満たしてない
      dp[i + 1][j + 1][0][j + 1 >= M] += dp[i][j][0][0] * P;
      dp[i + 1][1][1][0] += dp[i][j][0][0] * (1 - P);

      // (k, l) = (0, 1): 表が連続, 条件を満たした
      dp[i + 1][j + 1][0][1] += dp[i][j][0][1] * P;
      dp[i + 1][1][1][1] += dp[i][j][0][1] * (1 - P);

      // (k, l) = (1, 0): 裏が連続, 条件を満たしてない
      dp[i + 1][1][0][0] += dp[i][j][1][0] * P;
      dp[i + 1][j + 1][1][j + 1 >= M] += dp[i][j][1][0] * (1 - P);

      // (k, l) = (1, 1): 裏が連続, 条件を満たした
      dp[i + 1][1][0][1] += dp[i][j][1][1] * P;
      dp[i + 1][j + 1][1][1] += dp[i][j][1][1] * (1 - P);
    }
  }

  // 最終的な答えは, dp[N][j][k][1] (0 <= j < N, k = 0, 1)の和
  double ans = 0;
  for(int j = 0; j < N + 1; j++) {
    ans += dp[N][j][0][1];
    ans += dp[N][j][1][1];
  }
  printf("%.15lf\n", ans);
}

実行した結果「コインを100回投げて表か裏が10回連続で出る確率」は0.086659044348362となりました. 山本さんがブログの中で知人の数学者の方に「エレガントな解法」で計算してもらったという値と一致してます, やったね.

まとめ

というわけで, こんな感じでこの問題は答えが出せます. これぐらいならまだなんとか確率求められますが, 難しくなると圧倒的にエレファントな解法に分がありそうですね(無理やりまとめた).

ちなみに, 頑張るとO(NM)で答え出せるはずです. この解法はchokudaiさんが書いておられますね......プログラムみじかw