In [6]:

    

# 任意の関数 f(x) の [a, b] での積分をモンテカルロ法で求める
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform
import scipy.integrate

a, b = -1, 1
f = lambda x: x ** 2

# 被積分関数をプロット
x = np.linspace(-1.5, 1.5, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)

# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 0.666666666667
モンテカルロ積分: 0.667585521897






    

In [11]:

    

a, b = 0, 1
f = lambda x: x * np.exp(x**2)

# 被積分関数をプロット
x = np.linspace(0, 2, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)

# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 0.85914091423
モンテカルロ積分: 0.855417623185






    

In [8]:

    

a, b = 0, 1
f = lambda x: x**2 * np.exp(2*x)

# 被積分関数をプロット
x = np.linspace(0, 1.2, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)

# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 1.59726402473
モンテカルロ積分: 1.60222532927






    

In [52]:

    

a, b = 1, np.exp(1)
f = lambda x: np.log(x) / x

# 被積分関数をプロット
x = np.linspace(1, 5, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)

# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 0.5
モンテカルロ積分: 0.499921834232






    

In [59]:

    

a, b = 1, 2
f = lambda x: x**3 * np.log(x)

# 被積分関数をプロット
x = np.linspace(1, 3, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)

# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 1.83508872224
モンテカルロ積分: 1.8328253236






    

In [65]:

    

a, b = np.exp(1), np.exp(2)
f = lambda x: np.log(x) ** 2

# 被積分関数をプロット
x = np.linspace(2, 8, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
gca().add_patch(poly)
                
# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 12.0598303694
モンテカルロ積分: 12.0453421225






    

In [5]:

    

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform
import scipy.integrate

# 例3.3

a, b = 0, 1
f = lambda x: (np.cos(50 * x) + np.sin(20 * x)) ** 2

# 被積分関数をプロット
x = np.linspace(-0.2, 1.2, 1000)
y = f(x)
plt.plot(x, y)

# 積分範囲を色付け
ix = np.arange(a, b, 0.001)
iy = f(ix)
verts = [(a, 0)] + list(zip(ix, iy)) + [(b,0)]
poly = plt.Polygon(verts, facecolor='0.8', edgecolor='k')
plt.gca().add_patch(poly)
                
# scipy.integrateで積分を計算
I = scipy.integrate.quad(f, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分
N = 100000
x = uniform(loc=a, scale=b-a).rvs(size=N)
I = (b - a) * np.mean(f(x))
print u"モンテカルロ積分:", I


    

    




scipy.integrate: 0.96520093605
モンテカルロ積分: 0.966150133367






    

In [24]:

    

# 練習問題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


    

    




scipy.integrate: 3.43506155523
モンテカルロ積分(1): 3.4456557745
モンテカルロ積分(2): 3.43354612482

In [28]:

    

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform
import scipy.integrate

# モンテカルロ積分の収束テスト
# 例3.3の場合

N = 10000

a, b = 0, 1
h = lambda x: (np.cos(50 * x) + np.sin(20 * x)) ** 2

# scipy.integrateで積分を計算
I = scipy.integrate.quad(h, a, b)[0]
print "scipy.integrate:", I

# モンテカルロ積分の収束テスト
x = h(uniform(loc=a, scale=b-a).rvs(size=N))

# サンプル数1のh_1からサンプル数Nのh_Nまで推定値をまとめて求める
estint = np.cumsum(x) / np.arange(1, N + 1)

# サンプル数1のsqrt(v_1)からサンプル数Nのsqrt(v_N)まで標準偏差をまとめて求める
esterr = np.sqrt(np.cumsum((x - estint) ** 2)) / np.arange(1, N + 1)

plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='gray')
plt.plot(estint - 2 * esterr, color='gray')
plt.ylim((0, 2))
plt.show()


    

    




scipy.integrate: 0.96520093605






    

In [67]:

    

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cauchy, norm
import scipy.integrate

# モンテカルロ積分の収束テスト
# 練習問題3.1

N = 1000
x = 4

# scipy.integrateでの積分
h1 = lambda t: t * norm(loc=x).pdf(t) * cauchy.pdf(t)
h2 = lambda t: norm(loc=x).pdf(t) * cauchy.pdf(t)
num = scipy.integrate.quad(h1, -Inf, Inf)[0]
den = scipy.integrate.quad(h2, -Inf, Inf)[0]
I = num / den
print "scipy.integrate:", I

# (1) コーシー分布からサンプリングするモンテカルロ積分の収束テスト
# 分子も分母も同じサンプルを使用すると仮定
theta = cauchy.rvs(size=N)
num = theta * norm(loc=x).pdf(theta)
den = norm(loc=x).pdf(theta)

# 分母に0がくるサンプルを削除
num = num[den != 0]
den = den[den != 0]
Ndash = len(num)

y = num / den
estint = (np.cumsum(num) / np.arange(1, Ndash + 1)) / (np.cumsum(den) / np.arange(1, Ndash + 1))
esterr = np.sqrt(np.cumsum((y - estint) ** 2)) / np.arange(1, Ndash + 1)

plt.subplot(2, 1, 1)
plt.plot(estint, color='red', linewidth=2)
plt.plot(estint + 2 * esterr, color='pink', linewidth=1)
plt.plot(estint - 2 * esterr, color='pink', linewidth=1)
plt.title('sampling from cauchy distribution')
plt.ylim((0, 6))

# (2) 正規分布からサンプリングするモンテカルロ積分の収束テスト
theta = norm(loc=x).rvs(size=N)
num = theta * cauchy.pdf(theta)
den = cauchy.pdf(theta)

num = num[den != 0]
den = den[den != 0]
Ndash = len(num)

y = num / den
estint = (np.cumsum(num) / np.arange(1, Ndash + 1)) / (np.cumsum(den) / np.arange(1, Ndash + 1))
esterr = np.sqrt(np.cumsum((y - estint) ** 2)) / np.arange(1, Ndash + 1)

plt.subplot(2, 1, 2)
plt.plot(estint, color='blue', linewidth=2)
plt.plot(estint + 2 * esterr, color='cyan', linewidth=1)
plt.plot(estint - 2 * esterr, color='cyan', linewidth=1)
plt.title('sampling from normal distribution')
plt.ylim((0, 6))

plt.tight_layout()
plt.show()


    

    




scipy.integrate: 3.43506155523