"""
ランダムな値の分配による指数分布（ボルツマン分布）への収束シミュレーション。
概要:
多数の参加者間でランダムにリソース（値）を分配し、システム全体の統計的な分布が
いかにして指数分布（Boltzmann分布）に近づくかをリアルタイムで可視化します。
詳細説明:
1. 初期状態では全員が同じ値（平均値）を持っています。
2. 各ステップで、参加者にランダムな値を割り振り直し、全体の総和が一定になるよう正規化します。
3. 反復計算を繰り返す中で、以下の3つのグラフを表示します：
- 個々の参加者の現在の値（ランダムな揺らぎ）。
- 値を昇順にソートしたプロファイル。
- 頻度分布（ヒストグラム）と理論的な指数分布 $f(x) = A e^{-x/\langle x \rangle}$ の比較。
4. 統計力学における「最大エントロピー」の状態がいかに指数分布として現れるかを体験的に理解できます。
関連リンク: :doc:`randomdistribution_usage`
"""
import sys
import random
import numpy as np
from math import exp
import matplotlib.pyplot as plt
#=============================
# デフォルトパラメータ
#=============================
NP_DEF = 50 # 参加人数
VAL_DEF = 50 # 初期値（および平均値）
NITER_DEF = 500 # 最大反復回数
NPLOT_DEF = 100 # グラフ更新間隔
NFX_DEF = 21 # 分布関数のメッシュ数
[ドキュメント]
def main():
"""
シミュレーションの主実行ルーチン。
"""
# 引数処理
argv = sys.argv
npersons = int(argv[1]) if len(argv) >= 2 else NP_DEF
avg_val = int(argv[2]) if len(argv) >= 3 else VAL_DEF
niter = int(argv[3]) if len(argv) >= 4 else NITER_DEF
nplotinterval = int(argv[4]) if len(argv) >= 5 else NPLOT_DEF
nfx = int(argv[5]) if len(argv) >= 6 else NFX_DEF
# 初期化
x_index = np.arange(npersons)
v = np.full(npersons, float(avg_val))
# 理論的な指数分布の計算
# f(x) = (N/avg) * exp(-x/avg)
beta = 1.0 / avg_val
coeff_a = npersons * beta
x_max_theory = avg_val * 10.0
ftx = np.linspace(0, x_max_theory, 101)
fty = [coeff_a * exp(-beta * x) for x in ftx]
print(f"\nSimulation started: n={npersons}, avg={avg_val}")
print(f"Theoretical distribution: f(x) = {coeff_a:.4g} * exp(-{beta:.4g} * x)")
# 描画の準備
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 10))
plt.subplots_adjust(hspace=0.4)
# シミュレーションループ
for i in range(niter):
# 1. ランダムな割り振り
raw_values = np.random.rand(npersons) * avg_val
# 2. 総和一定のための正規化
current_sum = np.sum(raw_values)
v = raw_values * (avg_val * npersons / current_sum)
# ソート（分布確認用）
v_sorted = np.sort(v)
# 指定間隔でのプロット更新
if i == 0 or (i + 1) % nplotinterval == 0 or i == niter - 1:
ax1.cla()
ax2.cla()
ax3.cla()
# グラフ1: 生の分配データ
ax1.bar(x_index, v, color='skyblue', alpha=0.7)
ax1.axhline(avg_val, color='red', linestyle='--', linewidth=0.8, label='Average')
ax1.set_title(f'Step {i+1}/{niter}: Random Resource Allocation')
ax1.set_ylabel('Value')
ax1.set_xlim(-1, npersons)
# グラフ2: ソートされたデータ（ローレンツ曲線に近い表現）
ax2.bar(x_index, v_sorted, color='steelblue', alpha=0.7)
ax2.axhline(avg_val, color='red', linestyle='--', linewidth=0.8)
ax2.set_title('Sorted Distribution Profile')
ax2.set_ylabel('Value (Sorted)')
ax2.set_xlim(-1, npersons)
# グラフ3: 分布関数（ヒストグラム）と理論曲線
counts, bins = np.histogram(v, bins=nfx, range=(0, max(v)), density=False)
bin_centers = (bins[:-1] + bins[1:]) / 2.0
bin_width = bins[1] - bins[0]
# 頻度を「値の幅」で割って分布関数にする
y_dist = counts / bin_width
ax3.step(bin_centers, y_dist, where='mid', label='Simulation', color='black')
ax3.plot(ftx, fty, label='Theory (Boltzmann)', color='red', linestyle='dashed')
ax3.set_title('Probability Distribution Function')
ax3.set_xlabel('Value')
ax3.set_ylabel('Frequency Density')
ax3.set_xlim(0, max(v) * 1.1)
ax3.legend()
plt.pause(0.05)
print("\nSimulation completed.")
plt.show(block=False)
input("Press ENTER to exit>>")
if __name__ == '__main__':
main()