重要性采样

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

def f(x):
    return np.sin(x) * x;

def intf(x1, x2):
    return (np.sin(x2)- x2 * np.cos(x2)) - (np.sin(x1) - x1 * np.cos(x1));

def p(x, mu, sig):
    return (1 / np.sqrt(2 * np.pi * sig ** 2)) * np.exp(-(x - mu) ** 2 /(2.0 * sig ** 2));

def intp(x1, x2, mu, sig):
    return (norm.cdf(x2 - mu, scale = sig) - norm.cdf(x1 - mu, scale = sig));

def plot(xmin, xmax, mu, sig):
    x = np.linspace(xmin, xmax, 1000);
    plt.subplot(1,2,1);
    plt.plot(x, f(x), 'b', label = u'Original $x\sin(x)$');
    plt.plot(x, p(x, mu, sig), 'r', label = u'Importance Sampling Function: Normal');
    plt.xlabel('x');
    plt.legend();
    plt.show();

def demo():
    # Example: Calculate ∫sin(x)xdx
    mu = 2;
    sig =.7;
    xmax = np.pi 
    xmin = 0
    N = 1000

    plot(xmin, xmax, mu, sig);
    answer = intf(xmin, xmax);
    
    x = np.random.uniform(low = xmin, high = xmax, size = N)
    monteCarlo = (xmax - xmin) * np.mean(f(x))

    xis = mu + sig * np.random.randn(N, 1);
    xis = xis[(xis < xmax) & (xis > xmin)];#概率密度函数在采样区间[0 pi]上的积分需要等于1
    normal = intp(xmin, xmax, mu, sig)
    importanceSampling =np.mean(f(xis) / p(xis, mu, sig)) * normal  # 因此,此处需要除一个系数即p(x)在[0 pi]上的积分
    
    print(monteCarlo);
    print(importanceSampling);
    print(answer);

if __name__ == '__main__':
    demo();

先上段代码，占个坑。过几天补上文字说明；

结合新章节深入；

https://www.jb51.net/article/130511.htm