python 练习问题3.1コーシー·ベイズ推定量をモンテカルロ积分で计算する

Posted 2021-05-10

# 練習問題3.1 コーシー・ベイズ推定量をモンテカルロ積分で計算する
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cauchy, norm
import scipy.integrate

x = 4

# 分子の被積分関数
func1 = lambda t: t * norm(loc=x).pdf(t) * cauchy.pdf(t)

# 分母の被積分関数
func2 = lambda t: norm(loc=x).pdf(t) * cauchy.pdf(t)

# scipy.integrateで積分計算
I_nume = scipy.integrate.quad(func1, -Inf, Inf)[0]
I_denom = scipy.integrate.quad(func2, -Inf, Inf)[0]
I = I_nume / I_denom
print "scipy.integrate:", I

# モンテカルロ積分 (1)
# pをコーシー分布、fを正規分布とした場合
# コーシー分布からサンプリングする

x = 4
N = 100000

# 分子の積分をモンテカルロ法で計算
T = cauchy.rvs(size=N)
I_nume = np.mean(T * norm(loc=x).pdf(T))

# 分母の積分をモンテカルロ法で計算
T = cauchy.rvs(size=N)
I_denom = np.mean(norm(loc=x).pdf(T))

I = I_nume / I_denom
print u"モンテカルロ積分(1):", I

# モンテカルロ積分 (2)
# pを正規分布、fをコーシー分布とした場合
# 正規分布からサンプリングする

# 分子
T = norm(loc=x).rvs(size=N)
I_nume = np.mean(T * cauchy.pdf(T))

# 分母
T = norm(loc=x).rvs(size=N)
I_denom = np.mean(cauchy.pdf(T))

I = I_nume / I_denom
print u"モンテカルロ積分(2):", I