0

我想将先前的分布信息输入到函数中。我可以通过修改函数体手动输入单个分布，但我正在寻找一种通用的方法？例如，我想绘制给定先验分布的函数的后验分布。

set.seed(1)
n <- 10
pars <- runif(n)
y <- NA
for (i in 1:n)
  y[i] <- rbinom(1,1, prob=pars[i])

plotPosterior <- function(pars,y,mean=0,vari=4)
{  
  x <- seq(-3,3,by = .1)
  logLik <- NA
  for (i in seq(along.with=x))
    logLik[i] <- sum(y*log(1 + exp(pars-x[i])) - (y-1)*log(1 + exp(x[i]-pars)))

  posterior <- logLik * dnorm(x,mean=mean,sd=sqrt(vari))
  plot(x,posterior,type="l")
}
plotPosterior(pars,y,0,4)

我可以输入正态分布的均值方差参数。但是，如果我想使用例如 beta 发行版，我必须重写该函数。相反，我想要一种输入诸如“ dnorm(mean=xx,sd=yy)”或“ dbeta(shape1=xx, shape2=yy)”之类的分布的方法……
我看到的唯一可行的方法是dnorm(x,mean=mean,sd=sqrt(vari))将函数作为输入输入。但我不想预先指定x。有没有其他方法可以做到这一点？

4

1 回答 1

1

为了清楚起见，这是从评论中提取的工作解决方案，

set.seed(1)
n <- 10
pars <- runif(n)
y <- NA
for (i in 1:n)
  y[i] <- rbinom(1,1, prob=pars[i])

plotPosterior <- function(pars,y, fun = dnorm, 
                          params.fun = list(mean=0, sd=2))
{  
  x <- seq(-3,3,by = .1)
  logLik <- NA
  for (i in seq(along.with=x))
    logLik[i] <- sum(y*log(1 + exp(pars-x[i])) - (y-1)*log(1 + exp(x[i]-pars)))

  posterior <- logLik * do.call(fun, c(list(x), params.fun))
  plot(x,posterior,type="l")
}

plotPosterior(pars, y) # default params and function
plotPosterior(pars, y, fun = dbeta, params.fun = list(shape1=2, shape2=3))
于 2013-06-18T21:32:19.237 回答