estim_jags generates inputs for JAGS (data list and model string or external file)
for sampling from the posterior distribution,
given samples of the bivariate chi-bar-squared distribution.
estim_jags(samples, d, dimC = d, linC = 0, v_prior = NA, prior_sample_size = 1, filename = NA, overwrite = FALSE)
| samples | N-by-2 matrix representing independent samples from the bivariate chi-bar-squared distribution of a convex cone |
|---|---|
| d | the dimension of the ambient space |
| dimC | the dimension of the cone |
| linC | the lineality of the cone |
| v_prior | a prior estimate of the vector of intrinsic volumes (NA by default) |
| prior_sample_size | the sample size for the prior estimate (1 by default -> noninformative) |
| filename | filename for output (NA by default, in which case the return is a string) |
| overwrite | logical; determines whether the output should overwrite an existing file |
If filename==NA then the output of estim_jags is a list containing the following elements:
model: a string that forms the description of the JAGS model,
which can be directly used as input (via text connection
or external file) for creating a JAGS model object;
post_samp(n) returns an n-by-(dimC+1)
matrix whose rows form a set of n independent samples
of the posterior distribution,
data: a data list containing the prepared data to be used
with the model string for defining a JAGS model object,
variable.names: the single string "V" to be used as additional
parameter when creating samples from the JAGS object to
indicate that only this vector should be tracked.
If filename!=NA then the model string will be written to the file with
the specified name and the output will only contain the elements data
and variable.names.
For the JAGS models enforcing log-concavity is not supported; use the
analogous Stan model instead, provided by estim_stan.
See the example below for details on how to use the outputs with rjags functions; see this vignette for further info.
Package: conivol
# NOT RUN { library(tidyverse) library(rjags) options(mc.cores = parallel::detectCores()) # defining the example cone D <- c(5,8) alpha <- c( asin(sqrt(0.9)) , asin(sqrt(0.8))) v_exact <- circ_ivols(D, alpha, product = TRUE) d <- sum(D) # getting the sample data N <- 10^3 set.seed(1234) m_samp <- rbichibarsq(N,v_exact) # compute initial guess est <- estim_statdim_var(d, m_samp) v0 <- init_ivols(d,init_mode=1,delta=est$delta,var=est$var) # obtain input data for JAGS model; use v0 as prior in_jags <- estim_jags(m_samp, d, v_prior=v0) # create JAGS model model_connection <- textConnection(in_jags$model) mod <- jags.model(model_connection , data = in_jags$data , n.chains = 4 , n.adapt = 500) close(model_connection) # update(mod, 1e3) # simulate posterior distribution and display trace plots and summaries mod_sim <- coda.samples(model=mod, variable.names=in_jags$variable.names, n.iter=1e4) # plot(mod_sim, ask=TRUE) mod_csim <- as.mcmc(do.call(rbind, mod_sim)) tib_plot <- tibble( x=0:d, y_true=v_exact, y_est=v0, y_bayes_mean=summary(mod_csim)$statistics[ ,'Mean'], y_bayes_quant25=summary(mod_csim)$quantiles[ ,'25%'], y_bayes_quant50=summary(mod_csim)$quantiles[ ,'50%'], y_bayes_quant75=summary(mod_csim)$quantiles[ ,'75%'] ) # plot true values, the estimate v0, and marginals of the posterior samples ggplot(data=tib_plot) + geom_line(aes(x=x, y=y_true), color="red") + geom_line(aes(x=x, y=y_est), color="black") + geom_line(aes(x=x, y=y_bayes_mean), color="blue") + geom_line(aes(x=x, y=y_bayes_quant25), color="green") + geom_line(aes(x=x, y=y_bayes_quant50), color="green") + geom_line(aes(x=x, y=y_bayes_quant75), color="green") + theme(axis.title.x=element_blank(), axis.title.y=element_blank()) # }