From 140329dd009e9f9a924c3d8d9421d6942436ee2a Mon Sep 17 00:00:00 2001 From: Pierre Gloaguen Date: Wed, 21 Aug 2024 10:26:01 +0200 Subject: [PATCH 1/4] script fonctionnel --- script_dZInormal.R | 82 +++++++++++++++++++++------------------------- 1 file changed, 37 insertions(+), 45 deletions(-) diff --git a/script_dZInormal.R b/script_dZInormal.R index 85a7c59..01b08ec 100644 --- a/script_dZInormal.R +++ b/script_dZInormal.R @@ -1,17 +1,7 @@ -# my_function <- function(x, prob, mu, sigma){ -# pos_indexes <- which(x != 0) -# zero_indexes <- which(x == 0) -# n_zeros <- length(x) - length(pos_indexes) -# p_term <- sum(log(prob[zero_indexes])) + -# sum(log(1 - prob[pos_indexes])) -# mu_term <- 0 -# if(length(pos_indexes) > 0){ -# mu_term <- mixtools::logdmvnorm(x[pos_indexes], -# mu[pos_indexes], -# sigma[pos_indexes, pos_indexes]) -# } -# p_term + mu_term -# } + + +library(nimble) +library(ggmcmc) dZInormal <- nimbleFunction( run = function(x = double(1), @@ -35,21 +25,14 @@ dZInormal <- nimbleFunction( } log_output <- p_term + mu_term if(log){ - return(exp(log_output)) + return(log_output) } else{ - return(log_output) + return(exp(log_output)) } } ) -# -# sum(apply(Y, 1, -# function(y) -# dZInormal(y, rep(0.4, 5), mu = 0:4, sigma = diag(1, 5)))) -# -# sum(apply(Y, 1, -# function(y) -# dZInormal(y, rep(0.4, 5), mu = rep(0, 5), sigma = diag(1, 5)))) + my_code <- nimbleCode({ for(j in 1:p){ @@ -63,38 +46,47 @@ my_code <- nimbleCode({ }) registerDistributions(list( - dZInormal = list(BUGSdist = "dZInormal(prob, mu, sigma)", - discrete = FALSE, pqAvail = FALSE, - types = c('value = double(1)', 'prob = double(1)', - 'mu = double(1)', 'sigma = double(2)')) + dZInormal = list(BUGSdist = "dZInormal(prob, mu, sigma)", # How to call in nimble + discrete = FALSE, # Distribution is not discrete + pqAvail = FALSE, # CDF and quantile function are not available + types = c('value = double(1)', # The random variable is a vector + 'prob = double(1)', # a vector + 'mu = double(1)', # vector + 'sigma = double(2)')) # double(2) is a matrix )) -n_obs <- 5000 + +# Generating data +set.seed(123) +n_obs <- 1000 n_species <- 3 -my_sigma <- diag(1, n_species) -U <- mixtools::rmvnorm(n = n_obs, mu = 0:(n_species - 1), sigma = my_sigma) +# Values +U <- mixtools::rmvnorm(n = n_obs, + mu = 0:(n_species - 1), + sigma = diag(1, n_species)) +# Mask (matrix of zeros and ones) Z <- rbinom(n = n_obs * n_species, size = 1, prob = .8) %>% matrix(nrow = n_obs) +# Observations Y <- round(U * Z, 9) -my_model <- nimbleModel(code = my_code, - constants = list(p = n_species, - n = n_obs, + + +my_model <- nimbleModel(code = my_code, + constants = list(p = n_species, + n = n_obs, Ip = diag(1, n_species)), data = list(Y = Y), - inits = list(mu = 0:n_species, + inits = list(mu = rep(0, n_species), prob = rep(0.5, n_species), sigma = diag(1, n_species))) - results <- nimbleMCMC(my_model, - nchains = 1, niter = 10000, nburnin = 1000, thin = 10) -hist(results[,"prob[]"]) - + samplesAsCodaMCMC = TRUE, + nchains = 2, niter = 10000, + nburnin = 1000, thin = 10) -y <- c(0, 0, 2, 1) +ggs(results) %>% + ggplot(aes(x = Iteration, y = value, color = factor(Chain))) + + facet_wrap(~Parameter, scales = "free") + + geom_line() -ps <- runif(length(y)) -mus <- rnorm(length(y)) -sigma <- rWishart(n = 1, Sigma = diag(1, length(y)), df = length(y))[,,1] -my_function(y, ps, mus, sigma) -dZInormal(y, ps, mus, sigma) \ No newline at end of file From f2552b9a1e4514b81e07095a734c1f3091ccc4e9 Mon Sep 17 00:00:00 2001 From: Pierre Gloaguen Date: Wed, 21 Aug 2024 11:31:05 +0200 Subject: [PATCH 2/4] virer dependence en mixtools --- 01_nimble.qmd | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/01_nimble.qmd b/01_nimble.qmd index 39cfcb1..da5c68e 100644 --- a/01_nimble.qmd +++ b/01_nimble.qmd @@ -703,13 +703,15 @@ And now, let's run it! ```{r} +#| label: fit-ZI-normal +#| cache: true # Generating data set.seed(123) n_obs <- 1000 n_species <- 3 # Values -U <- mixtools::rmvnorm(n = n_obs, - mu = 0:(n_species - 1), +U <- mvtnorm::rmvnorm(n = n_obs, + mean = 0:(n_species - 1), sigma = diag(1, n_species)) # Mask (matrix of zeros and ones) Z <- rbinom(n = n_obs * n_species, size = 1, prob = .8) %>% From 5f4178bce877ed0e5d5c985918ac45f021b6709b Mon Sep 17 00:00:00 2001 From: Pierre Gloaguen Date: Wed, 21 Aug 2024 11:32:38 +0200 Subject: [PATCH 3/4] renommage pour ignorer --- script_dZInormal.R | 92 ---------------------------------------------- 1 file changed, 92 deletions(-) delete mode 100644 script_dZInormal.R diff --git a/script_dZInormal.R b/script_dZInormal.R deleted file mode 100644 index 01b08ec..0000000 --- a/script_dZInormal.R +++ /dev/null @@ -1,92 +0,0 @@ - - -library(nimble) -library(ggmcmc) - -dZInormal <- nimbleFunction( - run = function(x = double(1), - prob = double(1), - mu = double(1), - sigma = double(2), - log = integer(0, default = 0)){ - returnType(double(0)) - non_nul_indexes <- which(x!=0) - nul_indexes <- which(x == 0) - p_term <- sum(log(prob[nul_indexes])) + - sum(log(1 - prob[non_nul_indexes])) - mu_term <- 0 - if(length(non_nul_indexes) > 0){ - chol_mat <- chol(sigma[non_nul_indexes, non_nul_indexes]) - restricted_x = x[non_nul_indexes] - restricted_mu = mu[non_nul_indexes] - mu_term <- dmnorm_chol(restricted_x, - restricted_mu, - chol_mat, prec_param = FALSE, log = TRUE) - } - log_output <- p_term + mu_term - if(log){ - return(log_output) - } - else{ - return(exp(log_output)) - } - } -) - - -my_code <- nimbleCode({ - for(j in 1:p){ - prob[j] ~ dunif(0, 1) - mu[j] ~ dnorm(0, 1) - } - sigma[1:p, 1:p] ~ dwish(Ip[1:p, 1:p], p) - for(i in 1:n){ - Y[i, 1:p] ~ dZInormal(prob[1:p], mu[1:p], sigma[1:p, 1:p]) - } -}) - -registerDistributions(list( - dZInormal = list(BUGSdist = "dZInormal(prob, mu, sigma)", # How to call in nimble - discrete = FALSE, # Distribution is not discrete - pqAvail = FALSE, # CDF and quantile function are not available - types = c('value = double(1)', # The random variable is a vector - 'prob = double(1)', # a vector - 'mu = double(1)', # vector - 'sigma = double(2)')) # double(2) is a matrix -)) - - -# Generating data -set.seed(123) -n_obs <- 1000 -n_species <- 3 -# Values -U <- mixtools::rmvnorm(n = n_obs, - mu = 0:(n_species - 1), - sigma = diag(1, n_species)) -# Mask (matrix of zeros and ones) -Z <- rbinom(n = n_obs * n_species, size = 1, prob = .8) %>% - matrix(nrow = n_obs) -# Observations -Y <- round(U * Z, 9) - - - -my_model <- nimbleModel(code = my_code, - constants = list(p = n_species, - n = n_obs, - Ip = diag(1, n_species)), - data = list(Y = Y), - inits = list(mu = rep(0, n_species), - prob = rep(0.5, n_species), - sigma = diag(1, n_species))) -results <- nimbleMCMC(my_model, - samplesAsCodaMCMC = TRUE, - nchains = 2, niter = 10000, - nburnin = 1000, thin = 10) - -ggs(results) %>% - ggplot(aes(x = Iteration, y = value, color = factor(Chain))) + - facet_wrap(~Parameter, scales = "free") + - geom_line() - From 0d6ab3f15e9f4e2421196f92691a35539bb3f274 Mon Sep 17 00:00:00 2001 From: Marie-Pierre Etienne Date: Wed, 21 Aug 2024 10:10:05 +0200 Subject: [PATCH 4/4] Update pr_action.yml triggered at every changes --- .github/workflows/pr_actions.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/pr_actions.yml b/.github/workflows/pr_actions.yml index 069c7c9..3970a31 100644 --- a/.github/workflows/pr_actions.yml +++ b/.github/workflows/pr_actions.yml @@ -3,7 +3,7 @@ name: pr_check on: # Triggers the workflow on push or pull request events but only for the master branch pull_request: - types: [opened, reopened, edited] + branches: [main, master] # A workflow run is made up of one or more jobs that can run sequentially or in parallel