使用Python進行離散拒絕採樣Monte Carlo

我正在使用python來使用拒絕接受方法來採樣離散MC分佈。由於曲線類似於power law，我決定在它周圍設置一個簡單的信封（在x = 77處）以使代碼更快。代碼不按預期執行，不過，因爲它示出了用於與在整個區域上的簡單的矩形圖中，相比於包絡：here 使用Python進行離散拒絕採樣Monte Carlo

x是data.rank，0-5000 之間y是數據.freq

任何人都可以發現代碼有問題嗎？兩個直方圖的輸出應該相等。謝謝！

#!/usr/bin/env python       

import numpy as np 
import matplotlib.pyplot as plt 
import random 
import pandas as pd 

# Read data 
data = pd.read_csv('data.csv') 

# Rejection-sampling MC (with envelope function gx at rank=77) 
N = 1000 
M = 1.0001 
cutoff = 77 
gx = np.ones(cutoff) * max(data['freq']) 
gx = np.append(gx, np.ones(len(data['rank'])-cutoff) * data['freq'][cutoff-1]) 
histx = [] 
while N > 0: 
    rx = random.randint(0,len(gx)-1) 
    ry = random.uniform(0,1) 
    if ry < data['freq'][rx]/(M*gx[rx]): 
     histx.append(rx) 
     N += -1 
plt.hist(histx, bins=100, histtype='stepfilled', color='b',label='Enveloped (Fast)') 

# Rejection-sampling MC (with envelope function gx at rank=77) 
N2 = 1000 
histy = [] 
while N2 > 0: 
    rx2 = random.randint(0,len(gx)-1) 
    ry2 = random.uniform(0,max(data['freq'])) 
    if ry2 < data['freq'][rx2]: 
     histy.append(rx2) 
     N2 += -1 
plt.hist(histy, bins=100, histtype='stepfilled', color='r', alpha=0.5, label='Normal (Slow)') 
plt.legend() 
plt.show()

來源

2015-10-14 GermanoMosconi

這些都是不同的採樣

在第二種情況下的，這樣就

ry2 = random.uniform(0.0, 1.0) 
if ry2 < data['freq'][rx2]/max(data['freq']): 
    accept

在第一種情況下，你有

ry = random.uniform(0.0, 1.0) 
if ry < data['freq'][rx]/gx[rx]: 
    accept

哪裏gx[rx]不等於max(data['freq']) ，你會有差異

一些建議：儘量使用U（0,1）rng呼叫，更容易發現，替換或更改rng。其次，在MC性能（以及正確性之後）中最重要的是，嘗試在主環路外部計算max(data['freq'])或len(gx)之類的東西

來源

2015-10-14 18:01:15

謝謝您的回答和建議。我仍然對此感到困惑。使用「包絡」應該是一種可靠的方法來對相同的曲線進行採樣，如下所示：https：//en.wikipedia.org/wiki/Ziggurat_algorithm。這：http://imgur.com/PNx0XzM是我想要實現的。我也在考慮在日誌空間中計算MC來加快速度，但我也被困在了這個問題上。這是一個分散的問題嗎？ – GermanoMosconi

不，離散不是問題，請參閱更新建議的解決方案 –

將在幾個小時內添加更多，請諒解延遲 –

使用Python進行離散拒絕採樣Monte Carlo

回答

相關問題