---
title: "Using priorsense with NIMBLE"
vignette: >
  %\VignetteIndexEntry{Using priorsense with NIMBLE}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---


```{r}
#| include: false
ggplot2::theme_set(bayesplot::theme_default(base_family = "sans"))

options(priorsense.plot_help_text = FALSE)
```


```{r}
#| message: false
#| warning: false
library(nimble)
library(priorsense)
```


`priorsense` is compatible with models fit with `nimble`

Consider the univariate normal model with unknown mu and sigma available
via`example_powerscale_model("univariate_normal")`. In NIMBLE, `lprior` and
`log_lik` variables can be defined as below. By also defining separate
`lprior_mu` and `lprior_sigma` variables, it will be possible to check the
sensitivity for each prior separately.

```{r}
model <- example_powerscale_model(language = "nimble")
```

```{r}
#| echo: false
#| results: asis
cat("```\n")
model$model_code
cat("```")
```

Instantiate and compile the model as usual.
```{r}
#| message: false
#| warning: false
#| results: false

inits <- list(
  mu = 0,
  sigma = 1
)

model <- nimbleModel(
  model$model_code,
  data = model$data,
  inits = inits,
  constants = list(N = model$data$N)
)


cmodel <- compileNimble(model)

mcmc <- buildMCMC(
  model,
  monitors = c(
    "mu",
    "sigma",
    "lprior",
    "lprior_mu",
    "lprior_sigma",
    "log_lik"
  )
)

cmcmc <- compileNimble(mcmc, project = model)
```

For `priorsense` compatibility, it is easiest to save the draws as a
`coda::mcmc.list` object, by specifying `samplesAsCodaMCMC = TRUE`. Otherwise,
the resulting object can be transformed to a `posterior::draws_df` object, with
`posterior::as_draws_df`.

```{r}
#| message: false
#| warning: false
#| results: false

fit <- runMCMC(
  cmcmc,
  niter = 20000,
  nburnin = 5000,
  nchains = 4,
  thin = 5,
  setSeed = c(123, 456, 789, 101112),
  samplesAsCodaMCMC = TRUE # alternatively, coerce the output using `posterior::as_draws_df`
)
```

Then the `priorsense` functions will work as usual.

```{r}
powerscale_sensitivity(fit)
```

```{r}
powerscale_sensitivity(fit, prior_selection = "sigma")
```

```{r}
powerscale_sensitivity(fit, prior_selection = "mu")
```


```{r}
#| message: false
#| warning: false
#| fig-width: 6
#| fig-height: 4
powerscale_plot_dens(fit)
```