如何解决从后验预测分布中采样stan vs inla
我正在尝试在bayesplot
对象上通过INLA
包实现函数,并且有点不确定如何从后验预测分布中得出。我想我差不多了,但是rstan
的抽奖比INLA
的抽奖更具可变性。
在rstan
中,使用bayesplot
vignette的简化示例,我可以:
library(bayesplot)
library(ggplot2)
library(rstanarm)
library(ggpubr)
library(tidyverse)
#rstan model set up
roaches$roach100 <- roaches$roach1 / 100 # pre-treatment number of roaches (in 100s)
fit_poisson <- stan_glm(y ~ roach100 + treatment + senior,offset = log(exposure2),family = poisson(link = "log"),data = roaches,seed = 1111,refresh = 0)
#In order to use the PPC functions from the bayesplot package we need a vector y of outcome values:
y <- roaches$y
#and a matrix yrep of draws from the posterior predictive distribution,yrep_poisson <- posterior_predict(fit_poisson,draws = 500)
#then plot:
p1 <- bayesplot::ppc_dens_overlay(y,yrep_poisson[1:50,])
p1
我想将该图复制到INLA
对象上。根据{{1}}小插图,您可以执行此操作,因为它们提供了定义简单的bayesplot
方法的代码,该方法创建了例如类的拟合模型对象。 pp_check
:
foo
要使用pp_check.foo <- function(object,type = c("multiple","overlaid"),...) {
type <- match.arg(type)
y <- object[["y"]]
yrep <- object[["yrep"]]
stopifnot(nrow(yrep) >= 50)
samp <- sample(nrow(yrep),size = ifelse(type == "overlaid",50,5))
yrep <- yrep[samp,]
if (type == "overlaid") {
ppc_dens_overlay(y,yrep,...)
} else {
ppc_hist(y,...)
}
}
,我们只需列出包含pp_check.foo
和y
组件的列表,并将其赋予foo类:
yrep
INLA
x <- list(y = rnorm(200),yrep = matrix(rnorm(1e5),nrow = 500,ncol = 200))
class(x) <- "foo"
#create plot above:
pp_check(x,type = "overlaid")
#create same model but in inla:
library(INLA)
fit_poisson_inla <- inla(y ~ roach100 + treatment + senior,control.predictor = list(compute = T),family = "poisson")
返回每个inla_object_name$marginals.fitted.values
的后验预测分布:
y
我认为从中反复采样将满足我的需要,但是实际上创建每个分布时,每个观察值只能使用75个值(fit_poisson_inla$marginals.fitted.values
#so to get distribution for first oberservation:
fitted.Predictor.1 <- fit_poisson_inla$marginals.fitted.values[[1]]
),而实际上我希望从所有范围的值中采样。我认为我们可以通过使用线性预测变量使用dim(fitted.Predictor.1)
来做到这一点(第4.3 here节):
inla.tmarginal
我的问题是,如何从这个fitted_dist <- fit_poisson_inla$marginals.linear.predictor
#should i have used "inla.rmarginal(n,marginal)"?
marginal_dist <- lapply(fitted_dist,function(y) inla.tmarginal(function(x) {exp(x)},y)) %>% map(~ as.data.frame(.) %>% rename(.,xx = x))
#resample 500 times
yrep_poisson_inla <- as.matrix(bind_rows(rerun(500,lapply(marginal_dist,function(x) sample(x$xx,1)) %>% as.data.frame())))
#convert to class foo for pp_check
x <- list(y = y,yrep = yrep_poisson_inla[1:50,])
class(x) <- "foo"
p2 <- pp_check(x,type = "overlaid")
#plot
ggarrange(p1,p2,ncol = 1,nrow = 2,labels = c("rstan","inla sample"))
(inla
)对象的后验预测分布中正确得出平局矩阵并传递到fit_poisson_inla
? pp_check
产生离散值,而yrep_poisson
产生连续值。 yrep_poisson_inla
抽奖中的变化远远大于rstan
(第二个图)。我做的是否正确,这只是一些抽样问题,还是不同方法的产物?在更复杂的示例中,差异可能很大。
谢谢
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。