前書き
私の記事では幾度となく出てくる、ザ・ゴール漫画版の話である。
これは、ケントベックが書いたエクストリームプログラミングの元ネタの一冊である。
エクストリームプログラミングには制約理論(TOC)という章があるが、その制約理論に関して説明しているのが「ザ・ゴール」である。
このザ・ゴールに「サイコロゲーム」というちょっと変わったゲームが出てくる。そのゲームに関する数理解析である。
依存的事象と統計的変動
主人公は、上層部から閉鎖の圧力を受けている工場の責任者である。この主人公が、アドバイザーであるジョナとの会話の中で、「バランスのとれた工場」について話す描写がある。「バランスのとれた工場」とは
すべてのリソースの生産能力が市場の需要と完全にバランスが取れている工場だ
とされている。市場の需要より生産能力が劣っていた場合、機会損失をしていることになるし、市場の需要より生産能力が多い場合は、過剰投資をしていることになる。
そう考えた場合、生産能力を需要により増減させることができればよいのではないか。と考える。しかし、それはできない。とジョナが言う。その制御できない理由が、
- 依存的事象
- 統計的変動
- 依存的事象・統計的事象の2つが同時に起こったときの効果
と説明されている。
これらを説明すると、例えば、「近くのラーメン屋でラーメンを食べる」ということを考える。その時、
- 家からラーメン屋に行く
- ラーメンを注文する
- ラーメンを食べる
- ラーメン屋から家に帰る
というステップを踏むだろう。この、1,2,3,4という順番通りに行わなければならない手順（ステップ・工程）を「依存的事象」と呼んでいる。
「家からラーメン屋に行く」ことを考えたとき、その間には信号があるだろう。そのとき、信号が赤であれば止まる必要があるし、運よく青であれば、そのまま通過することができるだろう。このようにラーメン屋に行くだけでも、間のランダムな要素によって、その所要時間に幅がある。また、「ラーメンを注文する」でも、お店の混み具合によっては中々オーダーが通らなかったり、提供に時間がかかる。ここにも所要時間に幅がある。これを「統計的変動」と呼んでいる。
「依存的事象」と「統計的変動」の2つの事象が組み合わさったときに、人間の単純な思考が外れるようになる。 というのが、ザ・ゴールのポイントの1つとなる。
サイコロゲーム
第５章 窮余の一策では、サイコロゲームというものが出てくる。前章の、「依存的事象」「統計的変動」の2つを理解するために主人公が作ったゲームだ。
これはマッチ棒が何本移動できるか？を競うゲームである。
準備するものは4つである。
- 5人のプレイヤー
- 5つのカップ
- 5つのサイコロ
- 60本のマッチ
ゲームのルール
①プレイヤー5人に対し、1つずつカップを配布する。
②カップ1に60本のマッチ棒を入れる。
③プレイヤー1がサイコロを振り、出た目の数だけプレイヤー2にマッチ棒を渡す。
④プレイヤー2がサイコロを振り、出た目の数だけプレイヤー3にマッチ棒を渡す。ただし、自分のカップ内にあるマッチ棒の数以上は渡せない。例えば、カップ内にマッチ棒が3本で、サイコロの出た目が6のときは、隣の人に3本だけマッチ棒を渡す。
⑤④の工程をプレイヤー3～5まで行う。プレイヤー5が出したマッチ棒の数をスコアとする。
⑥③～⑤の工程を10周行う。
⑦プレイヤー5が出したマッチ棒の総数が最終スコアとなる。
主人公の思考と実際の乖離
主人公は以下のように考えた。
「サイコロの出る目の期待値は3.5本である。ということは1週でプレイヤー5が出すスコアも3.5本であろう。それを10周すること言うことは、最終スコアは35本に近くなるだろう。」
このゲームを行った結果、書籍上では20本しかスコアが出なかった。これはどのようなことが起こっていたのだろうか。
サイコロゲームの数理解析
この記事で示しているルールだが、書籍内の内容とは少し変更がある。それは、「カップ1に入れる本数を60本」としている点だ。これは、解析を行う際に都合がよかったために行った(本編では、とりあえず手持ちのマッチ棒をたくさん入れておく。ぐらいの描写である)。サイコロの目の最大値は6である。これがゲームにおいて出続けた場合どうなるだろうか。すべての工程において6が出るため、1週あたりのプレイヤー5が出すマッチ棒は6本である。それを10周するため、最終に移動するマッチ棒は60本となる。これによりゲームにおける利用する最大のマッチ棒の数が分かった。
今回の数理解析であるが、確率分布を用いる。
p_{i,j}(n=m): プレイヤーiがj週目においてカップ内にm本のマッチを持っている確率
と定義する。
プレイヤー1は初期状態(1週目)として、60本のマッチを持った状態である。その時の状態が100%であるため、以下のような定義になる。
\begin{eqnarray}
p_{1,1}(n)
=
\begin{cases}
1.0 & ( n = 60 ) \\
0.0 & ( 0 \le n \lt 60 )
\end{cases}
\end{eqnarray}
また、プレイヤーが持っているマッチの数は、どんな時も必ず0本から60本のいずれかである。そのため、それらの総和は必ず1.0となる。
\sum_{n=0}^{60} p_{i,j}(n) = 1.0
次に、プレイヤーが渡すマッチの本数の確率分布について考える。
q_{i,j}(n=m): プレイヤーiがj週目において、m本のマッチを渡す確率
と定義する。この時、プレイヤーは必ず1～6本のマッチを渡すので、それらの総和は1.0となる。
\sum_{n=1}^6 q_{i,j}(n) = 1.0
ここでプレイヤー1が1週目に渡すマッチの本数について考える。プレイヤー1は初期状態として、60本持っているため、サイコロの出目が1～6どれが出ても、プレイヤー2にマッチを渡すことができる。そのため、すべての値が1/6となる。
q_{1,1}(n) = {1 \over 6}
ここから、プレイヤー1の2週目のカップ内のマッチの本数を考えると以下のようになる。
\begin{eqnarray}
p_{1,2}(n)
=
\begin{cases}
{1 \over 6} & ( n \in \lbrace 54,55,56,57,58,59 \rbrace ) \\
0.0 & ( それ以外)
\end{cases}
\end{eqnarray}
プレイヤー1の視点だと、もともとマッチが60本あり、そこから1～6本の間で、等確率でマッチの本数が減るので、上式のようになる。
またプレイヤー2の1週目のマッチの本数は以下のようになる。
\begin{eqnarray}
p_{2,1}(n)
=
\begin{cases}
{1 \over 6} & ( n \in \lbrace 1,2,3,4,5,6 \rbrace ) \\
0 & ( それ以外)
\end{cases}
\end{eqnarray}
プレイヤー2の視点だと、マッチを1～6本を等確率でもらうので、上式のようになる。
ここで、特殊なケースについて考える。カップのマッチの本数<サイコロの目の数である場合だ。ここでは本来存在しないが、わかりやすい例を考えてみる。仮にプレイヤー3がJ週目に、3本しかマッチを持っていないとする。
\begin{eqnarray}
p_{3,J}(n)
=
\begin{cases}
1.0 & ( n = 3 )\\
0.0 & ( n \neq 3)
\end{cases}
\end{eqnarray}
この時、渡すマッチの確率が以下のように書ける。
\begin{eqnarray}
q_{3,J}(n)
=
\begin{cases}
1 \over 6 & ( n \in \lbrace 1,2 \rbrace )\\
2 \over 3 & ( n = 3 ) \\
0 & (n \in \lbrace 4,5,6 \rbrace)
\end{cases}
\end{eqnarray}
サイコロを振り、それぞれの目がでる確率は1/6である。そのため、1,2が出る確率は1/6である。しかし、マッチが3本しか入っていないとき、目が4,5,6と出ても3本しか渡せない。そのため、その部分の確率が吸われ、3本出る確率が2/3となり、4,5,6本渡す可能性が0となる。
このようにカップの中にあるマッチの本数を条件に渡すマッチの本数の確率分布がぶれるため、数式としては表現しにくい。そのため、これらをプログラムにして計算を行った。
プログラムによる数理解析
プログラムにより解析した最終的なマッチ棒の分布は以下のようになります。
sympyによる解析プログラム
from itertools import combinations
import sympy as sp
import matplotlib
import numpy as np
from typing import List
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from filecache import filecache
one = sp.Rational(1, 1)
zero = sp.Rational(0, 1)
@filecache
def simulate_dice_game(players: int, round_limit: int, dice: List[int]) -> List[sp.Rational]:
dice_kind = len(dice)
dice_distribution = [sp.Rational(1, dice_kind)] * dice_kind
dice_map = dict(zip(dice, dice_distribution))
max_matches = max(dice) * round_limit
initial_bowl_matches_distribution = [zero] * (max_matches) + [one]
bowl_distributions = [initial_bowl_matches_distribution] + [[one] + [zero] * max_matches for _ in range(players)]
for n in range(round_limit):
for i in range(players):
bowl_distributions[i],next_handed_match_distribution = update_self_bowl_distribution(dice_map, bowl_distributions[i])
bowl_distributions[i + 1] = update_next_bowl_distribution(bowl_distributions[i+1], next_handed_match_distribution)
return bowl_distributions[-1]
def update_self_bowl_distribution(dice_map, target_distribution):
update_self_bowl_distribution = [zero] * len(target_distribution)
max_dice = max(dice_map.keys())
next_handed_match_distribution = dict(zip(range(1,max_dice+1), [zero] * max_dice))
#print(next_handed_match_distribution)
#print(target_distribution)
for j, prob in enumerate(target_distribution):
if prob == 0:
continue
for dice_value, dice_prob in dice_map.items():
if j >= dice_value:
update_self_bowl_distribution[j - dice_value] += prob * dice_prob
next_handed_match_distribution[dice_value] += prob * dice_prob
else:
update_self_bowl_distribution[0] += prob * dice_prob
next_handed_match_distribution[j] += prob * dice_prob
return update_self_bowl_distribution,next_handed_match_distribution
def update_next_bowl_distribution(target_distribution, next_handed_match_distribution):
bowl_distribution = [zero] * len(target_distribution)
for j, prob in enumerate(target_distribution):
if prob == 0:
continue
for dice_value, dice_prob in next_handed_match_distribution.items():
bowl_distribution[j + dice_value] += prob * dice_prob
return bowl_distribution
def plot_results(distribution: List[sp.Rational], title: str, filename: str) -> None:
matches_distribution_length = len(distribution)
fig, ax1 = plt.subplots(figsize=(12, 8))
color = 'tab:blue'
ax1.set_xlabel('Number of Matches')
ax1.set_ylabel('Probability', color=color)
ax1.bar(range(matches_distribution_length), distribution, color=color, alpha=0.6, label='Probability')
ax1.tick_params(axis='y', labelcolor=color)
cumulative_distribution = [sum(distribution[:i + 1]) for i in range(matches_distribution_length)]
ax2 = ax1.twinx()
color = 'tab:red'
ax2.set_ylabel('Cumulative Probability', color=color)
ax2.plot(range(matches_distribution_length), cumulative_distribution, color=color, marker='o', label='Cumulative Probability')
ax2.tick_params(axis='y', labelcolor=color)
ax2.yaxis.set_major_locator(plt.MultipleLocator(0.1))
mean, variance, mode, median = describe_statics(distribution)
textstr = '\n'.join((
f'Mean: {mean:.2f}',
f'Variance: {variance:.2f}',
f'Mode: {mode}',
f'Median: {median}'
))
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
ax1.text(0.05, 0.95, textstr, transform=ax1.transAxes, fontsize=14,
verticalalignment='top', bbox=props)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.title(title)
plt.grid(True)
plt.savefig(filename)
def describe_statics(distribution):
matches_distribution_length = len(distribution)
cumulative_distribution = [sum(distribution[:i + 1]) for i in range(matches_distribution_length)]
mean = np.average(range(matches_distribution_length), weights=distribution)
variance = np.average((np.array(range(matches_distribution_length)) - mean) ** 2, weights=distribution)
mode = np.argmax(distribution)
median = np.searchsorted(cumulative_distribution, 0.5)
return mean,variance,mode,median
if __name__ == "__main__":
players = 5
round_limit = 10
dice = [1,2,3,4,5,6]
distribution = simulate_dice_game(players, round_limit, dice)
plot_results(
distribution,
'Distribution and Cumulative Distribution of Matches in the Last Bowl (c5)',
'distribution_and_cumulative_distribution_c5.png'
)
横軸がマッチの本数です。棒グラフが、その確率を示し、左の軸に対応しています。赤い折れ線グラフは累積分布になっており、右の軸と対応しています。
解析結果を見ると、平均値が23.18本、中央値・最頻値ともに23本と結論付けられています。主人公は20本しか動かせなかった。というのは、それほど悪い結果ではないと思います。一方で、主人公の見立ての35本というのはとても非現実的であったことも分かります。累積分布で見ても35本以上は1%もありません。
理解のためのシミュレーション
主人公の主張が正しいと仮定します。それは、「1週の期待値が3.5であるから、10周すれば35本になるはずだ」ということです。では、「1週の期待値が3.5でさえあればいいのか？」と考えてみます。
では、ここでダイスの目が1,6しか出なかった場合のシミュレーションを行ってみます。この場合でも期待値は3.5のままです。
結果としていびつなグラフになり、平均値は18.66、最頻値が15、中央値が20となります。ここから、カップの中のマッチの残り数というものがボトルネックになり、最終的なマッチの本数が伸び悩むことが分かります。
では、一方でダイスの目が3,4だけの場合はどうでしょうか？
平均が31.73、最頻値が31、中央値が32です。となり、これほど条件を良くしても、期待値が35を超えることはありません。
直感的な理解
プログラムによる解析は厳密化もしれませんが、直感とずれてしまい、いまいち腹落ち感に欠けます、そこで簡易化したモデルについて考えてみます。
1週目の最後のプレイヤー5が期待値3.5を出すことを考えます。この時、プレイヤー5は4,5,6の目を出しているとよさそうです。そのとき、プレイヤー5がそれらの目を出す確率は1/2です。しかし、そのためには、プレイヤー1～4も4,5,6の目を出し続ける必要がありそうです。その確率は(1/2)^4になります。
ここで、プレイヤー5が4,5,6の目を出したときに動かすマッチの本数を6本、そうでないときを3本と過大に見積もって、その期待値を計算してみます。
6 \times ({1 \over 2}) ^ 4 \times {1 \over 2} + 3 \times ( 1 - ({1 \over 2})^4 ) \times {1 \over 2} \simeq 1.59
プレイヤー5が4,5,6の数字を出すためには、プレイヤー1～4も4,5,6を出す必要があるため(1/2)^4の項が必要になります。プレイヤー5が1,2,3の数字を出す時、プレイヤー1～4の挙動は、先ほどのプレイヤー1～4の振る舞いの余事象となるため(1-(1/2)^4)としています。プレイヤー5が4,5,6の目を出す時のマッチの本数は6ではありませんが、ここは過大に見積もっています。同様に、プレイヤー5が1,2,3の目を出す時、この時、プレイヤー5の目の最大値は3となります。そのため、ここも過大に見積もり3としています。
そのように過大に見積もったとしても、1週当たりのマッチの本数は1.59本しか渡せません。ここから10周を見積もっても、15.9本しか期待値として得られないのが分かります。
ただし、厳密なシミュレーションの平均が23本というところと乖離はあります。以上の概算は、確かに上振れするようにモデル化していますが、カップの中にマッチが残り、それにより渡せるマッチの本数が増える。という部分がありません。そのため、厳密なシミュレーションより値が下回っています。
依存的事象と統計的変動
ここでは、依存的事象と統計的事象の振る舞いを分離して解析を試みる。プレイヤー数を変動させつつ、サイコロの目の出る期待値を3.5のままいろいろなパターンをシミュレーションし、その平均値をプロットした。
プレイヤーの数が増えれば増えるほど、期待値が下がっていくことは直感的である。サイコロを振った際、少なくとも3以上の出目が出続けなければ、期待値が3.5を超えない。そう考えたとき、一人プレイヤーが増えるたびに期待値が3.5を超える確率は2/3ずつ厳しくなっていく。そう考えるとプレイヤーの数が増えると期待値が下がっていくことに違和感は少ない。
また、出目が3,4のときというのは安定的に期待値が安定的に3.5、少なくとも期待値3以上は続けて出るような印象がある。そのため、Dice(3,4)の期待値がとても高いのは納得できる。これは、先の解析と同じだが、1,6の目しか出ないDice(1,6)はかなり揺れ幅のある動きをする。そのため、プレイヤー数が増えるごとに期待値の減少具合が大きい。
一方で驚くべき動きをしているのはDice(1,2,5,6)である。これはプレイヤー5人のときのグラフも出力した。
はっきり言ってなぜかはわからない。プログラムのバグかもしれない。結果を正しいとするならば、この時にのみ奇跡的にバランスし、期待値が35本になる。ただし、振れのある答えでもある。というのは分散が42もあり、上下に裾の長い分布になっているため、30～40本である確率がおおむね60%である。先ほど解析した、出目が3,4のときの場合、30～40本である確率が100%である。そのため、出目が1,2,5,6のパターンは出目に幅を持たせた結果、期待値的には高いかもしれないが、平均的には低い値も出る可能性があるという感じである。一方で、プレイヤー数6人以上の挙動は本当に分がわからない。ほかのデータの傾向を見るに、こちらの方が振る舞いとして正しそうな雰囲気はある。
微分の解析
出目が[1,6],[2,5],[3,4]の場合を考えてみる。プレイヤー数を1～6人まで変更し、そのスコアの平均値を求める。それを
y = Ae^{Bx}
で近似した。結果、
Dice: Dice (1, 6)
Equation: y = 39.4843 * exp(-0.1527 * x)
R^2: 0.9809
Derivative: -6.02812380382797*exp(-0.152671312061542*x)
Dice: Dice (2, 5)
Equation: y = 36.6801 * exp(-0.0758 * x)
R^2: 0.9636
Derivative: -2.7800993242666*exp(-0.0757931175309534*x)
Dice: Dice (3, 4)
Equation: y = 35.3293 * exp(-0.0215 * x)
R^2: 0.9479
Derivative: -0.759272905429162*exp(-0.021491332858096*x)
となった。決定係数は0.9を超えているためほぼ近似できていると言えるだろう。Derivativeは近似式を微分したものである。この微分をグラフ化した。
(これは私が近似式を恣意的に選んでいる面もあるが、)サイコロの目の分散が大きいほど、期待値の下がり度合いが大きい。とは言えそうである。一方で、その下がり幅はプレイヤー数が増えると0に近づきそうである。これは、直感的にも理解できて、プレイヤー数が100人ぐらいになると、乱数が厳しくて、最終的にマッチ棒が1本も渡らなさそうである。そのあと、プレイヤー数が101人に増えても、マッチ棒は1本も渡らなさそうである。そう考えた場合、変化率は0に近づくと推察できる。
また、これらのグラフは交わらなさそうである。そのため、プレイヤー数が増えても、サイコロの目の分散が大きいほうが、下がり幅が大きいと考えられる。
これは、ある意味で依存的事象と統計的変動が組み合わさり、相乗効果的に期待値が変化してる証左ともいえるかもしれない。
まとめ
プレイヤーの数を変化させて実験した。これは、依存的事象の検証ととらえることができる。その結果、プレイヤーの数を増やすごとに、マッチの本数の期待値には現象が見られた。ここから単純な理解をすると、依存的なステップ・工程が増えるごとに、期待値に低下がみられる。
また、サイコロの目の分布を変化させて実験した。これは、統計的変動の検証ととらえることができる。その結果、サイコロの目が1,6だけのようにばらつきが大きい場合、その期待値の減少が大きい。一方で、サイコロの目が3,4だけのときのようにばらつきが小さい場合、その期待値の減少は小さい。
物語では主人公はサイコロの期待値が3.5だから、その10倍で・・・という風にサイコロの振る舞いを固定した値として取り扱っていた。しかし、同じ期待値が3.5だとしても、サイコロの出目のパターンによって、最終的なマッチの本数は異なる。また、平均値だけの振る舞いであるが、プレイヤーの数が増えれば増えるほど、その差は大きくなっていくことが確認できた。
また、サイコロの目が2つだけのパターンではあるが、サイコロの目の分散が大きくなり、なおかつ、プレイヤーの数が増えると、期待値が相乗効果で下がっていくという傾向はみられた。これは、依存的事象と統計的変動が相互作用している例といえるかもしれない。
感想
ちょっと消化不良気味である。1,2,5,6の謎がいまいちわからない。それ以外を除けば、そこそこ納得いく答えではある。依存的事象と統計的変動が組み合わさったとき、「直感」というか雑な期待値による近似がほぼ役に立たないということが分かる例であった。
タスクの進行については、完了時期の見立てをする。ということをしたくなる。その時には、ベロシティの平均値がー残ったベロシティがーとやりたくはなるが、その予測はほぼ外れることが分かる。仮に当てたいとするのであれば、ベロシティを安定させれば、そこそこ終了時期を当てられる。ともいえる。
一方で、この記事を書いていて思い出した話がある。それは、システム移行に関するバグを直していたときのことである。システム移行の完了は期限が決められていた。とは、いうものの移行に関してバグが直らないものは直らんのである。では、いつ直るんだ。という話になる。その時、上司に詰められたので、じゃあ予測してやろう。と思った。
毎日、システム移行完了したユーザーが加算されていく。1日目は1人、2日目は3人、3日目は4人・・・と、そのデータを記録していた。そうすると、1日当たり移行完了した人数の確率分布を得る。そして、今日残っているユーザー数が分かっている。ということは、移行完了した人数の確率分布からランダム抽出を行い、残っているユーザー数から引く。その行為をユーザー数が0になるまで行う。このシミュレーションを何度も行うことにより、何日で移行完了可能か？という確率分布が計算できる。その計算結果をもって、「まぁ、今までの実績から考えて、締め切りまでに終わる可能性は〇〇%っスね」という話をした記憶がある。個人的には、非常に科学的で論理的な姿勢を示したつもりではあったが、「ぐだぐだ言ってないで早いこと直せ」と一蹴された記憶がある。わからんではない。
この話の面倒なところは科学リテラシーとビジネスの軋轢の問題もある。仮に、上司がマッチの本数の期待値が3.5本なんだから、10周やれば35本だ！！と言い張ったとして、いや、それ違いますよ。依存的事象と統計的変動により、そうはならんのですよ。と言ったところで、何も変わらんというのが事実である。おおむねその時には手遅れで、35本の仕事を取ってきてしまっており、現場が踏ん張ることになるのである。これは根本的には上司の科学リテラシーの無さからくるビジネス上のミスであり、相手先に謝り、その判断をした責任において上司の首を飛ばすのが本質的には正しいが、そうはならないし、そこで踏ん張ってしまって現場が乗り越えるとそれが実績となり、それがスタンダードになる。だんだん書いてて辛くなってきたな。
先ほどの私の移行のケースも、上司はなんらか、その上の上司と移行に関する締め切りの取り決めをしていたのであろう。それは上司なりの理屈をもって決めていたのかもしれとないが、部下が、それ論理的に考えて無理ですよ。と言われたところで、自分が論理的に間違っていたことを認めるのは難しい。しかも、その上の上司が理解したり、ビジネスとして容認できるかは別問題である。それだと、とりあえず反論するな。早く終わらせるために仕事しろ。が私に言うべき言葉の最適解になるのだろう。
下手の考え休みに似たり。とはいうが、下手が間違った推論と理論をもって行動した結果、周りに迷惑かけるという結果を引き起こしているので、下手には行動せず休んでもらった方が、周りが平和になるのかもしれない。
一方で、これは人間の認知バグのような気もする。これは、ある意味で、直感の概算が全然合わない。という話になる。ここがずれるとスケジュールが火の車になる。これを本質的にひかない、直感の方法はあるのだろうか・・・