モンテカルロ法 — 10 万回走らせて確率を数える

Zig で書いた勝率シミュレーションを、確率モデル・推定量・誤差・アルゴリズムの順に完全に分解します。「なぜ 10 万回なのか」「なぜ 0.11 なのか」「1 着の確率を式で書くとどうなるか」まで答えられるようになります。

このページで学ぶこと

1 頭の「出来」を「平均 + 揺らぎ」で表す確率モデル
脚質とペースの相性が 1 次関数で表されていること
「1 着になる確率」を式で書くと積分になり、手では解けないこと
モンテカルロ推定の誤差 $\sqrt{p(1-p)/N}$ と、10 万回の精度
2 頭だけなら $\Phi$ で正確に解けること (検算の方法)
揺らぎの係数 0.11 の決め方と、モデルの限界

確率モデル — 出来 = 平均 + 揺らぎ

1 回のレースで、 $i$ 番目の馬の「その日の出来」 $X_i$ を次のようにモデル化します。

X_i = \mu_i + \varepsilon_i

$\mu_i$ … その馬の 平均的な出来 (毎回同じ値)
$\varepsilon_i$ … その日の 揺らぎ (毎回違う乱数)

そして 出来がいちばん大きい馬が 1 着、2 番目が 2 着、… と決めます。とても単純ですが、心理学で 1927 年に提案された Thurstone (サーストン) モデル と同じ形で、「比べて順位をつける」現象の基本モデルとして今も使われています。

平均的な出来 $\mu_i$

native/zig/montecarlo.zig

const base = (speed[i] / 100.0) * styleMultiplier(style[i], pace) * (0.85 + 0.15 * stamina[i]);

\mu_i = \underbrace{\frac{\mathrm{SI}_i}{100}}_{\text{能力}} \times \underbrace{m(\text{脚質}_i, \pi)}_{\text{展開}} \times \underbrace{(0.85 + 0.15\,s_i)}_{\text{スタミナ}}

$\mathrm{SI}_i$ … スピード指数 (サンプル出馬表) または能力値から作った値 (能力値と適性)
$\pi$ … 想定ペース (0.8 = スロー、1.0 = ミドル、1.2 = ハイ)
$s_i$ … スタミナ (0〜1)。調教評価から S = 1.0、A+ = 0.9、A = 0.8、B = 0.65、C = 0.5 に換算。掛け算の項は 0.85〜1.0 の範囲になります

脚質とペースの相性 $m$ — 1 次関数

native/zig/montecarlo.zig

fn styleMultiplier(s: u32, pace: f32) f32 {
    return switch (s) {
        0 => 1.04 - (pace - 1.0) * 0.30,
        1 => 1.02 - (pace - 1.0) * 0.12,
        2 => 1.00 + (pace - 1.0) * 0.15,
        else => 0.98 + (pace - 1.0) * 0.30,
    };
}

$\delta = \pi - 1$ (ペースのミドルからのずれ) とすると、4 本の直線です。

m_\text{逃}(\delta) = 1.04 - 0.30\,\delta,\quad m_\text{先}(\delta) = 1.02 - 0.12\,\delta,\quad m_\text{差}(\delta) = 1.00 + 0.15\,\delta,\quad m_\text{追}(\delta) = 0.98 + 0.30\,\delta

ペース	逃げ	先行	差し	追込
スロー (0.8)	1.100	1.044	0.970	0.920
ミドル (1.0)	1.040	1.020	1.000	0.980
ハイ (1.2)	0.980	0.996	1.030	1.040

傾きの符号が意味を持っています。逃げ・先行は傾きが負 (ペースが速いほど不利)、差し・追込は正 (速いほど有利) です。2 本の直線が交わる点では、両者の有利不利がちょうど入れ替わります。たとえば逃げと追込は

1.04 - 0.30\,\delta = 0.98 + 0.30\,\delta \;\Longrightarrow\; \delta = 0.1 \;(\pi = 1.1)

で互角になります。ミドルよりやや速いペースが分かれ目、という設計です。

揺らぎ $\varepsilon_i$

native/zig/montecarlo.zig

fn gaussianish() f32 {
    return (nextRand() + nextRand() + nextRand() + nextRand()) - 2.0;
}
// ...
const score = base + gaussianish() * 0.11;

確率の基本で計算したとおり、 $\varepsilon_i$ は平均 0・標準偏差

\sigma = 0.11 \times \sqrt{4 \times \tfrac{1}{12}} \approx 0.0635

の、ほぼ正規分布です。各馬の揺らぎは独立 (別々の乱数) です。

実際の値 — 有馬記念 (サンプル)

馬	SI	脚質	調教	$\mu_i$ (ミドル)
ドウデュース	118.5	差し	S	1.1850
レガレイラ	116.0	差し	S	1.1600
スターズオンアース	114.2	先行	A+	1.1474
アーバンシック	115.2	追込	S	1.1290
ダノンデサイル	109.8	先行	A+	1.1032
ブローザホーン	113.4	差し	A	1.1000

16 頭の $\mu_i$ は 0.981〜1.185 に入り、差は最大 0.204 です。揺らぎの標準偏差 0.0635 と比べると「1〜3 標準偏差くらいの差」で、これが「強い馬が勝ちやすいが、番狂わせもある」バランスを生んでいます。

1 着になる確率を式で書く

$i$ 番目の馬が 1 着になる確率は「 $X_i$ が他のすべての $X_j$ より大きい確率」です。

p_i = P\!\left(X_i > X_j \text{ がすべての } j \ne i \text{ で成り立つ}\right)

揺らぎを正規分布とみなすと、「 $X_i$ がちょうど $x$ になり、他の全員が $x$ より小さい」確率を、すべての $x$ について足し合わせる (積分する) ことになります。

p_i = \int_{-\infty}^{\infty} \underbrace{\frac{1}{\sigma}\,\varphi\!\left(\frac{x - \mu_i}{\sigma}\right)}_{X_i \text{ が } x \text{ になる密度}} \;\prod_{j \ne i} \underbrace{\Phi\!\left(\frac{x - \mu_j}{\sigma}\right)}_{X_j < x \text{ の確率}} \, dx

( $\varphi$ は標準正規分布の密度、 $\prod$ は「全部掛ける」記号。他の馬は独立なので確率を掛け算できます)

この積分は手では解けない

$\Phi$ 自体が初等関数で書けない上に、それが 15 個掛かった積を積分するので、答えを式の形で求めることはできません。数値積分 (細かく区切って足す) なら 1 頭ずつ計算できますが、2 着・3 着まで含めると複雑になります。そこで「実際に何回も走らせて数える」モンテカルロ法 を使います。

モンテカルロ推定 — 数えて割る

$N$ 回レースを走らせ、 $i$ 番目の馬が 1 着になった回数を $W_i$ とすると、

\hat p_i = \frac{W_i}{N}

を $p_i$ の推定値とします。各回の「1 着かどうか」は確率 $p_i$ で 1、それ以外で 0 になる確率変数 (ベルヌーイ変数) なので、

E[\hat p_i] = p_i \quad (\text{偏りがない}), \qquad \sigma_{\hat p_i} = \sqrt{\frac{p_i(1 - p_i)}{N}}

です。 $N = 100{,}000$ 、 $p \approx 0.362$ (ドウデュース) なら

\sigma_{\hat p} = \sqrt{\frac{0.362 \times 0.638}{100000}} \approx 0.00152

つまり 標準誤差 0.15 ポイント、95% の確率で真の値は $\pm 1.96\sigma \approx \pm 0.30$ ポイントの範囲にあります。画面の「36.2%」は、小数第 1 位まで意味のある精度で表示していることになります。

実際の収束

アプリと同じシード (20261109) で $N$ だけを変えた、ドウデュースの勝率です。

$N$	勝率	標準誤差の理論値
100	38.00%	±4.8 pt
1,000	36.20%	±1.5 pt
10,000	36.02%	±0.48 pt
100,000	36.22%	±0.15 pt

下のデモは、アプリとまったく同じ Zig の WebAssembly を動かしています。 $N$ を変えて、誤差の幅 (黒い線) が $1/\sqrt N$ で縮むのを確かめてください。

2 頭だけなら正確に解ける — 検算の方法

上の積分は解けませんが、2 頭だけ なら簡単です。 $X_A - X_B = (\mu_A - \mu_B) + (\varepsilon_A - \varepsilon_B)$ で、揺らぎの差の分散は $\sigma^2 + \sigma^2 = 2\sigma^2$ (独立な分散は足し算)。よって

P(X_A > X_B) = P\!\left(Z > -\frac{\Delta\mu}{\sigma\sqrt2}\right) = \Phi\!\left(\frac{\Delta\mu}{\sigma\sqrt{2}}\right), \qquad \Delta\mu = \mu_A - \mu_B

ドウデュースとレガレイラは $\Delta\mu = 0.025$ なので

\Phi\!\left(\frac{0.025}{0.0635 \times 1.414}\right) = \Phi(0.278) \approx 0.6096

アプリと同じ乱数で 20 万回比べると 0.6083 で、差は 0.0013 (シミュレーションの標準誤差 0.0011 の約 1.2 倍で、偶然の誤差として十分ありうる大きさ。加えて、実際の揺らぎは正規分布ではなく Irwin–Hall 分布なので、ごくわずかなずれは理論的にも生じる)。モデルの実装が式どおりに動いている、という検算になります。

さらに深く「1 対 1 の確率」から 16 頭の勝率は出せない

「ドウデュースは各馬に 6〜9 割勝つ」から 16 頭立ての勝率を出そうとして、1 対 1 の確率を全部掛けるのは誤りです。 $X_A > X_B$ と $X_A > X_C$ は、どちらも同じ $X_A$ に依存しているので独立ではないからです ( $X_A$ が大きい日は、両方に同時に勝ちやすい)。だからこそ積分 (または模擬実験) が必要になります。

アルゴリズム — 上位 3 頭だけを追う

native/zig/montecarlo.zig

var top = [3]u32{ 0, 0, 0 };
var topScore = [3]f32{ -1e9, -1e9, -1e9 };
for (0..n) |idx| {
    const i: u32 = @intCast(idx);
    const base = (speed[i] / 100.0) * styleMultiplier(style[i], pace) * (0.85 + 0.15 * stamina[i]);
    const score = base + gaussianish() * 0.11;
    if (score > topScore[0]) {
        top[2] = top[1]; topScore[2] = topScore[1];
        top[1] = top[0]; topScore[1] = topScore[0];
        top[0] = i; topScore[0] = score;
    } else if (score > topScore[1]) {
        top[2] = top[1]; topScore[2] = topScore[1];
        top[1] = i; topScore[1] = score;
    } else if (score > topScore[2]) {
        top[2] = i; topScore[2] = score;
    }
}
wins[top[0]] += 1;
for (top[0..@min(n, 3)]) |p| places[p] += 1;

知りたいのは「1 着」と「3 着以内」だけなので、全頭を並べ替える必要はありません。暫定の上位 3 頭を覚えておき、それより良い馬が来たら下にずらして入れる (挿入ソートの一部) だけです。

1 回のレースの計算量は $O(n)$ (並べ替えると $O(n \log n)$ )
$N = 10^5$ 、 $n = 16$ なら、出来の計算が $1.6 \times 10^6$ 回、乱数が $6.4 \times 10^6$ 個
ブラウザでの実行時間は数十ミリ秒

3 着以内率は各回 3 頭が数えられるので、全頭を足すと 300% になります (テストで確認しています)。

揺らぎの係数 0.11 の決め方

0.11 は、揺らぎの大きさを決める唯一の自由なつまみです。小さすぎると「強い馬が必ず勝つ」、大きすぎると「誰が勝つか分からない」モデルになります。有馬記念 (サンプル) で係数を変えたときの、ドウデュースの勝率です。

係数	揺らぎの標準偏差	ドウデュースの勝率	最も弱い馬の勝率
0	0	100%	0%
0.03	0.017	80.8%	0%
0.05	0.029	62.0%	0%
0.11	0.064	36.3%	0.01%
0.20	0.115	22.1%	0.56%
0.50	0.289	11.8%	2.90%

実際の競馬では、1 番人気の勝率は 3 割前後です。0.11 は「スピード指数 101〜118 くらいの実在の出馬表で、1 番人気の勝率が 3〜4 割になる」ように手で調整した値です。

0.11 はデータから推定した値ではない

本来は、過去のレース結果を使って「実際の着順が最も起こりやすくなる係数」を探す (最尤推定) べきものです。このアプリは過去のレース結果データを持っていないため、手で合わせています。本物の予想へで、推定と評価の方法を紹介します。

モデルの限界

揺らぎが独立: 実際は「その日の馬場が内有利」のように、何頭かに共通して効く要因があります (相関のある揺らぎ)
ペースが固定: 実際のペースは、走ってみないと分からない確率変数です
出来の差がそのまま着順: 不利・出遅れ・故障のような「まれだが大きい」出来事が入っていません (Irwin–Hall の揺らぎには上限 $\pm 0.22$ があります)
能力の数字の誤差: スピード指数そのものの誤差は考えていません

それでも、このモデルは 全頭の確率の合計が必ず 100% になり、同じ条件なら同じ結果を返し、2 頭の場合は解析解と一致する、という「確率モデルとしての最低限の性質」を満たしています。これが予想エンジンの点数との大きな違いです。

まとめ

出来 $X_i = \mu_i + \varepsilon_i$ (Thurstone モデル)。 $\mu_i$ = 能力 × 展開 × スタミナ、 $\varepsilon_i$ は標準偏差 0.0635
脚質の有利不利はペースの 1 次関数。逃げと追込はペース 1.1 で互角
1 着の確率は積分で書けるが解けない → $N$ 回走らせて数える
誤差は $\sqrt{p(1-p)/N}$ 。10 万回で ±0.15 ポイント
2 頭なら $\Phi(\Delta\mu/(\sigma\sqrt2))$ で厳密に計算でき、検算に使える
0.11 は 1 番人気の勝率が 3〜4 割になるよう手で合わせた値

練習問題何回走らせればよいか

勝率 10% 前後の馬について、95% 信頼区間の幅を ±0.1 ポイント (0.001) 以内にしたいとき、何回走らせればよいですか。

答えを見る

$1.96\sqrt{0.1 \times 0.9 / N} \le 0.001$ を解くと $N \ge 1.96^2 \times 0.09 / 0.001^2 \approx 345{,}744$ 。約 35 万回です。今の 10 万回の 3.5 倍で、誤差は $\sqrt{3.5} \approx 1.87$ 分の 1 になります。

練習問題ハイペースでの逃げ馬

ハイペース ( $\pi = 1.2$ ) のとき、スピード指数 110・スタミナ 0.8 の逃げ馬と、スピード指数 106・スタミナ 0.8 の追込馬は、どちらの $\mu$ が大きいですか。

答えを見る

スタミナの項はどちらも $0.85 + 0.15 \times 0.8 = 0.97$ 。逃げ: $1.10 \times 0.98 \times 0.97 = 1.0457$ 。追込: $1.06 \times 1.04 \times 0.97 = 1.0693$ 。スピード指数が 4 低い追込馬の方が、平均的な出来は上になります。

練習問題2 頭の勝率から μ の差を逆算する

2 頭の対戦で強い方が 75% 勝つとき、 $\Delta\mu$ はいくつですか ( $\sigma = 0.0635$ 、 $\Phi(0.674) = 0.75$ を使ってよい)。

答えを見る

$\Delta\mu / (\sigma\sqrt2) = 0.674$ より $\Delta\mu = 0.674 \times 0.0635 \times 1.414 \approx 0.0605$ 。スピード指数 (他の条件が同じなら) に直すと約 6 の差です。

モンテカルロ法 — 10 万回走らせて確率を数える

#確率モデル — 出来 = 平均 + 揺らぎ

#平均的な出来 μi\mu_iμi​

#脚質とペースの相性 mmm — 1 次関数

#揺らぎ εi\varepsilon_iεi​

#実際の値 — 有馬記念 (サンプル)

#1 着になる確率を式で書く

#モンテカルロ推定 — 数えて割る

#実際の収束

#2 頭だけなら正確に解ける — 検算の方法

#アルゴリズム — 上位 3 頭だけを追う

#揺らぎの係数 0.11 の決め方

#モデルの限界

確率モデル — 出来 = 平均 + 揺らぎ

平均的な出来 $\mu_i$

脚質とペースの相性 $m$ — 1 次関数

揺らぎ $\varepsilon_i$

実際の値 — 有馬記念 (サンプル)

1 着になる確率を式で書く

モンテカルロ推定 — 数えて割る

実際の収束

2 頭だけなら正確に解ける — 検算の方法

アルゴリズム — 上位 3 頭だけを追う

揺らぎの係数 0.11 の決め方

モデルの限界