如何使用其概率函数最好地模拟任意单变量随机变量？

Question

在R中,如果只有概率密度函数可用,那么模拟任意单变量随机变量的最佳方法是什么？

Answer 1

当您仅获得密度时,这是逆cdf方法的(慢)实现.

den<-dnorm #replace with your own density

#calculates the cdf by numerical integration
cdf<-function(x) integrate(den,-Inf,x)[[1]]

#inverts the cdf
inverse.cdf<-function(x,cdf,starting.value=0){
 lower.found<-FALSE
 lower<-starting.value
 while(!lower.found){
  if(cdf(lower)>=(x-.000001))
   lower<-lower-(lower-starting.value)^2-1
  else
   lower.found<-TRUE
 }
 upper.found<-FALSE
 upper<-starting.value
 while(!upper.found){
  if(cdf(upper)<=(x+.000001))
   upper<-upper+(upper-starting.value)^2+1
  else
   upper.found<-TRUE
 }
 uniroot(function(y) cdf(y)-x,c(lower,upper))$root
}

#generates 1000 random variables of distribution 'den'
vars<-apply(matrix(runif(1000)),1,function(x) inverse.cdf(x,cdf))
hist(vars)

Answer 2

澄清上面"使用Metropolis-Hastings"的答案:

假设ddist()是你的概率密度函数

就像是:

n <- 10000
cand.sd <- 0.1
init <- 0
vals <- numeric(n)
vals[1] <- init 
oldprob <- 0
for (i in 2:n) {
    newval <- rnorm(1,mean=vals[i-1],sd=cand.sd)
    newprob <- ddist(newval)
    if (runif(1)<newprob/oldprob) {
        vals[i] <- newval
    } else vals[i] <- vals[i-1]
   oldprob <- newprob
}

笔记:

1. 完全未经测试
2. 效率取决于候选人分布(即价值cand.sd).为了获得最大效率,调整cand.sd到25-40%的接受率
3. 结果将是自相关的...(虽然我猜你总是可以 sample()把结果加到他们身上,或者瘦身)
4. 如果您的起始值很奇怪,可能需要丢弃"老化"

这个问题的经典方法是拒绝抽样(参见例如Press et al Numerical Recipes)