python monte_carlo_integration2.py

Posted 2021-05-10

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