diff --git a/.gitignore b/.gitignore index 31e7331..dd718b2 100644 --- a/.gitignore +++ b/.gitignore @@ -52,5 +52,7 @@ rsconnect/ # For niimble branch _site/ +*_files/ +*_cache/ 01_nimble.html local_nimble*.R diff --git a/01_nimble.qmd b/01_nimble.qmd index b85f18f..9f3e74f 100644 --- a/01_nimble.qmd +++ b/01_nimble.qmd @@ -40,13 +40,13 @@ But, this remains limited and to have full flexibility, ones need to build its o ```{r} #| label: example1_data - -data_ex1 <- rnbinom(n = 10, prob = 0.4, size = 12) +data_ex1 <- rnbinom(n = 1000, prob = 0.4, size = 12) ``` ```{r} +#| label: code_neg_bin code_neg_bin <- nimbleCode({ # Observation model for(i in 1:n){# n is never defined before, it will be a constant @@ -77,11 +77,125 @@ model_neg_bin <- nimbleModel(code = code_neg_bin, ```{r} #| label: posterior_samples #| cache: true -posterior_samples_neg_bin <- nimbleMCMC(model_neg_bin) +posterior_samples_neg_bin <- nimbleMCMC(model_neg_bin, + nchains = 3, + niter = 10000, + thin = 10, + nburnin = 1000) ``` ## Exploring the results +```{r} +#| label: str_posterior_samples +str(posterior_samples_neg_bin) +``` + +```{r} +#| label: custom_formatted_results +formatted_results <- imap_dfr(posterior_samples_neg_bin, + function(x, nm){ + as.data.frame(x) %>% + rowid_to_column(var = "Iteration") %>% + mutate(Chain = str_remove(nm, "chain")) + }) %>% + pivot_longer(cols = -c("Iteration", "Chain"), + names_to = "Parameter", + values_to = "value") +``` + +```{r} +#| label: plot +ggplot(formatted_results) + + aes(x = Iteration, + y = value, color = Chain) + + facet_wrap(~Parameter, scales = "free") + + geom_line() + + labs(x = "Sample ID", y = "Parameter value", color = "") +``` + +## Package for automatic formatting + +```{r} +#| label: posterior_samples_coda +#| cache: true +posterior_samples_neg_bin <- nimbleMCMC(model_neg_bin, + nchains = 3, + niter = 10000, + thin = 10, + nburnin = 1000, + samplesAsCodaMCMC = TRUE) +``` + +```{r} +#| label: str_posterior_samples_coda +str(posterior_samples_neg_bin) +``` + +```{r} +#| label: coda_formatted_results +formatted_results <- ggs(posterior_samples_neg_bin) +formatted_results +``` + +## Defining a `nimbleFunction` + +```{r} +get_prob <- nimbleFunction( + run = function(log_mu = double(0), + log_theta = double(0)) { # type declarations + z <- log_theta - log_mu + output <- 1 / (1 + exp(-z)) + return(output) + returnType(double(0)) # return type declaration + }) +get_prob(log(18), log(12)) +``` + + +```{r} +#| label: code_neg_bin_alternatif +#| cache: true +code_alternatif <- nimbleCode({ + # Observation model + for(i in 1:n){# n is never defined before, it will be a constant + y[i] ~ dnbinom(prob, theta) + } + # Alternative vectorized formulation + # y[1:n] ~ dnbinom(prob, theta) + # PRIORS + log_mu ~ dnorm(0, 1) + theta ~ dexp(0.1) + # Quantites deterministes + log_theta <- log(theta) + prob <- get_prob(log_mu = log_mu, log_theta = log_theta) +}) + +model_alternatif <- nimbleModel(code = code_alternatif, + name = "Alternative negative binomial", + constants = list(n = length(data_ex1)), + data = list(y = data_ex1), + inits = list(mu = 0.5, theta = 1)) + +posterior_samples_alternatif <- nimbleMCMC(model_alternatif, + nchains = 3, + niter = 10000, + thin = 10, + nburnin = 1000, + monitors = c("prob", "theta"), + samplesAsCodaMCMC = TRUE) +``` + +```{r} +posterior_samples_alternatif %>% + ggs() %>% + ggplot() + + aes(x = Iteration, + y = value, color = factor(Chain)) + + facet_wrap(~Parameter, scales = "free") + + geom_line() + + labs(x = "Sample ID", y = "Parameter value", color = "") +``` ## Alternative MCMC sampler