From d9028a066a6ed9d1f537cf1088908c012cef1e6c Mon Sep 17 00:00:00 2001 From: Pierre Gloaguen Date: Mon, 19 Aug 2024 18:23:54 +0200 Subject: [PATCH 1/3] fin exemple basique --- 01_nimble.qmd | 57 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/01_nimble.qmd b/01_nimble.qmd index 42f70d8..cc84cb5 100644 --- a/01_nimble.qmd +++ b/01_nimble.qmd @@ -76,11 +76,66 @@ 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 = sample_id, + y = sample_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 +#| 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 +``` ## Alternative MCMC sampler From 334de99a761c039738daaa80b18f01be302e991e Mon Sep 17 00:00:00 2001 From: Pierre Gloaguen Date: Mon, 19 Aug 2024 18:58:02 +0200 Subject: [PATCH 2/3] ignore _files et _cache --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) 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 From cf9514a661eb6f3a696b658c2394e6da3faf9f5a Mon Sep 17 00:00:00 2001 From: Pierre Gloaguen Date: Mon, 19 Aug 2024 18:58:20 +0200 Subject: [PATCH 3/3] fin exemple avec fonction --- 01_nimble.qmd | 72 ++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 65 insertions(+), 7 deletions(-) diff --git a/01_nimble.qmd b/01_nimble.qmd index cc84cb5..e90fb9f 100644 --- a/01_nimble.qmd +++ b/01_nimble.qmd @@ -39,13 +39,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 @@ -100,15 +100,15 @@ formatted_results <- imap_dfr(posterior_samples_neg_bin, }) %>% pivot_longer(cols = -c("Iteration", "Chain"), names_to = "Parameter", - values_to = "Value") + values_to = "value") ``` ```{r} #| label: plot ggplot(formatted_results) + - aes(x = sample_id, - y = sample_value, color = chain) + - facet_wrap(~parameter, scales = "free") + + aes(x = Iteration, + y = value, color = Chain) + + facet_wrap(~Parameter, scales = "free") + geom_line() + labs(x = "Sample ID", y = "Parameter value", color = "") ``` @@ -116,7 +116,7 @@ ggplot(formatted_results) + ## Package for automatic formatting ```{r} -#| label: posterior_samples +#| label: posterior_samples_coda #| cache: true posterior_samples_neg_bin <- nimbleMCMC(model_neg_bin, nchains = 3, @@ -140,6 +140,64 @@ formatted_results ## Alternative MCMC sampler +## 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 = "") +```