PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 389 / 正多面体のサイコロ

1つの偏りのない四面体のサイコロを振り,出た目 T を記録する。
T 個の偏りのない六面体のサイコロを振り, 出た目の合計 C を記録する。
C 個の偏りのない八面体のサイコロを振り, 出た目の合計 O を記録する。
O 個の偏りのない十二面体のサイコロを振り, 出た目の合計 D を記録する。
D 個の偏りのない二十面体のサイコロを振り, 出た目の合計 I を記録する。
I の分散を四捨五入して小数点以下第4位まで求めよ。

Problem 389 - Project Euler

まずCの分散を求める

C の分散をどう求めるかを考えて,それを一般化しました。

C の確率分布

C の確率分布は一応求められます。

 P(C=c)=\displaystyle\sum_t P(T=t)P(C=c | T=t)

この式の P(C=c | T=t) は条件つき確率です。T=t の条件下で C=c となる確率を表します。P(後|前) のように書きます。

同じ方法で I の確率分布も表せますが,こういうダサい方法では解きたくないので Law of Total Variance(Eve's law)を使います。

Law of Total Variance(Eve's law)

 Var(C)=E(Var(C|T))+Var(E(C|T))

この公式は wikipedia にも載っています*1が,こちらのページがわかりやすいです。

www.probabilityformula.org

Var(C) の各項の計算に入りましょう。E(C|T) は T がわかっているときの C の期待値です。六面体のサイコロを T 個振ったときの目の和が C なので,1個だけ振ったときの目の T 倍になります。六面体のサイコロを1回振ったときの目を c とすると

 E(C|T)=E(c)T

E(c) は定数なので V(aX)=a^2 V(X) が使えて,Var(C) の右辺第2項は簡単な形になります。

 Var(E(C|T)) =Var(T) \left\{E(c)\right\}^2

次に Var(C) の右辺第1項を求めます。Var(C|T) はTがわかっているときの C の分散です。Cは六面体のサイコロを T 個振ったときの目の和で,各サイコロの出目は独立なので V(X1+X2)=V(X1)+V(X2) が使えます。

 Var(C|T)=Var(c)T
 \therefore E(Var(C|T))= Var(c) E(T)

以上を Var(C) の式に代入します。

 Var(C)=Var(c) E(T)+Var(T) \left\{E(c)\right\}^2

n 面体のサイコロでは E=(n+1)/2, V=(n^2-1)/12 となることを使うと,Var(C) が計算できます。

 Var(C)=\dfrac{6^2-1}{12}\cdot \dfrac{5}{2}+\dfrac{4^2-1}{12}\cdot \left(\dfrac{7}{2}\right)^2=\dfrac{1085}{48}

期待値も必要

上で導いた式を繰り返し使うには E(T), E(C), ……を求めておく必要があります。Adam's law というらしいのですが,この問題に使うのはちょっと大げさです。

 E(Y)=E(E(Y|X))

(サイコロを1個振ったときの目の平均値)*(サイコロを振る回数の期待値)で求めれば十分でしょう。

E(T)=(1+4)/2=5/2
E(C)=E(c)E(T)=(1+6)/2*5/2=35/4
E(O)=E(o)E(C)=315/8
E(D)=E(d)E(O)=4095/16
E(I)=E(i)E(D)=85995/32

解答

係数を全部決めてしまってハードコーディングすることも考えましたが,一応,漸化式にまとめました。
答えの厳密な値は 2464129395/1024 です。


統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)