来自自定义似然函数的样本

Question

我在一个相当复杂的模型中使用了以下似然函数（实际上是对数尺度）：

library(plyr)
dcustom=function(x,sd,L,R){
    R. = (log(R) - log(x))/sd
    L. = (log(L) - log(x))/sd
    ll = pnorm(R.) - pnorm(L.)
    return(ll) 
}

df=data.frame(Range=seq(100,500),sd=rep(0.1,401),L=200,U=400)
df=mutate(df, Likelihood = dcustom(Range, sd,L,U))

with(df,plot(Range,Likelihood,type='l'))
abline(v=200)
abline(v=400)

在这个函数中，sd是预先确定的，L和R是"observations"（很像均匀分布的端点），所以3个都给定了。如果模型估计 x（派生参数）在 L-R 范围之间，则上述函数提供了很大的可能性 (1)，在边界附近（其锐度取决于 sd）平滑的可能性减少（在 0 和 1 之间） , 如果外面太多则为 0。

这个函数非常适合获取 x 的估计值，但现在我想做相反的事情：从上面的函数中随机抽取一个 x。如果我多次这样做，我会生成一个遵循上面绘制的曲线形状的直方图。

最终目标是在 C++ 中执行此操作，但我认为如果我能先弄清楚如何在 R 中执行此操作对我来说会更容易。

网上有一些有用的信息可以帮助我开始 (http://matlabtricks.com/post-44/generate-random-numbers-with-a-given-distribution, https://stats.stackexchange.com/questions/88697/sample-from-a-custom-continuous-distribution-in-r)，但我仍然不完全确定如何做以及如何编码。

我想（完全不确定！）步骤是：

将似然函数转化为概率分布
计算累积分布函数
逆变换采样

这是否正确？如果正确，我该如何编码？谢谢。

Answer 1

一个想法可能是使用 Metropolis Hasting 算法从给定所有其他参数和您的可能性的分布中获取样本。

# metropolis hasting algorithm
 set.seed(2018)
 n_sample <- 100000
 posterior_sample <- rep(NA, n_sample)
 x <- 300 # starting value: I chose 300 based on your likelihood plot
 for (i in 1:n_sample){
     lik <- dcustom(x = x, sd = 0.1, L = 200, R =400)
     # propose a value for x (you can adjust the stepsize with the sd)
     x.proposed <-  x + rnorm(1, 0, sd = 20)
     lik.proposed <- dcustom(x = x.proposed, sd = 0.1, L = 200, R = 400)
     r <- lik.proposed/lik # this is the acceptance ratio
         # accept new value with probablity of ratio 
         if (runif(1) < r) {
         x <- x.proposed
     posterior_sample[i] <- x
     }
    }

 # plotting the density  
 approximate_distr <- na.omit(posterior_sample) 
 d <- density(approximate_distr)
 plot(d, main = "Sample from distribution")
 abline(v=200)
 abline(v=400)

# If you now want to sample just a few values (for example, 5) you could use 
 sample(approximate_distr,5)
#[1] 281.7310 371.2317 378.0504 342.5199 412.3302

来自自定义似然函数的样本

Sample from a custom likelihood function

r

log-likelihood