Growth model with t-distribution: merge samples

Load packages, read data and source custom scripts

rm(list = ls())
library(bamlss)

#> Loading required package: coda

#> Loading required package: colorspace

#> Loading required package: mgcv

#> Loading required package: nlme

#> This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

#> 
#> Attaching package: 'bamlss'

#> The following object is masked from 'package:mgcv':
#> 
#>     smooth.construct

library(gamlss.dist)

#> Loading required package: MASS

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_processed <- file.path(path_data, "processed")
path_modelled <- file.path(path_data, "modelled")

source(file.path(path_proj, "src", "47-mcmc-merge-samples.R"))

path_modelled_data <- file.path(path_modelled, "bw-20-growth-re-t.rds")
path_modelled_continue <- gsub("(\\.rds)$", "-continue\\1", path_modelled_data)
path_modelled_continue2 <- gsub("(\\.rds)$", "-continue2\\1", path_modelled_data)
path_modelled_merged <- gsub("(\\.rds)$", "-merged\\1", path_modelled_data)

model <- readRDS(path_modelled_data)
model_continue <- readRDS(path_modelled_continue)
model_continue2 <- readRDS(path_modelled_continue2)

Merge samples

model$samples <- merge_samples(model_continue$samples, model_continue2$samples)
system.time(saveRDS(model, file = path_modelled_merged))

#>    user  system elapsed 
#>  11.925   0.056  12.076

Maximum auto-correlation function (ACF)

par(mar = c(4, 4, 0.5, 0), mfrow = c(1, 3))
plot(model, model = "mu", which = "max-acf", spar = FALSE)
plot(model, model = "sigma", which = "max-acf", spar = FALSE)
plot(model, model = "nu", which = "max-acf", spar = FALSE)

Figure 1: Maximum ACF of samples for \(\mu\) (left), \(\sigma\) (center) and \(\nu\)(right)

MCMC convergence

par(mar = c(4, 4, 3, 1), mfrow = c(1, 2))
plot(model, which = "samples")

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#> 630.065   1.215 637.373