Greta_mcmc_list not working with saveRDS?

#1

I want to save a mcmc() output so I do not need to rerun it next time. However, using saveRDS() (or save()) seems to decouple the greta_mcmc_list object from other greta_arrays.

Example code:

First, run this:

library(greta)

alpha <- normal(0, 1)
beta <- normal(0, 1)
sigma <- lognormal(1, 0.1)
y <- as_data(iris$Petal.Width)
mu <- alpha + iris$Petal.Length * beta
distribution(y) <- normal(mu, sigma)
m <- model(alpha, beta, sigma)
draws <- mcmc(m, n_samples = 500)

# save draws as rds so do not need to rerun mcmc every time
saveRDS(draws, file = "draws.rds")

Then restart the session and load the greta_mcmc_list object. This will fail halfway:

library(greta)

draws <- readRDS(file = "draws.rds")

# the codes below breaks 
petal_length_plot <- seq(min(iris$Petal.Length),
                         max(iris$Petal.Length),
                         length.out = 100)
mu_plot <- alpha + petal_length_plot * beta
mu_plot_draws <- calculate(mu_plot, values = draws)
mu_est <- colMeans(mu_plot_draws[[1]])
plot(mu_est ~ petal_length_plot, type = "n",
     ylim = range(mu_plot_draws[[1]]))
apply(mu_plot_draws[[1]], 1, lines,
      x = petal_length_plot, col = grey(0.8))
lines(mu_est ~ petal_length_plot, lwd = 2)

You should get an error:

Error: the target greta arrays do not appear to be connected to those in the greta_mcmc_list object

Did I do something wrong?

#2

Hi all, I have a solution that seems to work, though not sure if it’s the best.

In the question above, I only saved the draws object. The model object m was not saved, so in a new session R doesn’t know what the target and/or visible arrays are.

So when I also saved the m object, and then load it together with draws.rds in the new session, it works.

#3

Hi there,

Thanks for posting an issue and following up with an answer!

I tried your example but noted that I needed to save all of the objects and then write them back out - it looks like this might be something you’ve already solved but in case it is helpful here is some code I wrote to reproduce the above.

# this breaks because only the "draws" object has been created, it needs
# to know about alpha, beta, and all the other objects.
library(greta)

alpha <- normal(0, 1)
beta <- normal(0, 1)
sigma <- lognormal(1, 0.1)
y <- as_data(iris$Petal.Width)
mu <- alpha + iris$Petal.Length * beta
distribution(y) <- normal(mu, sigma)
m <- model(alpha, beta, sigma)
draws <- mcmc(m, n_samples = 500)

library(readr)
# save draws as rds so do not need to rerun mcmc every time
write_rds(draws, file = "mcmc-list-not-saving/draws.rds")
write_rds(alpha, file = "mcmc-list-not-saving/alpha.rds")
write_rds(beta, file = "mcmc-list-not-saving/beta.rds")
write_rds(sigma, file = "mcmc-list-not-saving/sigma.rds")
write_rds(y, file = "mcmc-list-not-saving/y.rds")
write_rds(mu, file = "mcmc-list-not-saving/mu.rds")
write_rds(m, file = "mcmc-list-not-saving/m.rds")

Then restart R

library(greta)
library(readr)
draws <- read_rds(file = "mcmc-list-not-saving/draws.rds")
alpha <- read_rds(file = "mcmc-list-not-saving/alpha.rds")
beta <- read_rds(file = "mcmc-list-not-saving/beta.rds")
sigma <- read_rds(file = "mcmc-list-not-saving/sigma.rds")
y <- read_rds(file = "mcmc-list-not-saving/y.rds")
mu <- read_rds(file = "mcmc-list-not-saving/mu.rds")
m <- read_rds(file = "mcmc-list-not-saving/m.rds")

petal_length_plot <- seq(min(iris$Petal.Length),
                         max(iris$Petal.Length),
                         length.out = 100)

mu_plot <- alpha + petal_length_plot * beta

mu_plot_draws <- calculate(mu_plot, values = draws)
mu_est <- colMeans(mu_plot_draws[[1]])
plot(mu_est ~ petal_length_plot, type = "n",
     ylim = range(mu_plot_draws[[1]]))
apply(mu_plot_draws[[1]], 1, lines,
      x = petal_length_plot, col = grey(0.8))
lines(mu_est ~ petal_length_plot, lwd = 2)

image
image1510×977 26.8 KB

A small note on workflows.

I’m working on creating some example workflows using targets, which helps facilitate this process you went through - of running code once and then saving it and using it later. I’ll post an example onto this page when I’m done.

1 Like