蒙特卡洛采樣、重要性采樣

有一個概率密度函數p(x)，求解隨機變量x基於此概率下某個函數f(x)的期望，表示如下：

如果概率分布形式比較簡單的話，我們可以采用解析的方法：

如果f(x)過於復雜的話，直接求解就非常復雜，我們采用蒙特卡洛的方法。根據大數定理，當采樣數量足夠大的話，采樣樣本可以無限近似地表示原分布，我們可以得到：

 我們令待采樣的分布為p(x)，另一個簡單可采樣且定義域與p(x)相同的概率密度函數為：

我們只需要簡單地從分布中采樣到x，然后分別計算樣本在兩個分布中的概率和函數值。

　　下面我們進行試驗，我們假設要采樣的數據來自均值為1，標准差為1的高斯分布，希望用另一個高斯分布來近似這個分布。我們選取三個：均值為1，標准差分別為1、0.5、2的高斯分布來進行比較。為了更好地看出重要性采樣的效果，這里的函數將選擇一個比較簡單的形式f(x)=x：

import numpy as np
import math
import matplotlib.pyplot as plt

def gaussian(x,u,sigma):    
    """
    param x:要計算概率密度值的點
    param u:均值
    param sigma:方差
    return x的概率密度值
    """
    return math.exp(-(x-u)**2/(2*sigma*sigma))/math.sqrt(2*math.pi*sigma*sigma)

def importance_sampling_test(ori_sigma,sample_sigma):
    """
    param ori_sigma:原始分布p(x)的方差
    param sample_sigma:采樣分布p~(x)的方差
    return
    
    """
    origin = []
    for n in range(10):
        #進行10次計算
        Sum = 0
        for i in range(100000):
            a = np.random.normal(1.0,ori_sigma)
            Sum += a              
        origin.append(Sum)
        
    isample = []
    for n in range(10):
        Sum2 = 0
        for i in range(100000):
            a = np.random.normal(1.0,sample_sigma) #計算從正太分布采樣出來的x
            ua = gaussian(a,1.0,sample_sigma)    #計算采樣概率密度
            na = gaussian(a,1.0,ori_sigma)       #計算原始概率密度
            Sum2 += a*na/ua    
        isample.append(Sum2)
    
    origin = np.array(origin)
    isample = np.array(isample)
    
    print(np.mean(origin),np.std(origin))
    print(np.mean(isample),np.std(isample))
    

importance_sampling_test(1.0,1.0)
importance_sampling_test(1.0,0.5)
importance_sampling_test(1.0,2.0)

xs = np.linspace(-5,6,301)
y1 = [gaussian(x,1.0,1.0) for x in xs]
y2 = [gaussian(x,1.0,0.5) for x in xs]
y3 = [gaussian(x,1.0,2.0) for x in xs]

fig = plt.figure(figsize=(8,5))

plt.plot(xs,y1,label="sigma=1.0")
plt.plot(xs,y2,label="sigma=0.5")
plt.plot(xs,y3,label="sigma=2.0")
plt.legend()
plt.show()