Ccmmutty logo
Commutty IT
・0 pv・2 min read

NumPyで試す計算機統計学

Python3統計学モンテカルロ法NumPy
https://cdn.magicode.io/media/notebox/8d5d19ae-d460-4f93-8137-b3b7dde0736b.jpeg

やったこと

  1. 下記の書籍の例5.1 「単純モンテカルロ積分」と例5.2 「単純モンテカルロ積分（続き）」のコードをPython3で書き直した。

例5.1 「単純モンテカルロ積分」

モンテカルロ法により推定する積分の式

この問題で近似する積分の式は以下の通りである。この式でやっていることは、00から11の一様分布に従う確率変数XXに対して、期待値E(e−X)E(e^{-X})を計算しているのと同じである。
θ=∫01e−xdx\theta = \int_{0}^{1} e^{-x} dx
実際の答えは、1−e−11-e^{-1}でおよそ0.63210.6321程度になる。
python
# ライブラリのインポート
import numpy as np
python
m = 10000
x = np.random.rand(m)
theta_hat = np.mean(np.exp(-x))
# 推定値
print(theta_hat)
# 真の値
print(1 - np.exp(-1))
0.6323103709213596 0.6321205588285577

例5.2 「単純モンテカルロ積分（続き）」

モンテカルロ法により推定する積分の式

この問題で近似する積分の式は以下の通りである。この式でやっていることは、22から44の一様分布に従う確率変数XXに対して、期待値E(e−X)E(e^{-X})を計算しているのと同じである。例5.1の場合と違って確率密度関数の12\frac{1}{2}が考慮されていないので平均値をとったものに4−2=24-2=2倍する必要がある。
θ=∫24e−xdx\theta = \int_{2}^{4} e^{-x} dx
実際の答えは、e−2−e−4e^{-2}-e^{-4}でおよそ0.11700.1170程度になる。
python
m = 10000
x = np.random.rand(m) * 2 + 2
theta_hat = np.mean(np.exp(-x)) * (4 - 2)
# 推定値
print(theta_hat)
# 真の値
print(np.exp(-2) - np.exp(-4))
0.1167710927924694 0.11701964434787852

Discussion

コメントにはログインが必要です。