2

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

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 是“观察值”（非常类似于均匀分布的端点），因此给出了所有 3 个。如果模型估计 x（派生参数）在 LR 范围之间，则上述函数提供了很大的似然性（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）但我仍然不完全确定如何去做以及如何编码。

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

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

这是正确的，如果是这样，我该如何编码？谢谢你。

4

1 回答 1

2

一个想法可能是使用 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
于 2018-02-26T20:52:50.633 回答