Skip to content

Instantly share code, notes, and snippets.

@njtierney
Created March 21, 2025 02:14
Show Gist options
  • Save njtierney/7f39e6c7a231b3da23f1708b74a54817 to your computer and use it in GitHub Desktop.
Save njtierney/7f39e6c7a231b3da23f1708b74a54817 to your computer and use it in GitHub Desktop.
# Warming up the sampler and running extra samples on the model with the initial dataset (ie. not changing it each iteration) should tell us whether the sampler is correctly adapted.
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
# example code from test_posteriors_geweke
n <- 10
mu1 <- rnorm(1, 0, 3)
sd1 <- rlnorm(1)
sd2 <- rlnorm(1)

# prior (n draws)
p_theta <- function(n) {
  rnorm(n, mu1, sd1)
}

# likelihood
p_x_bar_theta <- function(theta) {
  rnorm(n, theta, sd2)
}

# define the greta model (single precision for slice sampler)
x <- as_data(rep(0, n))
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#> 
greta_theta <- normal(mu1, sd1)
distribution(x) <- normal(greta_theta, sd2)
model <- model(greta_theta)

draws <- mcmc(
  model,
  warmup = 2000,
  n_samples = 100,
  sampler = adaptive_hmc(),
  chains = 10
)
#> running 10 chains simultaneously on up to 8 CPU cores
#> 
#>   sampling                                            0/100 | eta:  ?s            sampling ===================                       50/100 | eta:  2s            sampling ======================================   100/100 | eta:  0s

extras <- extra_samples(draws = draws, n_samples = 10000)
#> running 10 chains simultaneously on up to 8 CPU cores
#>   sampling                                          0/10000 | eta:  ?s            sampling                                         50/10000 | eta:  7s            sampling                                        100/10000 | eta:  6s            sampling =                                      150/10000 | eta:  6s            sampling =                                      200/10000 | eta:  6s            sampling =                                      250/10000 | eta:  6s            sampling =                                      300/10000 | eta:  6s            sampling =                                      350/10000 | eta:  6s            sampling =                                      400/10000 | eta:  6s            sampling ==                                     450/10000 | eta:  6s            sampling ==                                     500/10000 | eta:  6s            sampling ==                                     550/10000 | eta:  6s            sampling ==                                     600/10000 | eta:  6s            sampling ==                                     650/10000 | eta:  6s            sampling ===                                    700/10000 | eta:  6s            sampling ===                                    750/10000 | eta:  6s            sampling ===                                    800/10000 | eta:  6s            sampling ===                                    850/10000 | eta:  6s            sampling ===                                    900/10000 | eta:  6s            sampling ===                                    950/10000 | eta:  6s            sampling ====                                  1000/10000 | eta:  6s            sampling ====                                  1050/10000 | eta:  6s            sampling ====                                  1100/10000 | eta:  6s            sampling ====                                  1150/10000 | eta:  6s            sampling ====                                  1200/10000 | eta:  6s            sampling ====                                  1250/10000 | eta:  6s            sampling =====                                 1300/10000 | eta:  6s            sampling =====                                 1350/10000 | eta:  6s            sampling =====                                 1400/10000 | eta:  6s            sampling =====                                 1450/10000 | eta:  6s            sampling =====                                 1500/10000 | eta:  6s            sampling ======                                1550/10000 | eta:  5s            sampling ======                                1600/10000 | eta:  5s            sampling ======                                1650/10000 | eta:  5s            sampling ======                                1700/10000 | eta:  5s            sampling ======                                1750/10000 | eta:  5s            sampling ======                                1800/10000 | eta:  5s            sampling =======                               1850/10000 | eta:  5s            sampling =======                               1900/10000 | eta:  5s            sampling =======                               1950/10000 | eta:  5s            sampling =======                               2000/10000 | eta:  5s            sampling =======                               2050/10000 | eta:  5s            sampling ========                              2100/10000 | eta:  5s            sampling ========                              2150/10000 | eta:  5s            sampling ========                              2200/10000 | eta:  5s            sampling ========                              2250/10000 | eta:  5s            sampling ========                              2300/10000 | eta:  5s            sampling ========                              2350/10000 | eta:  5s            sampling =========                             2400/10000 | eta:  5s            sampling =========                             2450/10000 | eta:  5s            sampling =========                             2500/10000 | eta:  5s            sampling =========                             2550/10000 | eta:  5s            sampling =========                             2600/10000 | eta:  5s            sampling ==========                            2650/10000 | eta:  5s            sampling ==========                            2700/10000 | eta:  5s            sampling ==========                            2750/10000 | eta:  5s            sampling ==========                            2800/10000 | eta:  5s            sampling ==========                            2850/10000 | eta:  5s            sampling ==========                            2900/10000 | eta:  5s            sampling ===========                           2950/10000 | eta:  5s            sampling ===========                           3000/10000 | eta:  5s            sampling ===========                           3050/10000 | eta:  5s            sampling ===========                           3100/10000 | eta:  5s            sampling ===========                           3150/10000 | eta:  5s            sampling ============                          3200/10000 | eta:  5s            sampling ============                          3250/10000 | eta:  4s            sampling ============                          3300/10000 | eta:  4s            sampling ============                          3350/10000 | eta:  4s            sampling ============                          3400/10000 | eta:  4s            sampling ============                          3450/10000 | eta:  4s            sampling =============                         3500/10000 | eta:  4s            sampling =============                         3550/10000 | eta:  4s            sampling =============                         3600/10000 | eta:  4s            sampling =============                         3650/10000 | eta:  4s            sampling =============                         3700/10000 | eta:  4s            sampling ==============                        3750/10000 | eta:  4s            sampling ==============                        3800/10000 | eta:  4s            sampling ==============                        3850/10000 | eta:  4s            sampling ==============                        3900/10000 | eta:  4s            sampling ==============                        3950/10000 | eta:  4s            sampling ==============                        4000/10000 | eta:  4s            sampling ===============                       4050/10000 | eta:  4s            sampling ===============                       4100/10000 | eta:  4s            sampling ===============                       4150/10000 | eta:  4s            sampling ===============                       4200/10000 | eta:  4s            sampling ===============                       4250/10000 | eta:  4s            sampling ===============                       4300/10000 | eta:  4s            sampling ================                      4350/10000 | eta:  4s            sampling ================                      4400/10000 | eta:  4s            sampling ================                      4450/10000 | eta:  4s            sampling ================                      4500/10000 | eta:  4s            sampling ================                      4550/10000 | eta:  4s            sampling =================                     4600/10000 | eta:  4s            sampling =================                     4650/10000 | eta:  4s            sampling =================                     4700/10000 | eta:  3s            sampling =================                     4750/10000 | eta:  3s            sampling =================                     4800/10000 | eta:  3s            sampling =================                     4850/10000 | eta:  3s            sampling ==================                    4900/10000 | eta:  3s            sampling ==================                    4950/10000 | eta:  3s            sampling ==================                    5000/10000 | eta:  3s            sampling ==================                    5050/10000 | eta:  3s            sampling ==================                    5100/10000 | eta:  3s            sampling ===================                   5150/10000 | eta:  3s            sampling ===================                   5200/10000 | eta:  3s            sampling ===================                   5250/10000 | eta:  3s            sampling ===================                   5300/10000 | eta:  3s            sampling ===================                   5350/10000 | eta:  3s            sampling ===================                   5400/10000 | eta:  3s            sampling ====================                  5450/10000 | eta:  3s            sampling ====================                  5500/10000 | eta:  3s            sampling ====================                  5550/10000 | eta:  3s            sampling ====================                  5600/10000 | eta:  3s            sampling ====================                  5650/10000 | eta:  3s            sampling =====================                 5700/10000 | eta:  3s            sampling =====================                 5750/10000 | eta:  3s            sampling =====================                 5800/10000 | eta:  3s            sampling =====================                 5850/10000 | eta:  3s            sampling =====================                 5900/10000 | eta:  3s            sampling =====================                 5950/10000 | eta:  3s            sampling ======================                6000/10000 | eta:  3s            sampling ======================                6050/10000 | eta:  3s            sampling ======================                6100/10000 | eta:  3s            sampling ======================                6150/10000 | eta:  2s            sampling ======================                6200/10000 | eta:  2s            sampling ======================                6250/10000 | eta:  2s            sampling =======================               6300/10000 | eta:  2s            sampling =======================               6350/10000 | eta:  2s            sampling =======================               6400/10000 | eta:  2s            sampling =======================               6450/10000 | eta:  2s            sampling =======================               6500/10000 | eta:  2s            sampling ========================              6550/10000 | eta:  2s            sampling ========================              6600/10000 | eta:  2s            sampling ========================              6650/10000 | eta:  2s            sampling ========================              6700/10000 | eta:  2s            sampling ========================              6750/10000 | eta:  2s            sampling ========================              6800/10000 | eta:  2s            sampling =========================             6850/10000 | eta:  2s            sampling =========================             6900/10000 | eta:  2s            sampling =========================             6950/10000 | eta:  2s            sampling =========================             7000/10000 | eta:  2s            sampling =========================             7050/10000 | eta:  2s            sampling ==========================            7100/10000 | eta:  2s            sampling ==========================            7150/10000 | eta:  2s            sampling ==========================            7200/10000 | eta:  2s            sampling ==========================            7250/10000 | eta:  2s            sampling ==========================            7300/10000 | eta:  2s            sampling ==========================            7350/10000 | eta:  2s            sampling ===========================           7400/10000 | eta:  2s            sampling ===========================           7450/10000 | eta:  2s            sampling ===========================           7500/10000 | eta:  2s            sampling ===========================           7550/10000 | eta:  2s            sampling ===========================           7600/10000 | eta:  2s            sampling ============================          7650/10000 | eta:  1s            sampling ============================          7700/10000 | eta:  1s            sampling ============================          7750/10000 | eta:  1s            sampling ============================          7800/10000 | eta:  1s            sampling ============================          7850/10000 | eta:  1s            sampling ============================          7900/10000 | eta:  1s            sampling =============================         7950/10000 | eta:  1s            sampling =============================         8000/10000 | eta:  1s            sampling =============================         8050/10000 | eta:  1s            sampling =============================         8100/10000 | eta:  1s            sampling =============================         8150/10000 | eta:  1s            sampling ==============================        8200/10000 | eta:  1s            sampling ==============================        8250/10000 | eta:  1s            sampling ==============================        8300/10000 | eta:  1s            sampling ==============================        8350/10000 | eta:  1s            sampling ==============================        8400/10000 | eta:  1s            sampling ==============================        8450/10000 | eta:  1s            sampling ===============================       8500/10000 | eta:  1s            sampling ===============================       8550/10000 | eta:  1s            sampling ===============================       8600/10000 | eta:  1s            sampling ===============================       8650/10000 | eta:  1s            sampling ===============================       8700/10000 | eta:  1s            sampling ================================      8750/10000 | eta:  1s            sampling ================================      8800/10000 | eta:  1s            sampling ================================      8850/10000 | eta:  1s            sampling ================================      8900/10000 | eta:  1s            sampling ================================      8950/10000 | eta:  1s            sampling ================================      9000/10000 | eta:  1s            sampling =================================     9050/10000 | eta:  1s            sampling =================================     9100/10000 | eta:  1s            sampling =================================     9150/10000 | eta:  1s            sampling =================================     9200/10000 | eta:  1s            sampling =================================     9250/10000 | eta:  0s            sampling =================================     9300/10000 | eta:  0s            sampling ==================================    9350/10000 | eta:  0s            sampling ==================================    9400/10000 | eta:  0s            sampling ==================================    9450/10000 | eta:  0s            sampling ==================================    9500/10000 | eta:  0s            sampling ==================================    9550/10000 | eta:  0s            sampling ===================================   9600/10000 | eta:  0s            sampling ===================================   9650/10000 | eta:  0s            sampling ===================================   9700/10000 | eta:  0s            sampling ===================================   9750/10000 | eta:  0s            sampling ===================================   9800/10000 | eta:  0s            sampling ===================================   9850/10000 | eta:  0s            sampling ====================================  9900/10000 | eta:  0s            sampling ====================================  9950/10000 | eta:  0s            sampling ==================================== 10000/10000 | eta:  0s

plot(extras)

Created on 2025-03-21 with reprex v2.1.1

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.4.2 (2024-10-31)
#>  os       macOS Sequoia 15.3.2
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Australia/Hobart
#>  date     2025-03-21
#>  pandoc   3.2.1 @ /opt/homebrew/bin/ (via rmarkdown)
#>  quarto   1.6.36 @ /usr/local/bin/quarto
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date (UTC) lib source
#>  abind         1.4-8      2024-09-12 [1] CRAN (R 4.4.1)
#>  backports     1.5.0      2024-05-23 [1] CRAN (R 4.4.0)
#>  base64enc     0.1-3      2015-07-28 [1] CRAN (R 4.4.0)
#>  callr         3.7.6      2024-03-25 [1] CRAN (R 4.4.0)
#>  cli           3.6.4      2025-02-13 [1] CRAN (R 4.4.1)
#>  coda          0.19-4.1   2024-01-31 [1] CRAN (R 4.4.0)
#>  codetools     0.2-20     2024-03-31 [2] CRAN (R 4.4.2)
#>  crayon        1.5.3      2024-06-20 [1] CRAN (R 4.4.0)
#>  curl          6.2.1      2025-02-19 [1] CRAN (R 4.4.1)
#>  digest        0.6.37     2024-08-19 [1] CRAN (R 4.4.1)
#>  evaluate      1.0.3      2025-01-10 [1] CRAN (R 4.4.1)
#>  fastmap       1.2.0      2024-05-15 [1] CRAN (R 4.4.0)
#>  fs            1.6.5      2024-10-30 [1] CRAN (R 4.4.1)
#>  future        1.34.0     2024-07-29 [1] CRAN (R 4.4.0)
#>  globals       0.16.3     2024-03-08 [1] CRAN (R 4.4.0)
#>  glue          1.8.0      2024-09-30 [1] CRAN (R 4.4.1)
#>  greta       * 0.5.0.9000 2025-03-21 [1] local
#>  hms           1.1.3      2023-03-21 [1] CRAN (R 4.4.0)
#>  htmltools     0.5.8.1    2024-04-04 [1] CRAN (R 4.4.0)
#>  jsonlite      1.9.1      2025-03-03 [1] CRAN (R 4.4.1)
#>  knitr         1.50       2025-03-16 [1] CRAN (R 4.4.1)
#>  lattice       0.22-6     2024-03-20 [2] CRAN (R 4.4.2)
#>  lifecycle     1.0.4      2023-11-07 [1] CRAN (R 4.4.0)
#>  listenv       0.9.1      2024-01-29 [1] CRAN (R 4.4.0)
#>  magrittr      2.0.3      2022-03-30 [1] CRAN (R 4.4.0)
#>  Matrix        1.7-3      2025-03-11 [2] CRAN (R 4.4.1)
#>  parallelly    1.42.0     2025-01-30 [1] CRAN (R 4.4.1)
#>  pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.4.0)
#>  png           0.1-8      2022-11-29 [1] CRAN (R 4.4.0)
#>  prettyunits   1.2.0      2023-09-24 [1] CRAN (R 4.4.0)
#>  processx      3.8.6      2025-02-21 [1] CRAN (R 4.4.1)
#>  progress      1.2.3      2023-12-06 [1] CRAN (R 4.4.0)
#>  ps            1.9.0      2025-02-18 [1] CRAN (R 4.4.1)
#>  R6            2.6.1      2025-02-15 [1] CRAN (R 4.4.1)
#>  Rcpp          1.0.14     2025-01-12 [1] CRAN (R 4.4.1)
#>  reprex        2.1.1      2024-07-06 [1] CRAN (R 4.4.0)
#>  reticulate    1.41.0.1   2025-03-09 [1] CRAN (R 4.4.1)
#>  rlang         1.1.5      2025-01-17 [1] CRAN (R 4.4.1)
#>  rmarkdown     2.29       2024-11-04 [1] CRAN (R 4.4.1)
#>  rstudioapi    0.17.1     2024-10-22 [1] CRAN (R 4.4.1)
#>  sessioninfo   1.2.3      2025-02-05 [1] CRAN (R 4.4.1)
#>  tensorflow    2.16.0     2024-04-15 [1] CRAN (R 4.4.0)
#>  tfautograph   0.3.2      2021-09-17 [1] CRAN (R 4.4.0)
#>  tfruns        1.5.3      2024-04-19 [1] CRAN (R 4.4.0)
#>  vctrs         0.6.5      2023-12-01 [1] CRAN (R 4.4.0)
#>  whisker       0.4.1      2022-12-05 [1] CRAN (R 4.4.0)
#>  withr         3.0.2      2024-10-28 [1] CRAN (R 4.4.1)
#>  xfun          0.51       2025-02-19 [1] CRAN (R 4.4.1)
#>  xml2          1.3.8      2025-03-14 [1] CRAN (R 4.4.1)
#>  yaml          2.3.10     2024-07-26 [1] CRAN (R 4.4.0)
#> 
#>  [1] /Users/nick/Library/R/arm64/4.4/library
#>  [2] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
#>  * ── Packages attached to the search path.
#> 
#> ─ 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.10.dylib
#>  pythonhome:     /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2:/Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2
#>  version:        3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:51:49) [Clang 16.0.6 ]
#>  numpy:          /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/python3.10/site-packages/numpy
#>  numpy_version:  1.26.4
#>  tensorflow:     /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/python3.10/site-packages/tensorflow
#>  
#>  NOTE: Python version was forced by use_python() function
#> 
#> ──────────────────────────────────────────────────────────────────────────────
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment