10. Python を使った数値シミュレーション
前回は統計の基礎について学びました。今回は、その知識を活かしてコンピュータでシミュレーションを行う方法を学んでいきましょう。
目次
- 10.1 シミュレーションとは?やり方とその意義
- 10.2 乱数の基礎と random モジュールの使い方
- 10.3 確率分布とは(一様分布・二項分布・正規分布)
- 10.4 大数の法則:試行回数を増やすと平均値が一定に近づく
- 10.5 モンテカルロ法:乱数を使って円周率を求めてみよう
- 10.6 色々なシミュレーションをしてみよう
- Ex.1 コインを 100 回投げて表と裏の回数を数えるプログラム
- Ex.2 サイコロを繰り返し振って、各目の出る確率が 1/6 に近づく様子をグラフ化するプログラム
この講座で使用する Google Colab の URL
10. Python を使った数値シミュレーション (Google Colab)
演習課題
Ex.10. Python を使った数値シミュレーション (Google Colab)
この講座で使用する Python, Jupyter Notebook のファイルと実行環境
Lesson 10: high-school-python-code (GitHub)
10.1 シミュレーションとは?やり方とその意義
シミュレーションとは、現実の事象をコンピュータで再現することです。例えば、サイコロを 1000 回振ったらどうなるか、雨の日が続く確率はどのくらいかなど、実際に試すのが難しい場合でもコンピュータで簡単に実験できます。
シミュレーションの基本例
import random
# サイコロを100回振るシミュレーション
results = []
for _ in range(100):
dice = random.randint(1, 6)
results.append(dice)
# 結果の集計
for number in range(1, 7):
count = results.count(number)
percentage = count / len(results) * 100
print(f"{number}の目: {count}回 ({percentage:.1f}%)")
10.2 乱数の基礎と random モジュールの使い方
random モジュールを使って、様々な種類の乱数を生成できます:
random モジュールの基本的な使い方
import random
# 整数の乱数(1から6まで)
dice = random.randint(1, 6)
# 0から1までの小数
probability = random.random()
# 範囲を指定した小数
temperature = random.uniform(20.0, 30.0)
# リストからランダムな要素を選択
colors = ["赤", "青", "緑", "黄"]
chosen = random.choice(colors)
# リストをシャッフル
cards = list(range(1, 53))
random.shuffle(cards)
10.3 確率分布とは(一様分布・二項分布・正規分布)
代表的な確率分布をシミュレーションで体験します:
様々な確率分布のシミュレーション
import random
import matplotlib.pyplot as plt
import numpy as np
# 一様分布(サイコロのような、全ての値が同じ確率で出る)
uniform_numbers = [random.uniform(0, 1) for _ in range(1000)]
# 二項分布(コイン投げを10回試行する実験を1000回)
def coin_flip_experiment():
flips = [random.choice(['表', '裏']) for _ in range(10)]
return flips.count('表')
binomial_results = [coin_flip_experiment() for _ in range(1000)]
# 正規分布(身長や体重のような、平均値を中心に分布する)
normal_numbers = np.random.normal(170, 5, 1000) # 平均170、標準偏差5
# ヒストグラムの描画
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.hist(uniform_numbers, bins=30)
plt.title('一様分布')
plt.subplot(1, 3, 2)
plt.hist(binomial_results, bins=11)
plt.title('二項分布')
plt.subplot(1, 3, 3)
plt.hist(normal_numbers, bins=30)
plt.title('正規分布')
plt.tight_layout()
plt.show()
10.4 大数の法則:試行回数を増やすと平均値が一定に近づく
コイン投げを例に、大数の法則を確認します:
大数の法則のシミュレーション
import random
import matplotlib.pyplot as plt
def coin_flip_simulation(trials):
heads = 0
heads_ratio = []
for i in range(1, trials + 1):
if random.random() < 0.5: # コインを投げる
heads += 1
heads_ratio.append(heads / i) # 表の割合を記録
return heads_ratio
# 異なる試行回数でシミュレーション
trials_1000 = coin_flip_simulation(1000)
trials_10000 = coin_flip_simulation(10000)
# グラフの描画
plt.figure(figsize=(12, 6))
plt.plot(trials_1000, label='1000回試行')
plt.plot(trials_10000, label='10000回試行')
plt.axhline(y=0.5, color='r', linestyle='--', label='理論値(0.5)')
plt.xlabel('試行回数')
plt.ylabel('表の出る確率')
plt.title('コイン投げシミュレーション:大数の法則の確認')
plt.legend()
plt.grid(True)
plt.show()
10.5 モンテカルロ法:乱数を使って円周率を求めてみよう
モンテカルロ法を使って円周率を計算します:
モンテカルロ法で円周率を計算
import random
import matplotlib.pyplot as plt
def estimate_pi(points):
inside_circle = 0
x_inside, y_inside = [], []
x_outside, y_outside = [], []
for _ in range(points):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
# 点が単位円の中にあるか判定
if x*x + y*y <= 1:
inside_circle += 1
x_inside.append(x)
y_inside.append(y)
else:
x_outside.append(x)
y_outside.append(y)
# 円周率の計算
pi_estimate = 4 * inside_circle / points
return pi_estimate, x_inside, y_inside, x_outside, y_outside
# シミュレーション実行
points = 10000
pi_estimate, x_in, y_in, x_out, y_out = estimate_pi(points)
# 結果の表示
print(f"計算された円周率: {pi_estimate}")
print(f"実際の円周率との誤差: {abs(pi_estimate - 3.14159265359)}")
# 点のプロット
plt.figure(figsize=(8, 8))
plt.scatter(x_in, y_in, c='blue', s=1, alpha=0.6, label='円の内側')
plt.scatter(x_out, y_out, c='red', s=1, alpha=0.6, label='円の外側')
plt.axis('equal')
plt.title(f'モンテカルロ法による円周率の計算\n推定値: {pi_estimate:.6f}')
plt.legend()
plt.grid(True)
plt.show()
10.6 色々なシミュレーションをしてみよう
実践的なシミュレーション例を見てみましょう:
じゃんけんのシミュレーション
import random
from collections import Counter
def play_janken():
hands = ['グー', 'チョキ', 'パー']
rules = {
('グー', 'チョキ'): 'グー',
('チョキ', 'パー'): 'チョキ',
('パー', 'グー'): 'パー',
}
player1 = random.choice(hands)
player2 = random.choice(hands)
if player1 == player2:
return '引き分け'
# 勝敗判定
if (player1, player2) in rules:
if rules[(player1, player2)] == player1:
return 'プレイヤー1の勝ち'
else:
if rules[(player2, player1)] == player2:
return 'プレイヤー2の勝ち'
# 1000回のじゃんけんをシミュレーション
results = []
for _ in range(1000):
result = play_janken()
results.append(result)
# 結果の集計
counter = Counter(results)
for outcome, count in counter.items():
percentage = count / len(results) * 100
print(f"{outcome}: {count}回 ({percentage:.1f}%)")
10.7 グラフで結果を可視化して考察する
シミュレーション結果を視覚的に表現します:
サイコロの期待値のシミュレーション
import random
import matplotlib.pyplot as plt
import numpy as np
def dice_simulation(trials):
# サイコロを振って平均値を記録
averages = []
running_sum = 0
for i in range(1, trials + 1):
roll = random.randint(1, 6)
running_sum += roll
average = running_sum / i
averages.append(average)
return averages
# 複数回のシミュレーション
num_simulations = 5
trials = 1000
plt.figure(figsize=(12, 6))
# 各シミュレーションを実行してプロット
for i in range(num_simulations):
averages = dice_simulation(trials)
plt.plot(averages, alpha=0.5, label=f'シミュレーション{i+1}')
# 理論値(3.5)の線を追加
plt.axhline(y=3.5, color='r', linestyle='--', label='理論値(3.5)')
plt.xlabel('試行回数')
plt.ylabel('平均値')
plt.title('サイコロの平均値の収束')
plt.legend()
plt.grid(True)
plt.show()
Ex.1 コインを 100 回投げて表と裏の回数を数えるプログラム
コイン投げシミュレーション
import random
from collections import Counter
import matplotlib.pyplot as plt
def coin_flip_experiment(flips):
# コインを指定回数投げる
results = []
for _ in range(flips):
result = 'タ' if random.random() < 0.5 else 'オ'
results.append(result)
# 結果を集計
count = Counter(results)
# 結果を表示
print(f"\n{flips}回コインを投げた結果:")
for face, num in count.items():
percentage = (num / flips) * 100
print(f"{face}: {num}回 ({percentage:.1f}%)")
# グラフの作成
plt.figure(figsize=(8, 6))
plt.bar(count.keys(), count.values())
plt.title(f'コイン投げの結果 ({flips}回)')
plt.ylabel('回数')
plt.show()
return count
# プログラムの実行
result = coin_flip_experiment(100)
# チャレンジ:
# - 複数回実験を行い、結果の平均を取る
# - 試行回数を変えて実験する
# - アニメーションで結果をリアルタイム表示する
Ex.2 サイコロを繰り返し振って、各目の出る確率が 1/6 に近づく様子をグラフ化するプログラム
サイコロの確率収束シミュレーション
import random
import matplotlib.pyplot as plt
import numpy as np
def dice_probability_simulation(trials):
# 結果を記録する配列
counts = {i: 0 for i in range(1, 7)}
probabilities = {i: [] for i in range(1, 7)}
# サイコロを振るシミュレーション
for i in range(1, trials + 1):
roll = random.randint(1, 6)
counts[roll] += 1
# 各目の確率を計算して記録
for num in range(1, 7):
prob = counts[num] / i
probabilities[num].append(prob)
# グラフの描画
plt.figure(figsize=(12, 6))
# 各目の確率の推移をプロット
for num in range(1, 7):
plt.plot(probabilities[num],
label=f'{num}の目',
alpha=0.7)
# 理論値(1/6)の線を追加
plt.axhline(y=1/6, color='r', linestyle='--',
label='理論値(1/6)')
plt.xlabel('試行回数')
plt.ylabel('確率')
plt.title('サイコロの確率収束シミュレーション')
plt.legend()
plt.grid(True)
plt.show()
# 最終的な確率を表示
print("\n最終的な確率:")
for num in range(1, 7):
final_prob = probabilities[num][-1]
print(f"{num}の目: {final_prob:.3f}")
# プログラムの実行
dice_probability_simulation(1000)
# チャレンジ:
# - 複数のサイコロを同時に投げる
# - 目の合計の分布を調べる
# - 特定の目が出る確率を予測する
まとめ
この章で学んだこと:
- シミュレーションの基本的な考え方
- random モジュールの使い方
- 様々な確率分布の特徴
- 大数の法則の実践的な理解
- モンテカルロ法による数値計算
- データの可視化と結果の考察方法
次の章では、Pythonの便利な活用方法について学び、より実践的なプログラミングスキルを身につけていきます!