コイントスのギャンブルのやつ

gistfile1.txt

"""
ギャンブルのモンテカルロシミュレーション

ルール:
  - 表/裏が出る確率はそれぞれ 50%
  - 現在資金の n% を賭ける
  - 表が出たら、賭けた分が 1.5 倍になる
  - 裏が出たら、賭けた分が 0.6 倍になる
    （n=100% のとき、資産全体が 1.5 倍 / 0.6 倍 になる）

  1ゲーム = 100 回プレイ、それを 1000 回繰り返す（モンテカルロ試行）。
  n を 0%〜100% まで 1% 刻みで振り、最終資産の
  平均値・中央値・5%タイル・95%タイル を箱ひげ図で表示する。
"""

import numpy as np
import matplotlib
matplotlib.use("Agg")  # ウィンドウを開かず PNG に保存する
import matplotlib.pyplot as plt

# ---- 日本語フォント（環境に無ければ無視される） ----
plt.rcParams["font.family"] = ["MS Gothic", "Yu Gothic", "Meiryo", "sans-serif"]
plt.rcParams["axes.unicode_minus"] = False

# ---- パラメータ ----
W0 = 1000.0         # 初期資金
N_ROUNDS = 1000     # 1ゲームあたりのプレイ回数
N_SIMS = 1000       # モンテカルロ試行回数
UP, DOWN = 1.5, 0.6 # 賭けた分の倍率（表 / 裏）
fractions = np.arange(0, 101) / 100.0  # 0%〜100% を 1% 刻み

rng = np.random.default_rng(42)

# ---- シミュレーション本体 ----
# 毎ラウンド「その時点の現在資金 wealth の f 割」を賭ける（複利／固定比率ベット）。
#   stake  = f * wealth                       … 賭ける額（現在資金に比例）
#   表: wealth = (wealth - stake) + stake * UP    = wealth * (1 + f*(UP-1))
#   裏: wealth = (wealth - stake) + stake * DOWN  = wealth * (1 - f*(1-DOWN))
#
# 全 n で同じコインの出目を共有し、公平に比較できるようにする。
# coin[s, t] : 試行 s のラウンド t が表なら True
coin = rng.random((N_SIMS, N_ROUNDS)) < 0.5

results = {}  # f -> 最終資産の配列 (N_SIMS,)
for f in fractions:
    up_mult = 1.0 + f * (UP - 1.0)      # 表のときの資産倍率
    down_mult = 1.0 - f * (1.0 - DOWN)  # 裏のときの資産倍率
    wealth = np.full(N_SIMS, W0)        # 各試行の現在資金（初期資金からスタート）
    for t in range(N_ROUNDS):
        # stake = f * wealth を賭け、表/裏に応じて現在資金を更新
        wealth = np.where(coin[:, t], wealth * up_mult, wealth * down_mult)
    results[f] = wealth

# ---- 統計量の計算 ----
percents = fractions * 100
mean = np.array([results[f].mean() for f in fractions])
median = np.array([np.median(results[f]) for f in fractions])
p05 = np.array([np.percentile(results[f], 5) for f in fractions])
p95 = np.array([np.percentile(results[f], 95) for f in fractions])

# ---- 箱ひげ図（箱 = 5%〜95%、中央線 = 中央値、点 = 平均値） ----
fig, ax = plt.subplots(figsize=(16, 7))

bxp_stats = []
for i, f in enumerate(fractions):
    bxp_stats.append({
        "whislo": p05[i], "q1": p05[i],
        "med": median[i],
        "q3": p95[i], "whishi": p95[i],
        "mean": mean[i],
        "label": f"{int(round(f * 100))}",
    })

ax.bxp(
    bxp_stats,
    showmeans=True, meanline=False, showfliers=False,
    positions=percents,
    widths=0.7,
    boxprops=dict(color="steelblue"),
    medianprops=dict(color="crimson", linewidth=1.5),
    meanprops=dict(marker="D", markerfacecolor="orange",
                   markeredgecolor="orange", markersize=4),
    whiskerprops=dict(color="steelblue"),
    capprops=dict(color="steelblue"),
)

ax.set_yscale("log")
ax.set_xlabel("賭ける割合 n (%)")
ax.set_ylabel("最終資産（初期資金=1000、対数スケール）")
ax.set_title(
    f"最終資産の分布（{N_ROUNDS}回プレイ × {N_SIMS}試行）\n"
    "箱=5%〜95%タイル / 赤線=中央値 / 橙◆=平均値"
)
ax.axhline(W0, color="gray", linestyle="--", linewidth=1, label="初期資金（±0）")

# x 軸の目盛りを 10% ごとに間引き
tick_idx = np.arange(0, 101, 10)
ax.set_xticks(percents[tick_idx])
ax.set_xticklabels([f"{p:.0f}" for p in percents[tick_idx]])

ax.legend(loc="upper left")
ax.grid(True, which="both", axis="y", alpha=0.3)
fig.tight_layout()

# ---- 参考: 中央値・平均値が最大となる n を表示 ----
best_median_n = percents[np.argmax(median)]
best_mean_n = percents[np.argmax(mean)]
print(f"中央値が最大となる n = {best_median_n:.0f}%  (中央値 = {median.max():.1f}  = {median.max()/W0:.3f} 倍)")
print(f"平均値が最大となる n = {best_mean_n:.0f}%  (平均値 = {mean.max():.3g}  = {mean.max()/W0:.3g} 倍)")

fig.savefig("gamble_montecarlo.png", dpi=120)
print("グラフを gamble_montecarlo.png に保存しました。")