# -*- coding: utf-8 -*-
"""
重要性采樣
題目: 計算y = f(x) = x**2 + 4*x*np.sin(x) 在【1，5】([a,b])之間的定積分
均勻分布: g(x) = 1/(b-a)
令 f(x) = (f(x)/g(x))*g(x)
則∫f(x)dx = ∫(f(x)/g(x))*g(x)dx = (1/N)∑(f(x)/g(x))
"""

import numpy as np
import matplotlib.pyplot as plt


def f(x):  # 求f(x)在2-3之間的定積分
    return (x ** 2 + 4 * x * np.sin(x))


def intf(x):  # 原函數
    return x ** 3 / 3.0 + 4.0 * np.sin(x) - 4.0 * x * np.cos(x)


def g(x):
    return 1 / 4


def z(x):  # f(x)/g(x)
    return f(x) / g(x)


a = 1
b = 5
N = 10000  # 一萬次實驗
X = np.random.uniform(low=a, high=b, size=N)  # 利用均勻分布隨機生成N個區間為【a,b】的數
Y = z(X)  # 帶入函數f(x)計算對應y值

Imc = np.sum(Y) / N  # 計算N次結果的均值
print("Monte Carlo estimation=", Imc)

exactVal = intf(b) - intf(a)  # 精確值，用原函數求解
print("Exactly number=", exactVal)

# 繪圖
error = abs((np.cumsum(Y) / np.arange(1, 10001)) - exactVal)  # 誤差絕對值
plt.plot(np.arange(1, 10001), error, alpha=0.7)
plt.xlabel("N")
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.ylabel("誤差絕對值")
plt.show()