# -*- 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()