Open giannilmbd opened 2 years ago
njtierney
commented
2 years ago Hi there!
Thanks for posting here, I was wondering if you can make the following file accessible so I can reproduce it?
source("~/Dropbox/Econometrics/R/Functions/LD_distribution/ld_distro.r")
One thing that can be really helpful with checking if it is a standalone reproducible example is to run the code using reprex: https://github.com/tidyverse/reprex
So you could write:
reprex({
# <your code goes here>
})
giannilmbd
commented
2 years ago Hi Nicholas, sure, sorry, I forgot to add that function. Here it is (it is part of the greta example I found on the web)
library (R6) distrib <- .internals$nodes$constructors$distrib distribution_node <- .internals$nodes$node_classes$distribution_node check_dims <- .internals$utils$checks$check_dims as.greta_array <- .internals$greta_arrays$as.greta_array is_scalar <- .internals$utils$misc$is_scalar fl <- .internals$utils$misc$fl
ld <- function (i, M, sigma, dim = 1) { distrib("ld", i, M, sigma, dim) }
ld_distribution <- R6Class ( "ld_distribution", inherit = distribution_node, public = list(
initialize = function (i, M, sigma, dim) {
# check if (i, M, dim) are in counting set
# (i)
if (length(i) > 1 ||
i <= 0 ||
!is.finite(i) ||
i != floor(i)) {
stop ("i must be a scalar positive integer, but was: ",
capture.output(dput(i)),
call. = FALSE)
}
# (M)
if (length(M) > 1 ||
M <= 0 ||
!is.finite(M) ||
M != floor(M)) {
stop ("M must be a scalar positive integer, but was: ",
capture.output(dput(M)),
call. = FALSE)
}
# (dim)
if (length(dim) > 1 ||
dim <= 0 ||
!is.finite(dim) ||
dim != floor(dim)) {
stop ("dim must be a scalar positive integer, but was: ",
capture.output(dput(dim)),
call. = FALSE)
}
# check if i > M
if (i > M) {
stop ("i can be no larger than M",
call. = FALSE)
}
# check if sigma is scalar
# (sigma)
if (length(sigma) > 1) {
stop ("sigma must be a scalar positive integer, but was: ",
capture.output(dput(sigma)),
call. = FALSE)
}
i <- as.greta_array(i)
M <- as.greta_array(M)
sigma <- as.greta_array(sigma)
self$bounds <- c(0, Inf)
super$initialize("ld", dim, truncation = c(0, Inf))
self$add_parameter(i, "i")
self$add_parameter(M, "M")
self$add_parameter(sigma, "sigma")
},
tf_distrib = function (parameters, dag) {
i <- parameters$i
M <- parameters$M
sigma <- parameters$sigma
# log pdf(x | i, M, sigma)
log_prob = function (x) {
(M - i) * tf$log(x) - x ^ fl(2) / (fl(2) * sigma)
}
list(log_prob = log_prob, cdf = NULL, log_cdf = NULL)
},
# no CDF for discrete distributions
tf_cdf_function = NULL,
tf_log_cdf_function = NULL) )
On Wed, Mar 23, 2022 at 9:44 AM Nicholas Tierney @.***> wrote:
Hi there!
Thanks for posting here, I was wondering if you can make the following file accessible so I can reproduce it?
source("~/Dropbox/Econometrics/R/Functions/LD_distribution/ld_distro.r")
One thing that can be really helpful with checking if it is a standalone reproducible example is to run the code using reprex: https://github.com/tidyverse/reprex
So you could write:
reprex({#
}) — Reply to this email directly, view it on GitHub https://github.com/greta-dev/greta/issues/509#issuecomment-1076095637, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHLXB32NPDELRWRNPQ5MVSTVBLKXTANCNFSM5RHF3L2A . You are receiving this because you authored the thread.Message ID: @.***>
BillyChen123
commented
1 year ago Hi! @giannilmbd
Sorry to trouble you. I am running the same code just now. Then when I run the code
head(calculate(ZZ_est,values = list(ZZ_est_raw_d = rep(1, sum(idx_d)), ZZ_est_raw_s = rep(2, sum(idx_s)))))
The result shows that the ld distribution is not yet sampling. Is that the problem of function R6? What should I do next?
giannilmbd
commented
1 year ago Hi, unfortunately I was simply trying to use that function (apparently unsuccessfully). I have no clue about the details. Best Gianni
On Thu, 20 Oct 2022 at 04:24, BillyChen123 @.***> wrote:
Hi! @giannilmbd https://github.com/giannilmbd Sorry to trouble you. I am running the same code just now. Then when I run the code head(calculate(ZZ_est,values = list(ZZ_est_raw_d = rep(1, sum(idx_d)), ZZ_est_raw_s = rep(2, sum(idx_s))))) The result shows that the ld distribution is not yet sampling. Is that the problem of function R6? What should I do next?
— Reply to this email directly, view it on GitHub https://github.com/greta-dev/greta/issues/509#issuecomment-1284820347, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHLXB33VJKBYFBERDZ6Y2OLWECUMTANCNFSM5RHF3L2A . You are receiving this because you were mentioned.Message ID: @.***>
-- Giovanni Lombardo Bank for International Settlements Basel - CH www.giovannilombardo.ch @.***
njtierney
commented
1 year ago Hi there @giannilmbd
I got this to reproduce, however there was an issue doing this in parallel, which I think is related to the distribution not being packaged up.
I suspect if we add the distribution to greta, perhaps in the new greta.distributions package, this might help resolve this?
For your use, here is the reprex:
library(greta)
#>
#> Attaching package: 'greta'
#> The following objects are masked from 'package:stats':
#>
#> binomial, cov2cor, poisson
#> The following objects are masked from 'package:base':
#>
#> %*%, apply, backsolve, beta, chol2inv, colMeans, colSums, diag,
#> eigen, forwardsolve, gamma, identity, rowMeans, rowSums, sweep,
#> tapply
# using latest version of greta, installed with
# remotes::install_github("njtierney/greta@tf2-poke-tf-fun")
# or
# remotes::install_github("https://github.com/njtierney/greta/tree/tf2-poke-tf-fun")
library(tensorflow)
library(MASS)
library(viridis)
#> Loading required package: viridisLite
# We need to define a helper function to fill in the diagonal of the prior for the loadings matrix Z with the appropriate distribution (see below).
zd <- function(M, sigma) {
zd <- zeros(M)
for (i in 1:M) {
zd[i] <- ld(i, M = M, sigma = sigma, dim = 1)
}
return(zd)
}
NN <- 30
TT <- 30
MM <- 3
## simulation
# latent factors
set.seed(123)
## MM x TT matrix of innovations
ww <- matrix(rnorm(MM * TT, 0, 1), MM, TT)
ww[, 1] <- rnorm(MM, 0, sqrt(5))
## MM x TT matrix of scaled latent trends
xx <- t(scale(apply(ww, 1, cumsum)))
head(xx)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] -1.43843703 -1.404926 -1.1858634 -1.3976761 -1.2071988 -0.3579210
#> [2,] 1.31267329 1.396602 0.5753723 1.3699997 1.4418507 1.7650363
#> [3,] -0.09420896 1.101622 0.6227131 0.8735937 0.4860328 -0.8851932
#> [,7] [,8] [,9] [,10] [,11] [,12]
#> [1,] -0.02458304 -0.1281814 -0.4252479 -0.3523532 -0.14966481 0.2676914
#> [2,] 1.45811809 0.7920750 -0.3028622 -1.0416974 -1.23324661 -0.6999074
#> [3,] -1.62973443 -2.1379547 -1.5538067 -0.6795829 -0.05545549 0.4246997
#> [,13] [,14] [,15] [,16] [,17] [,18]
#> [1,] 0.5309554 0.35012651 -0.2512866 -0.78507373 -0.41437465 -0.4279423
#> [2,] -0.7400981 -1.19107549 0.2169283 -0.04460925 -0.09872927 -0.1265591
#> [3,] 0.2113669 0.06639633 0.9086492 0.58327322 0.75989983 1.7141593
#> [,19] [,20] [,21] [,22] [,23] [,24]
#> [1,] -0.5352459 -0.2573927 -0.0769590 -0.56106393 -0.34804036 0.6263167
#> [2,] 0.8578759 0.9382774 0.6121881 -0.08357802 -0.04916969 -0.3679285
#> [3,] 0.6342897 0.7848552 0.5525260 0.76416164 1.40721366 -0.2028569
#> [,25] [,26] [,27] [,28] [,29] [,30]
#> [1,] 1.1043206 1.591751 1.677920 1.861035 1.7562424 1.9630740
#> [2,] -0.8283146 -1.013178 -1.103341 -1.343960 -1.1285796 -1.3401620
#> [3,] -0.6825717 -1.533718 -1.529699 -1.080407 -0.3156349 0.4853725
dim(xx)
#> [1] 3 30
matplot(1:dim(xx)[2], t(xx), "l")
# loading matrix
ZZ <- matrix(runif(NN * MM, -1, 1), NN, MM)
diag(ZZ) <- rev(sort(abs(diag(ZZ))))
ZZ[upper.tri(ZZ)] <- 0
ZZ <- round(ZZ, 2)
# observed series
## obs var
obs_var <- 0.2^2
## obs errors
ee <- t(mvrnorm(TT, matrix(0, NN, 1), diag(obs_var, NN, NN)))
## NN x TT matrix of observed data
yy <- ZZ %*% xx + ee
matplot(1:dim(yy)[2], t(yy), "l")
## ESTIMATIONlibrary(R6)
distrib <- .internals$nodes$constructors$distrib
distribution_node <- .internals$nodes$node_classes$distribution_node
check_dims <- .internals$utils$checks$check_dims
as.greta_array <- .internals$greta_arrays$as.greta_array
is_scalar <- .internals$utils$misc$is_scalar
fl <- .internals$utils$misc$fl
# function to call
ld <- function(i, M, sigma, dim = 1) {
distrib("ld", i, M, sigma, dim)
}
# Leung & Drton distn for prior of diag(Z)
ld_distribution <- R6Class(
"ld_distribution",
inherit = distribution_node,
public = list(
initialize = function(i, M, sigma, dim) {
# check if (i, M, dim) are in counting set
# (i)
if (length(i) > 1 ||
i <= 0 ||
!is.finite(i) ||
i != floor(i)) {
stop(
"i must be a scalar positive integer, but was: ",
capture.output(dput(i)),
call. = FALSE
)
}
# (M)
if (length(M) > 1 ||
M <= 0 ||
!is.finite(M) ||
M != floor(M)) {
stop(
"M must be a scalar positive integer, but was: ",
capture.output(dput(M)),
call. = FALSE
)
}
# (dim)
if (length(dim) > 1 ||
dim <= 0 ||
!is.finite(dim) ||
dim != floor(dim)) {
stop(
"dim must be a scalar positive integer, but was: ",
capture.output(dput(dim)),
call. = FALSE
)
}
# check if i > M
if (i > M) {
stop(
"i can be no larger than M",
call. = FALSE
)
}
# check if sigma is scalar
# (sigma)
if (length(sigma) > 1) {
stop(
"sigma must be a scalar positive integer, but was: ",
capture.output(dput(sigma)),
call. = FALSE
)
}
i <- as.greta_array(i)
M <- as.greta_array(M)
sigma <- as.greta_array(sigma)
self$bounds <- c(0, Inf)
super$initialize("ld", dim, truncation = c(0, Inf))
self$add_parameter(i, "i")
self$add_parameter(M, "M")
self$add_parameter(sigma, "sigma")
},
tf_distrib = function(parameters, dag) {
i <- parameters$i
M <- parameters$M
sigma <- parameters$sigma
# log pdf(x | i, M, sigma)
log_prob <- function(x) {
(M - i) * tf$math$log(x) - x^fl(2) / (fl(2) * sigma)
}
list(log_prob = log_prob, cdf = NULL, log_cdf = NULL)
},
# no CDF for discrete distributions
tf_cdf_function = NULL,
tf_log_cdf_function = NULL
)
)
## PRIORS
## empty loadings matrix
ZZ_est <- zeros(NN, MM)
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#>
## define sigma
sigma_est <- normal(0, 2, truncation = c(0, Inf))
## diagonal
idx_d <- row(ZZ_est) == col(ZZ_est)
ZZ_est_raw_d <- zd(MM, sigma_est)
ZZ_est[idx_d] <- ZZ_est_raw_d
## sub-diagonal
idx_s <- lower.tri(ZZ_est)
ZZ_est_raw_s <- normal(
0,
sigma_est,
dim = sum(idx_s),
truncation = c(
-1,
1
)
)
ZZ_est[idx_s] <- ZZ_est_raw_s
# check that everything works
# head(calculate(ZZ_est, values=list(ZZ_est_raw_d = rep(1, sum(idx_d)),
# ZZ_est_raw_s = rep(2, sum(idx_s)))))
# COVARIANCE MATRIX
# diagonal
RR_est <- inverse_gamma(alpha = 1, beta = 3, 1, truncation = c(0, 1))
## inital factor
X0 <- normal(mean = 0, sd = sqrt(10), dim = c(MM, 1))
## Likerlihood function
## factors for t = 2-TT
XX <- normal(mean = 0, sd = 1, dim = c(MM, TT))
## cumsum over proc errs for the random walk
xx_est <- t(apply(cbind(X0, XX), 1, "cumsum"))[, -1]
# prepare observed data
## scale data
yy_z <- t(scale(t(yy), scale = FALSE))
## vectorize data
yy_vec <- as_data(matrix(yy_z, 1, NN * TT))
## vectorize mean
Zx_vec <- t(c(ZZ_est %*% xx_est))
## define likelihood
distribution(yy_vec) <- normal(Zx_vec, RR_est)
### FIT THE MODEL WITH MCMC
# define the model for greta
mod_defn <- model(xx_est, ZZ_est, RR_est, sigma_est)
#### mcmc
# require(future) # all options of plan crash (except sequential)
# cores=detectCores()-1
# cl1 <- makeCluster(cores)
# # registerDoParallel(cl1)
# plan(multisession) # the cluster plan conflicts with TensorFlow
# model works for me
mod_fit <- mcmc(
mod_defn,
verbose = TRUE,
# n_cores = 8,
sampler = hmc(Lmin = 1, Lmax = 30, epsilon = 0.001, diag_sd = 1),
warmup = 10,
n_samples = 10,
# thin = 10,
chains = 1,
initial_values = initials(
RR_est = runif(1, 0.1, 1),
sigma_est = runif(1, 0.1, 1)
)
)
#>
#> warmup 0/10 | eta: ?s warmup ========================================== 10/10 | eta: 0s
#> sampling 0/10 | eta: ?s sampling ========================================== 10/10 | eta: 0s
## now trying with future - it fails, most likely due to an issue with
## having the distribution loaded properly
library(future)
plan(multisession)
mod_fit_future <- mcmc(
mod_defn,
verbose = TRUE,
# n_cores = 8,
sampler = hmc(Lmin = 1, Lmax = 30, epsilon = 0.001, diag_sd = 1),
warmup = 10,
n_samples = 10,
# thin = 10,
chains = 2,
initial_values = initials(
RR_est = runif(1, 0.1, 1),
sigma_est = runif(1, 0.1, 1)
)
)
#> only one set of initial values was provided, and was used for all chains
#> running 2 samplers in parallel, on up to 4 CPU cores
#> warmup 0/10 | eta: ?s ...
#> Error: greta hit a tensorflow error:
#> Error in eval(expr, p): RuntimeError: RuntimeError: RuntimeError: object 'tf'
#> not foundCreated on 2023-04-19 with reprex v2.0.2
giannilmbd
commented
1 year ago Hi Nicholas, thanks for following up on this. I‘ll give it a try once I get back to the problem. Cheers Gianni
On Wed, 19 Apr 2023 at 03:58, Nicholas Tierney @.***> wrote:
Hi there @giannilmbd https://github.com/giannilmbd
I got this to reproduce, however there was an issue doing this in parallel, which I think is related to the distribution not being packaged up.
I suspect if we add the distribution to greta, perhaps in the new greta.distributions package, this might help resolve this?
For your use, here is the reprex:
library(greta)#> #> Attaching package: 'greta'#> The following objects are masked from 'package:stats':#> #> binomial, cov2cor, poisson#> The following objects are masked from 'package:base':#> #> %*%, apply, backsolve, beta, chol2inv, colMeans, colSums, diag,#> eigen, forwardsolve, gamma, identity, rowMeans, rowSums, sweep,#> tapply# using latest version of greta, installed with# @.***")# or # remotes::install_github("https://github.com/njtierney/greta/tree/tf2-poke-tf-fun") library(tensorflow) library(MASS) library(viridis)#> Loading required package: viridisLite
We need to define a helper function to fill in the diagonal of the prior for the loadings matrix Z with the appropriate distribution (see below).
zd <- function(M, sigma) { zd <- zeros(M) for (i in 1:M) { zd[i] <- ld(i, M = M, sigma = sigma, dim = 1) } return(zd) }NN <- 30TT <- 30MM <- 3## simulation# latent factors
set.seed(123)## MM x TT matrix of innovationsww <- matrix(rnorm(MM * TT, 0, 1), MM, TT)ww[, 1] <- rnorm(MM, 0, sqrt(5))## MM x TT matrix of scaled latent trendsxx <- t(scale(apply(ww, 1, cumsum))) head(xx)#> [,1] [,2] [,3] [,4] [,5] [,6]#> [1,] -1.43843703 -1.404926 -1.1858634 -1.3976761 -1.2071988 -0.3579210#> [2,] 1.31267329 1.396602 0.5753723 1.3699997 1.4418507 1.7650363#> [3,] -0.09420896 1.101622 0.6227131 0.8735937 0.4860328 -0.8851932#> [,7] [,8] [,9] [,10] [,11] [,12]#> [1,] -0.02458304 -0.1281814 -0.4252479 -0.3523532 -0.14966481 0.2676914#> [2,] 1.45811809 0.7920750 -0.3028622 -1.0416974 -1.23324661 -0.6999074#> [3,] -1.62973443 -2.1379547 -1.5538067 -0.6795829 -0.05545549 0.4246997#> [,13] [,14] [,15] [,16] [,17] [,18]#> [1,] 0.5309554 0.35012651 -0.2512866 -0.78507373 -0.41437465 -0.4279423#> [2,] -0.7400981 -1.19107549 0.2169283 -0.04460925 -0.09872927 -0.1265591#> [3,] 0.2113669 0.06639633 0.9086492 0.58327322 0.75989983 1.7141593#> [,19] [,20] [,21] [,22] [,23] [,24]#> [1,] -0.5352459 -0.2573927 -0.0769590 -0.56106393 -0.34804036 0.6263167#> [2,] 0.8578759 0.9382774 0.6121881 -0.08357802 -0.04916969 -0.3679285#> [3,] 0.6342897 0.7848552 0.5525260 0.76416164 1.40721366 -0.2028569#> [,25] [,26] [,27] [,28] [,29] [,30]#> [1,] 1.1043206 1.591751 1.677920 1.861035 1.7562424 1.9630740#> [2,] -0.8283146 -1.013178 -1.103341 -1.343960 -1.1285796 -1.3401620#> [3,] -0.6825717 -1.533718 -1.529699 -1.080407 -0.3156349 0.4853725 dim(xx)#> [1] 3 30 matplot(1:dim(xx)[2], t(xx), "l")
loading matrixZZ <- matrix(runif(NN * MM, -1, 1), NN, MM)
diag(ZZ) <- rev(sort(abs(diag(ZZ))))ZZ[upper.tri(ZZ)] <- 0ZZ <- round(ZZ, 2)
observed series## obs varobs_var <- 0.2^2## obs errorsee <- t(mvrnorm(TT, matrix(0, NN, 1), diag(obs_var, NN, NN)))## NN x TT matrix of observed datayy <- ZZ %*% xx + ee
matplot(1:dim(yy)[2], t(yy), "l")
ESTIMATION
library(R6)distrib <- .internals$nodes$constructors$distribdistribution_node <- .internals$nodes$node_classes$distribution_nodecheck_dims <- .internals$utils$checks$check_dimsas.greta_array <- .internals$greta_arrays$as.greta_arrayis_scalar <- .internals$utils$misc$is_scalarfl <- .internals$utils$misc$fl
function to callld <- function(i, M, sigma, dim = 1) {
distrib("ld", i, M, sigma, dim) }
Leung & Drton distn for prior of diag(Z)ld_distribution <- R6Class(
"ld_distribution", inherit = distribution_node, public = list( initialize = function(i, M, sigma, dim) {
check if (i, M, dim) are in counting set
# (i) if (length(i) > 1 || i <= 0 || !is.finite(i) || i != floor(i)) { stop( "i must be a scalar positive integer, but was: ", capture.output(dput(i)), call. = FALSE ) } # (M) if (length(M) > 1 || M <= 0 || !is.finite(M) || M != floor(M)) { stop( "M must be a scalar positive integer, but was: ", capture.output(dput(M)), call. = FALSE ) } # (dim) if (length(dim) > 1 || dim <= 0 || !is.finite(dim) || dim != floor(dim)) { stop( "dim must be a scalar positive integer, but was: ", capture.output(dput(dim)), call. = FALSE ) } # check if i > M if (i > M) { stop( "i can be no larger than M", call. = FALSE ) } # check if sigma is scalar # (sigma) if (length(sigma) > 1) { stop( "sigma must be a scalar positive integer, but was: ", capture.output(dput(sigma)), call. = FALSE ) } i <- as.greta_array(i) M <- as.greta_array(M) sigma <- as.greta_array(sigma) self$bounds <- c(0, Inf) super$initialize("ld", dim, truncation = c(0, Inf)) self$add_parameter(i, "i") self$add_parameter(M, "M") self$add_parameter(sigma, "sigma") }, tf_distrib = function(parameters, dag) { i <- parameters$i M <- parameters$M sigma <- parameters$sigma
# log pdf(x | i, M, sigma) log_prob <- function(x) { (M - i) * tf$math$log(x) - x^fl(2) / (fl(2) * sigma) } list(log_prob = log_prob, cdf = NULL, log_cdf = NULL) }, # no CDF for discrete distributions tf_cdf_function = NULL, tf_log_cdf_function = NULL) )
PRIORS## empty loadings matrixZZ_est <- zeros(NN, MM)#> ℹ Initialising python and checking dependencies, this may take a moment.#> ✔ Initialising python and checking dependencies ... done!#> ## define sigmasigma_est <- normal(0, 2, truncation = c(0, Inf))## diagonalidx_d <- row(ZZ_est) == col(ZZ_est)ZZ_est_raw_d <- zd(MM, sigma_est)ZZ_est[idx_d] <- ZZ_est_raw_d## sub-diagonalidx_s <- lower.tri(ZZ_est)ZZ_est_raw_s <- normal(
0, sigma_est, dim = sum(idx_s), truncation = c( -1, 1 ) )ZZ_est[idx_s] <- ZZ_est_raw_s
check that everything works# head(calculate(ZZ_est, values=list(ZZ_est_raw_d = rep(1, sum(idx_d)),# ZZ_est_raw_s = rep(2, sum(idx_s)))))
COVARIANCE MATRIX# diagonalRR_est <- inverse_gamma(alpha = 1, beta = 3, 1, truncation = c(0, 1))
inital factorX0 <- normal(mean = 0, sd = sqrt(10), dim = c(MM, 1))## Likerlihood function## factors for t = 2-TTXX <- normal(mean = 0, sd = 1, dim = c(MM, TT))## cumsum over proc errs for the random walkxx_est <- t(apply(cbind(X0, XX), 1, "cumsum"))[, -1]
prepare observed data## scale datayy_z <- t(scale(t(yy), scale = FALSE))## vectorize datayy_vec <- as_data(matrix(yy_z, 1, NN TT))## vectorize meanZx_vec <- t(c(ZZ_est %% xx_est))## define likelihood
distribution(yy_vec) <- normal(Zx_vec, RR_est)
FIT THE MODEL WITH MCMC# define the model for gretamod_defn <- model(xx_est, ZZ_est, RR_est, sigma_est)#### mcmc# require(future) # all options of plan crash (except sequential)# cores=detectCores()-1# cl1 <- makeCluster(cores)# # registerDoParallel(cl1)# plan(multisession) # the cluster plan conflicts with TensorFlow# model works for memod_fit <- mcmc(
mod_defn, verbose = TRUE,
n_cores = 8,
sampler = hmc(Lmin = 1, Lmax = 30, epsilon = 0.001, diag_sd = 1), warmup = 10, n_samples = 10,
thin = 10,
chains = 1, initial_values = initials( RR_est = runif(1, 0.1, 1), sigma_est = runif(1, 0.1, 1) ) )#> #> warmup 0/10 | eta: ?s warmup ========================================== 10/10 | eta: 0s #> sampling 0/10 | eta: ?s sampling ========================================== 10/10 | eta: 0s
now trying with future - it fails, most likely due to an issue with## having the distribution loaded properly
library(future) plan(multisession) mod_fit_future <- mcmc( mod_defn, verbose = TRUE,
n_cores = 8,
sampler = hmc(Lmin = 1, Lmax = 30, epsilon = 0.001, diag_sd = 1), warmup = 10, n_samples = 10,
thin = 10,
chains = 2, initial_values = initials( RR_est = runif(1, 0.1, 1), sigma_est = runif(1, 0.1, 1) ) )#> only one set of initial values was provided, and was used for all chains#> running 2 samplers in parallel, on up to 4 CPU cores#> warmup 0/10 | eta: ?s ...#> Error: greta hit a tensorflow error:#> Error in eval(expr, p): RuntimeError: RuntimeError: RuntimeError: object 'tf'#> not found
Created on 2023-04-19 with reprex v2.0.2 https://reprex.tidyverse.org Session info
sessioninfo::session_info()#> ─ Session info ───────────────────────────────────────────────────────────────#> setting value#> version R version 4.2.3 (2023-03-15)#> os macOS Ventura 13.2#> system aarch64, darwin20#> ui X11#> language (EN)#> collate en_US.UTF-8#> ctype en_US.UTF-8#> tz Australia/Brisbane#> date 2023-04-19#> pandoc 2.19.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)#> #> ─ Packages ───────────────────────────────────────────────────────────────────#> package version date (UTC) lib source#> abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)#> backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)#> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.2.0)#> callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.0)#> cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.0)#> coda 0.19-4 2020-09-30 [1] CRAN (R 4.2.0)#> codetools 0.2-19 2023-02-01 [1] CRAN (R 4.2.3)#> colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.0)#> crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.0)#> curl 5.0.0 2023-01-12 [1] CRAN (R 4.2.0)#> digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.0)#> dplyr 1.1.1 2023-03-22 [1] CRAN (R 4.2.0)#> evaluate 0.20 2023-01-17 [1] CRAN (R 4.2.0)#> fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.0)#> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.0)#> fs 1.6.1 2023-02-06 [1] CRAN (R 4.2.0)#> future 1.32.0 2023-03-07 [1] CRAN (R 4.2.0)#> generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)#> ggplot2 3.4.2 2023-04-03 [1] CRAN (R 4.2.0)#> globals 0.16.2 2022-11-21 [1] CRAN (R 4.2.1)#> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)#> greta 0.4.3.9000 2023-04-19 [1] local#> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.2.0)#> gtable 0.3.3 2023-03-21 [1] CRAN (R 4.2.0)#> here 1.0.1 2020-12-13 [1] CRAN (R 4.2.0)#> highr 0.10 2022-12-22 [1] CRAN (R 4.2.0)#> hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.0)#> htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.0)#> httr 1.4.5 2023-02-24 [1] CRAN (R 4.2.0)#> jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.0)#> knitr 1.42 2023-01-25 [1] CRAN (R 4.2.0)#> lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.0)#> lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.0)#> listenv 0.9.0 2022-12-16 [1] CRAN (R 4.2.0)#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)#> MASS 7.3-58.3 2023-03-07 [1] CRAN (R 4.2.0)#> Matrix 1.5-4 2023-04-04 [1] CRAN (R 4.2.0)#> mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)#> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)#> parallelly 1.35.0 2023-03-23 [1] CRAN (R 4.2.0)#> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.0)#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)#> png 0.1-8 2022-11-29 [1] CRAN (R 4.2.0)#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.0)#> processx 3.8.0 2022-10-26 [1] CRAN (R 4.2.0)#> progress 1.2.2 2019-05-16 [1] CRAN (R 4.2.0)#> ps 1.7.4 2023-04-02 [1] CRAN (R 4.2.0)#> purrr 1.0.1 2023-01-10 [1] CRAN (R 4.2.0)#> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.0)#> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.0)#> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.0)#> R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.2.0)#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)#> Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.0)#> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.0)#> reticulate 1.28 2023-01-27 [1] CRAN (R 4.2.0)#> rlang 1.1.0 2023-03-14 [1] CRAN (R 4.2.0)#> rmarkdown 2.21 2023-03-26 [1] CRAN (R 4.2.0)#> rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.2.0)#> rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.0)#> scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)#> styler 1.9.1 2023-03-04 [1] CRAN (R 4.2.0)#> tensorflow 2.11.0 2022-12-19 [1] CRAN (R 4.2.0)#> tfautograph 0.3.2 2021-09-17 [1] CRAN (R 4.2.0)#> tfruns 1.5.1 2022-09-05 [1] CRAN (R 4.2.0)#> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.2.0)#> tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.0)#> utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.0)#> vctrs 0.6.1 2023-03-22 [1] CRAN (R 4.2.0)#> viridis 0.6.2 2021-10-13 [1] CRAN (R 4.2.0)#> viridisLite 0.4.1 2022-08-22 [1] CRAN (R 4.2.0)#> whisker 0.4.1 2022-12-05 [1] CRAN (R 4.2.0)#> withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)#> xfun 0.38 2023-03-24 [1] CRAN (R 4.2.0)#> xml2 1.3.3 2021-11-30 [1] CRAN (R 4.2.0)#> yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.0)#> #> [1] /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library#> #> ─ Python configuration ───────────────────────────────────────────────────────#> python: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/bin/python#> libpython: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/libpython3.8.dylib#> pythonhome: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2:/Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2#> version: 3.8.15 | packaged by conda-forge | (default, Nov 22 2022, 08:49:06) [Clang 14.0.6 ]#> numpy: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/python3.8/site-packages/numpy#> numpy_version: 1.23.2#> tensorflow: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/python3.8/site-packages/tensorflow#> #> NOTE: Python version was forced by use_python function#> #> ──────────────────────────────────────────────────────────────────────────────
— Reply to this email directly, view it on GitHub https://github.com/greta-dev/greta/issues/509#issuecomment-1514027015, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHLXB3ZIZLNO2L4IJVZTF2TXB5BDTANCNFSM5RHF3L2A . You are receiving this because you were mentioned.Message ID: @.***>
-- Giovanni Lombardo Bank for International Settlements Basel - CH www.giovannilombardo.ch @.***
Hi, Thanks a lot for this amazing package!
I'm trying to use greta's mcmc on a Linux (Ubuntu) machine in RStudio 2022.01.0 Build 205 and R 4.1.3. To check if everything works I've copied the example in https://mdscheuerell.github.io/gretaDFA/
When I try the following command:
future::plan(cluster) # or future::plan(multisession) mod_fit <- mcmc(mod_defn, verbose = TRUE,n_cores=8, sampler = hmc(Lmin = 1, Lmax = 30, epsilon = 0.001, diag_sd = 1), warmup = 2000, n_samples = 5000, thin = 10, chains = 1, initial_values = initials(RR_est = runif(1, 0.1, 1), sigma_est = runif(1, 0.1, 1) ) )
I get the error
Loaded Tensorflow version 1.14.0 Error in py_call_impl(callable, dots$args, dots$keywords) : RuntimeError: Evaluation error: object 'tf' not found.
Detailed traceback: File "/home/gianni/.local/share/r-miniconda/envs/greta-env/lib/python3.7/site-packages/tensorflow_probability/python/mcmc/sample.py", line 326, in sample_chain previous_kernel_results = kernel.bootstrap_results(current_state)
Any suggestion on how to fix this problem? Thanks Gianni
Reproducible chunk