2016-12-05 51 views
1

給定一個拉普拉斯分佈提案:使用重要性蒙特卡洛積分抽檢給出的提議功能

g(x) = 1/2*e^(-|x|) 

和樣本大小n = 1000,我想進行蒙特卡洛(MC)的集成用於估計θ:

enter image description here

通過重要性採樣。最終,我想在R達到那個時候,計算R中MC估計值的均值和標準偏差。


編輯(下面的答案後,遲到了)

這就是我對我的R代碼裏面至今:

library(VGAM) 
n = 1000 
x = rexp(n,0.5) 
hx = mean(2*exp(-sqrt(x))*(sin(x))^2) 
gx = rlaplace(n, location = 0, scale = 1) 
+0

我使用了rlaplace函數的VGAM包。 – Chris95

+0

非常感謝您對後期編輯感到非常抱歉! – Chris95

回答

2

enter image description here

現在我們可以寫一個簡單的R函數從拉普拉斯分佈中抽樣:

## `n` is sample size 
rlaplace <- function (n) { 
    u <- runif(n, 0, 1) 
    ifelse(u < 0.5, log(2 * u), -log(2* (1 - u))) 
    } 

另外寫一個函數爲拉普拉斯分佈的密度:

g <- function (x) ifelse(x < 0, 0.5 * exp(x), 0.5 * exp(-x)) 

現在,您的積是:

f <- function (x) { 
    ifelse(x > 0, exp(-sqrt(x) - 0.5 * x) * sin(x)^2, 0) 
    } 

現在我們估計整體使用1000個樣本(set.seed可重複性):

set.seed(0) 
x <- rlaplace(1000) 
mean(f(x)/g(x)) 
# [1] 0.2648853 

還與使用quad的數值積分rature:

integrate(f, lower = 0, upper = Inf) 
# 0.2617744 with absolute error < 1.6e-05