# this didn't work in MCMC, we'll have to work out why.
# simulate species count data over 10 timesteps
n_times <- 10
times <- seq_len(n_times)
truth <- 100 * exp(sin(times))
y <- rpois(n_times, truth)
plot(y ~ times, type = "l")
library(greta.gam)
#> Loading required package: 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
library(greta.gp)
eta <- smooths(~ 1 + s(time), data = data.frame(time = times))
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#>
lambda <- exp(eta)
distribution(y) <- poisson(lambda)
len <- lognormal(0, 1)
var <- normal(0, 1, truncation = c(0, Inf))
k <- rbf(lengthscales = len, variance = var)
eta <- gp(times, k)
lambda <- exp(eta)
distribution(y) <- poisson(lambda)
m <- model(eta)
draws <- mcmc(m)
#> running 4 chains simultaneously on up to 8 CPU cores
#>
#> warmup 0/1000 | eta: ?s warmup == 50/1000 | eta: 26s warmup ==== 100/1000 | eta: 15s warmup ====== 150/1000 | eta: 11s warmup ======== 200/1000 | eta: 9s warmup ========== 250/1000 | eta: 7s warmup =========== 300/1000 | eta: 6s warmup ============= 350/1000 | eta: 5s warmup =============== 400/1000 | eta: 5s warmup ================= 450/1000 | eta: 4s warmup =================== 500/1000 | eta: 4s warmup ===================== 550/1000 | eta: 3s warmup ======================= 600/1000 | eta: 3s warmup ========================= 650/1000 | eta: 2s warmup =========================== 700/1000 | eta: 2s warmup ============================ 750/1000 | eta: 2s warmup ============================== 800/1000 | eta: 1s warmup ================================ 850/1000 | eta: 1s warmup ================================== 900/1000 | eta: 1s warmup ==================================== 950/1000 | eta: 0s warmup ====================================== 1000/1000 | eta: 0s
#> sampling 0/1000 | eta: ?s sampling == 50/1000 | eta: 3s sampling ==== 100/1000 | eta: 3s sampling ====== 150/1000 | eta: 3s sampling ======== 200/1000 | eta: 3s sampling ========== 250/1000 | eta: 2s sampling =========== 300/1000 | eta: 2s sampling ============= 350/1000 | eta: 2s sampling =============== 400/1000 | eta: 2s sampling ================= 450/1000 | eta: 2s sampling =================== 500/1000 | eta: 2s sampling ===================== 550/1000 | eta: 1s sampling ======================= 600/1000 | eta: 1s sampling ========================= 650/1000 | eta: 1s sampling =========================== 700/1000 | eta: 1s sampling ============================ 750/1000 | eta: 1s sampling ============================== 800/1000 | eta: 1s sampling ================================ 850/1000 | eta: 0s sampling ================================== 900/1000 | eta: 0s sampling ==================================== 950/1000 | eta: 0s sampling ====================================== 1000/1000 | eta: 0s
# plot credible intervals
coda::gelman.diag(draws, autoburnin = FALSE)
#> Potential scale reduction factors:
#>
#> Point est. Upper C.I.
#> eta[1,1] 1.00 1.00
#> eta[2,1] 1.01 1.03
#> eta[3,1] 1.03 1.08
#> eta[4,1] 1.01 1.04
#> eta[5,1] 1.04 1.10
#> eta[6,1] 1.01 1.02
#> eta[7,1] 1.01 1.04
#> eta[8,1] 1.00 1.00
#> eta[9,1] 1.00 1.01
#> eta[10,1] 1.01 1.02
#>
#> Multivariate psrf
#>
#> 1.11
sims <- calculate(lambda, nsim = 100, values = draws)
dim(sims$lambda)
#> [1] 100 10 1
quants <- apply(sims$lambda[, , 1], 2, quantile, c(0.025, 0.975))
plot(quants[2, ] ~ times,
ylim = c(0, max(c(quants, y))),
type = "l")
lines(quants[1, ] ~ times)
points(y ~ times)
Created on 2024-02-09 with reprex v2.1.0
Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.3.2 (2023-10-31)
#> os macOS Sonoma 14.0
#> system aarch64, darwin20
#> ui X11
#> language (EN)
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Australia/Sydney
#> date 2024-02-09
#> pandoc 3.1.1 @ /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.3.0)
#> backports 1.4.1 2021-12-13 [1] CRAN (R 4.3.0)
#> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.3.0)
#> callr 3.7.3 2022-11-02 [1] CRAN (R 4.3.0)
#> cli 3.6.2 2023-12-11 [1] CRAN (R 4.3.1)
#> coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.3.1)
#> codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.2)
#> crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0)
#> curl 5.2.0 2023-12-08 [1] CRAN (R 4.3.1)
#> digest 0.6.34 2024-01-11 [1] CRAN (R 4.3.1)
#> evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.1)
#> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0)
#> fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.0)
#> future 1.33.1 2023-12-22 [1] CRAN (R 4.3.1)
#> globals 0.16.2 2022-11-21 [1] CRAN (R 4.3.0)
#> glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.1)
#> greta * 0.4.4.9000 2024-02-06 [1] Github (njtierney/greta@db4154e)
#> greta.gam * 0.1.0 2024-02-07 [1] Github (greta-dev/greta.gam@3420beb)
#> greta.gp * 0.2.1 2024-01-29 [1] CRAN (R 4.3.1)
#> highr 0.10 2022-12-22 [1] CRAN (R 4.3.0)
#> hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
#> htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.1)
#> jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.1)
#> knitr 1.45 2023-10-30 [1] CRAN (R 4.3.1)
#> lattice 0.22-5 2023-10-24 [1] CRAN (R 4.3.1)
#> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.1)
#> listenv 0.9.1 2024-01-29 [1] CRAN (R 4.3.1)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0)
#> Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.1)
#> mgcv 1.9-1 2023-12-21 [1] CRAN (R 4.3.1)
#> nlme 3.1-164 2023-11-27 [1] CRAN (R 4.3.1)
#> parallelly 1.36.0 2023-05-26 [1] CRAN (R 4.3.0)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0)
#> png 0.1-8 2022-11-29 [1] CRAN (R 4.3.0)
#> prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.3.1)
#> processx 3.8.3 2023-12-10 [1] CRAN (R 4.3.1)
#> progress 1.2.3 2023-12-06 [1] CRAN (R 4.3.1)
#> ps 1.7.6 2024-01-18 [1] CRAN (R 4.3.1)
#> purrr 1.0.2 2023-08-10 [1] CRAN (R 4.3.0)
#> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.3.0)
#> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.0)
#> R.oo 1.26.0 2024-01-24 [1] CRAN (R 4.3.1)
#> R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.3.1)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0)
#> Rcpp 1.0.12 2024-01-09 [1] CRAN (R 4.3.1)
#> reprex 2.1.0 2024-01-11 [1] CRAN (R 4.3.1)
#> reticulate 1.35.0 2024-01-31 [1] CRAN (R 4.3.1)
#> rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.1)
#> rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1)
#> rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.0)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
#> styler 1.10.2 2023-08-29 [1] CRAN (R 4.3.0)
#> tensorflow 2.15.0 2024-01-31 [1] CRAN (R 4.3.1)
#> tfautograph 0.3.2 2021-09-17 [1] CRAN (R 4.3.0)
#> tfruns 1.5.2 2024-01-26 [1] CRAN (R 4.3.1)
#> vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.1)
#> whisker 0.4.1 2022-12-05 [1] CRAN (R 4.3.0)
#> withr 3.0.0 2024-01-16 [1] CRAN (R 4.3.1)
#> xfun 0.42 2024-02-08 [1] CRAN (R 4.3.1)
#> xml2 1.3.6 2023-12-04 [1] CRAN (R 4.3.1)
#> yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.1)
#>
#> [1] /Library/Frameworks/R.framework/Versions/4.3-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
#>
#> ──────────────────────────────────────────────────────────────────────────────