はじめに
単位立方体上の2点をランダムに選び、その間の距離を考える。その距離の平均値（期待値）は、解析的に0.66程度に収束することが知られている。これをロビンス定数という。そこで、今回は計算機を用いたモンテカルロ法でロビンス定数の近似値を算出する。
プログラム
・単位立方体上にnumpy randomによって生成されたランダムな2点を作成する。
・その距離を三平方の定理により計算する。
以下、プログラムを示す。
python Robbins_constant_cal.py
"""
ロビンス定数 (Robin's Constant) のモンテカルロシミュレーション
============================================================
【ロビンス定数とは】
単位立方体 [0,1]^3 内にランダムに選んだ2点間の距離の期待値。
1978年にDavid Robbinsが厳密な解析解を求めた。
【解析的な値（厳密解）】
Δ = (4 + 17√2 - 6√3 + 21·ln(1+√2) + 42·ln(2+√3) - 7π) / 105
≈ 0.66170...
【理論的背景】
- 単位区間 [0,1] の1次元版: 期待値 = 1/3 ≈ 0.3333
- 単位正方形 [0,1]^2 の2次元版: 期待値 ≈ 0.5214
- 単位立方体 [0,1]^3 の3次元版: 期待値 ≈ 0.6617 ← このプログラムが扱う値
このプログラムでは:
1. ランダムな線分の距離分布をヒストグラムで可視化
2. 試行回数を増やすと期待値がロビンス定数に収束する様子を確認
3. 単位立方体内の線分を3Dで可視化
"""
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
from mpl_toolkits.mplot3d import Axes3D
# ─────────────────────────────────────────
# ロビンス定数の厳密解を計算
# ─────────────────────────────────────────
def robbins_constant_exact():
"""
単位立方体内の2点間距離の期待値（ロビンス定数）の厳密解を返す。
解析解:
Δ = (4 + 17√2 - 6√3 + 21·ln(1+√2) + 42·ln(2+√3) - 7π) / 105
参考: D. Robbins (1978), "Average distance between two points in a box"
"""
sqrt2 = math.sqrt(2)
sqrt3 = math.sqrt(3)
numerator = (
4
+ 17 * sqrt2
- 6 * sqrt3
+ 21 * math.log(1 + sqrt2)
+ 42 * math.log(2 + sqrt3)
- 7 * math.pi
)
return numerator / 105
ROBBINS_CONST = robbins_constant_exact()
# ─────────────────────────────────────────
# モンテカルロ計算
# ─────────────────────────────────────────
def sample_random_segments(n, rng=None):
"""
単位立方体内でランダムな線分を n 本サンプリングし、各線分の長さを返す。
Parameters
----------
n : int — サンプル数
rng : np.random.Generator (省略可)
Returns
-------
x1, y1, z1 : 始点座標の配列 (各長さ n)
x2, y2, z2 : 終点座標の配列 (各長さ n)
d : 線分の長さの配列 (長さ n)
"""
if rng is None:
rng = np.random.default_rng()
x1, y1, z1 = rng.random(n), rng.random(n), rng.random(n)
x2, y2, z2 = rng.random(n), rng.random(n), rng.random(n)
d = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1 - z2) ** 2)
return x1, y1, z1, x2, y2, z2, d
def estimate_robbins_constant(num_samples, rng=None):
"""
モンテカルロ法でロビンス定数を推定する。
Parameters
----------
num_samples : int — サンプル数（多いほど精度が上がる）
Returns
-------
estimate : float — 推定値
std_err : float — 標準誤差
"""
_, _, _, _, _, _, d = sample_random_segments(num_samples, rng)
return np.mean(d), np.std(d) / np.sqrt(num_samples)
# ─────────────────────────────────────────
# メインの可視化処理
# ─────────────────────────────────────────
NUM_SEGMENTS = 1000 # 可視化に使うサンプル数
NUM_TRIALS = 500 # 収束グラフの試行ステップ数
rng = np.random.default_rng(seed=42) # 再現性のためシードを固定
# ── 厳密解を表示 ──────────────────────────
print("=" * 50)
print(f"ロビンス定数（厳密解）: {ROBBINS_CONST:.8f}")
estimate, stderr = estimate_robbins_constant(NUM_SEGMENTS * 10, rng)
print(f"モンテカルロ推定値 : {estimate:.8f} ±{stderr:.6f} (n={NUM_SEGMENTS * 10})")
print(f"誤差 : {abs(estimate - ROBBINS_CONST):.2e}")
print("=" * 50)
# ── 図1: 距離の分布（ヒストグラム） ──────
x1, y1, z1, x2, y2, z2, d = sample_random_segments(NUM_SEGMENTS, rng)
fig1, ax1 = plt.subplots(figsize=(8, 5))
ax1.hist(d, bins=30, color="steelblue", edgecolor="white", alpha=0.85, label="距離の分布")
ax1.axvline(ROBBINS_CONST, color="crimson", linewidth=2,
linestyle="--", label=f"ロビンス定数 Δ ≈ {ROBBINS_CONST:.4f}")
ax1.axvline(np.mean(d), color="orange", linewidth=2,
linestyle="-.", label=f"サンプル平均 = {np.mean(d):.4f}")
ax1.set_xlabel("2点間の距離")
ax1.set_ylabel("度数")
ax1.set_title(f"単位立方体内のランダムな2点間距離の分布 (n={NUM_SEGMENTS})")
ax1.legend()
plt.tight_layout()
plt.savefig(f"ロビンス定数_分布_{NUM_SEGMENTS}.png", dpi=150)
plt.show()
# ── 図2: 収束グラフ ───────────────────────
# 試行回数を 1 → NUM_TRIALS と増やしながら、距離の平均がロビンス定数に収束する様子を観察する
trial_indices = np.arange(1, NUM_TRIALS + 1)
d_ave_array = np.zeros(NUM_TRIALS)
for i in range(NUM_TRIALS):
_, _, _, _, _, _, d_i = sample_random_segments(i + 1, rng)
d_ave_array[i] = np.mean(d_i)
fig2, ax2 = plt.subplots(figsize=(9, 5))
ax2.plot(trial_indices, d_ave_array, color="steelblue", linewidth=1.0,
alpha=0.9, label="サンプル平均（モンテカルロ）")
ax2.axhline(ROBBINS_CONST, color="crimson", linewidth=2,
linestyle="--", label=f"ロビンス定数 Δ ≈ {ROBBINS_CONST:.5f}")
# 誤差帯：標準誤差 ±2σ（95%信頼区間の目安）
sigma_band = 2 * np.std(d) / np.sqrt(trial_indices)
ax2.fill_between(trial_indices,
ROBBINS_CONST - sigma_band,
ROBBINS_CONST + sigma_band,
color="crimson", alpha=0.12, label="±2σ 信頼帯（目安）")
ax2.set_xlabel("試行回数（サンプルサイズ）")
ax2.set_ylabel("線分の長さの平均")
ax2.set_title("モンテカルロ推定値のロビンス定数への収束")
ax2.legend()
plt.tight_layout()
plt.savefig(f"ロビンス定数_試行回数_{NUM_TRIALS}.png", dpi=150)
plt.show()
# ── 図3: 3D 可視化 ────────────────────────
# 最後のサンプルバッチをそのまま使って単位立方体内の線分を描画する
x1, y1, z1, x2, y2, z2, d = sample_random_segments(NUM_SEGMENTS, rng)
# 線分の長さに応じて色付け（短い=青、長い=赤）
norm = plt.Normalize(d.min(), d.max())
cmap = plt.cm.coolwarm
fig3 = plt.figure(figsize=(9, 7))
ax3 = fig3.add_subplot(111, projection="3d")
for i in range(NUM_SEGMENTS):
color = cmap(norm(d[i]))
ax3.plot([x1[i], x2[i]], [y1[i], y2[i]], [z1[i], z2[i]],
color=color, linewidth=0.5, alpha=0.5)
# カラーバーを追加
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
fig3.colorbar(sm, ax=ax3, shrink=0.6, label="線分の長さ")
ax3.set_xlabel("X")
ax3.set_ylabel("Y")
ax3.set_zlabel("Z")
ax3.set_title(f"単位立方体内のランダム線分 (n={NUM_SEGMENTS})\n"
f"平均長さ ≈ {np.mean(d):.4f} | 理論値（ロビンス定数）≈ {ROBBINS_CONST:.4f}")
plt.tight_layout()
plt.savefig(f"ロビンス定数_図示_{NUM_SEGMENTS}.png", dpi=150)
plt.show()
結果
実験の様子
試行回数が1000回の場合の実験の様子を以下に示す。
距離のヒストグラム
試行回数が1000回の場合のヒストグラムを以下に示す。
このように、かなり正規分布に近づく。
試行回数
距離の期待値の収束の様子を以下に示す。
ただし、試行回数を500まで行う。
このように、試行回数を増やすにつれて、解析解へ近づく。
まとめ
今回は、モンテカルロ法を用いてロビンス定数を推定した。ロビンス定数は、厳密な理論解析解が知られている。しかし、その算出はとても難しいので、計算機を用いて近似解を求め比較を行った。
参考文献