アクチュアリーのためのPython入門（アクチュアリー数学編第4回）
シミュレーションの実例
📚 アクチュアリーのためのPython入門
この記事はアクチュアリー数学編の一部です。
▶ 目次はこちら
▶ 逆引きガイドはこちら
はじめに
アクチュアリー数学ではシミュレーションについて学習します。
理論が理解できても、結果数値を見るだけでは、
いまいち納得できない部分があります。
それをPythonで実際に計算してみたら、
どのような結果になるのかというのが、ここでの内容です。
Excelでも十分にできる内容ですが、
Pythonの勉強のため、アクチュアリー数学の基本的なものを実践していきます。
中心極限定理
まずは基本的な中心極限定理についてです。
標本を大きくとると正規分布に近づいていくという理論です。
まず、区間(0,1)の一様乱数を12個作成します。
それを合計して6を引きます。
この結果が理論的には、
期待値が0、分散が1になるので、
これを繰り返した結果は標準正規分布に近づくはずです。
それをPythonで1000回繰り返してみましょう。
コードは次のとおりです。
# 中心極限定理
import numpy as np
np.random.seed(42)
# 計算結果のリスト
NormRmd = []
for i in range(1000):
total = 0
for j in range(12):
# 一様乱数U(0,1)
total += np.random.rand()
# 平均を引く
x = total - 6
NormRmd.append(x)
NormRmdに結果がリスト化していますが、
これをヒストグラムにすると、次のようになります。
少しゆがみが残りますが、もっと標本を増やせば
標準正規分布に近づきます。
棄却法
棄却法によって別の確率分布に従う乱数を作成できますが、
理論的には分かるものの、分かりにくい部分があります。
実際にはどうなるのか再現してみます。
アクチュアリー試験の数学の2024年度問題3(3)を再現します。
指数分布の乱数と一様分布の乱数から、
正の標準正規分布の乱数を作成します。
1000回繰り返します。
コードは次のとおりです。
# 棄却法
import numpy as np
np.random.seed(42)
# 計算結果のリスト
NormRmd = []
for i in range(1000):
U1, U2 = np.random.rand(2)
# 指数分布の乱数を逆関数法で作成
y = - np.log(1 - U1)
# f/cg >= u なら採択
if U2 <= np.exp(-0.5 * ((y - 1) ** 2)):
NormRmd.append(y)
リストには743個の数値が残りますが、
$$ \frac{1}{c} = \sqrt{\frac{\pi}{2e}} \fallingdotseq 0.76 $$
なので、理論値には近いです。
この結果をグラフに出力すると次のとおりです。
求める確率分布には近いでしょうか。
極座標法
一様分布の乱数から標準正規分布の乱数を作成する方法として、
他に極座標法がありますので、それも実践してみます。
2つの一様乱数から、正規分布の乱数を作成することができるので、
効率がいいです。
詳しくはアクチュアリー会のモデリングのテキストにありますので、
そちらを参照ください。
コードは次のとおりです。
# 極座標法
import numpy as np
np.random.seed(42)
NormRmd = []
for i in range(1000):
U1, U2, U3 = np.random.rand(3)
V1 = 2 * U1 - 1
V2 = 2 * U2 - 1
S = V1 ** 2 + V2 ** 2
if S <= 1 and S > 0:
sqS = (-2 * np.log(S) / S) ** 0.5
x = sqS * V1
y = sqS * V2
NormRmd.extend([x, y])
結果をグラフにしてみると次のとおりです。
まとめ
今回はアクチュアリー数学で勉強するシミュレーションのうち、
中心極限定理、棄却法、極座標法を使って、
実際に目的の関数にどれくらい近づくのか視覚化してみました。
サンプルでいいので計算してみることで、
シミュレーションンへの理解が深まるものと思います。
（補足）
今回、logとexpはnumpyの関数を使用しました。
乱数リストの対数と指数を計算するので、
numpyの方が計算しやすいという理由です。
単一の数値を計算するなら、
一般的にはmathのlogとexpを使用します。
📚 ナビゲーション
◀ 前の記事
予定死亡率の補整３
▶ 次の記事
複合分布
📚 目次
アクチュアリーのためのPython入門