Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bias in sampling LKJ distribution #733

Open
njtierney opened this issue Oct 18, 2024 · 0 comments
Open

Bias in sampling LKJ distribution #733

njtierney opened this issue Oct 18, 2024 · 0 comments

Comments

@njtierney
Copy link
Collaborator

njtierney commented Oct 18, 2024 •
edited
Loading 

devtools::load_all(".")
#> ℹ Loading greta
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> 
#> ✔ Initialising python and checking dependencies ... done!

eta <- 3
dim <- 2
x <- lkj_correlation(eta = eta, dimension = dim)[1, 2]
iid_lkj <- function(n) {
  rlkjcorr(n, eta = eta, dimension = dim)[, 1, 2]
}

m <- model(x)
n_draws <- 5000
draws <- mcmc(m, warmup = 2000, n_samples = n_draws, one_by_one = TRUE)
#> running 4 chains simultaneously on up to 8 CPU cores
#>     warmup                                           0/2000 | eta:  ?s              warmup =                                        50/2000 | eta: 49s              warmup ==                                      100/2000 | eta: 27s              warmup ===                                     150/2000 | eta: 20s              warmup ====                                    200/2000 | eta: 16s              warmup =====                                   250/2000 | eta: 14s              warmup ======                                  300/2000 | eta: 12s              warmup =======                                 350/2000 | eta: 11s              warmup ========                                400/2000 | eta: 10s              warmup =========                               450/2000 | eta:  9s              warmup ==========                              500/2000 | eta:  9s              warmup ==========                              550/2000 | eta:  8s              warmup ===========                             600/2000 | eta:  8s              warmup ============                            650/2000 | eta:  7s              warmup =============                           700/2000 | eta:  7s              warmup ==============                          750/2000 | eta:  6s              warmup ===============                         800/2000 | eta:  6s              warmup ================                        850/2000 | eta:  6s              warmup =================                       900/2000 | eta:  6s              warmup ==================                      950/2000 | eta:  5s              warmup ===================                    1000/2000 | eta:  5s              warmup ====================                   1050/2000 | eta:  5s              warmup =====================                  1100/2000 | eta:  4s              warmup ======================                 1150/2000 | eta:  4s              warmup =======================                1200/2000 | eta:  4s              warmup ========================               1250/2000 | eta:  4s              warmup =========================              1300/2000 | eta:  3s              warmup ==========================             1350/2000 | eta:  3s              warmup ===========================            1400/2000 | eta:  3s              warmup ============================           1450/2000 | eta:  3s              warmup ============================           1500/2000 | eta:  2s              warmup =============================          1550/2000 | eta:  2s              warmup ==============================         1600/2000 | eta:  2s              warmup ===============================        1650/2000 | eta:  2s              warmup ================================       1700/2000 | eta:  1s              warmup =================================      1750/2000 | eta:  1s              warmup ==================================     1800/2000 | eta:  1s              warmup ===================================    1850/2000 | eta:  1s              warmup ====================================   1900/2000 | eta:  0s              warmup =====================================  1950/2000 | eta:  0s              warmup ====================================== 2000/2000 | eta:  0s          
#>   sampling                                           0/5000 | eta:  ?s            sampling                                          50/5000 | eta: 17s            sampling =                                       100/5000 | eta: 18s            sampling =                                       150/5000 | eta: 17s            sampling ==                                      200/5000 | eta: 17s            sampling ==                                      250/5000 | eta: 17s            sampling ==                                      300/5000 | eta: 17s            sampling ===                                     350/5000 | eta: 17s            sampling ===                                     400/5000 | eta: 16s            sampling ===                                     450/5000 | eta: 16s            sampling ====                                    500/5000 | eta: 16s            sampling ====                                    550/5000 | eta: 16s            sampling =====                                   600/5000 | eta: 16s            sampling =====                                   650/5000 | eta: 16s            sampling =====                                   700/5000 | eta: 15s            sampling ======                                  750/5000 | eta: 15s            sampling ======                                  800/5000 | eta: 15s            sampling ======                                  850/5000 | eta: 15s            sampling =======                                 900/5000 | eta: 15s            sampling =======                                 950/5000 | eta: 14s            sampling ========                               1000/5000 | eta: 14s            sampling ========                               1050/5000 | eta: 14s            sampling ========                               1100/5000 | eta: 14s            sampling =========                              1150/5000 | eta: 14s            sampling =========                              1200/5000 | eta: 14s            sampling ==========                             1250/5000 | eta: 13s            sampling ==========                             1300/5000 | eta: 13s            sampling ==========                             1350/5000 | eta: 13s            sampling ===========                            1400/5000 | eta: 13s            sampling ===========                            1450/5000 | eta: 13s            sampling ===========                            1500/5000 | eta: 13s            sampling ============                           1550/5000 | eta: 12s            sampling ============                           1600/5000 | eta: 12s            sampling =============                          1650/5000 | eta: 12s            sampling =============                          1700/5000 | eta: 12s            sampling =============                          1750/5000 | eta: 12s            sampling ==============                         1800/5000 | eta: 11s            sampling ==============                         1850/5000 | eta: 11s            sampling ==============                         1900/5000 | eta: 11s            sampling ===============                        1950/5000 | eta: 11s            sampling ===============                        2000/5000 | eta: 11s            sampling ================                       2050/5000 | eta: 11s            sampling ================                       2100/5000 | eta: 11s            sampling ================                       2150/5000 | eta: 10s            sampling =================                      2200/5000 | eta: 10s            sampling =================                      2250/5000 | eta: 10s            sampling =================                      2300/5000 | eta: 10s            sampling ==================                     2350/5000 | eta: 10s            sampling ==================                     2400/5000 | eta: 10s            sampling ===================                    2450/5000 | eta:  9s            sampling ===================                    2500/5000 | eta:  9s            sampling ===================                    2550/5000 | eta:  9s            sampling ====================                   2600/5000 | eta:  9s            sampling ====================                   2650/5000 | eta:  9s            sampling =====================                  2700/5000 | eta:  8s            sampling =====================                  2750/5000 | eta:  8s            sampling =====================                  2800/5000 | eta:  8s            sampling ======================                 2850/5000 | eta:  8s            sampling ======================                 2900/5000 | eta:  8s            sampling ======================                 2950/5000 | eta:  8s            sampling =======================                3000/5000 | eta:  7s            sampling =======================                3050/5000 | eta:  7s            sampling ========================               3100/5000 | eta:  7s            sampling ========================               3150/5000 | eta:  7s            sampling ========================               3200/5000 | eta:  7s            sampling =========================              3250/5000 | eta:  7s            sampling =========================              3300/5000 | eta:  6s            sampling =========================              3350/5000 | eta:  6s            sampling ==========================             3400/5000 | eta:  6s            sampling ==========================             3450/5000 | eta:  6s            sampling ===========================            3500/5000 | eta:  6s            sampling ===========================            3550/5000 | eta:  5s            sampling ===========================            3600/5000 | eta:  5s            sampling ============================           3650/5000 | eta:  5s            sampling ============================           3700/5000 | eta:  5s            sampling ============================           3750/5000 | eta:  5s            sampling =============================          3800/5000 | eta:  5s            sampling =============================          3850/5000 | eta:  4s            sampling ==============================         3900/5000 | eta:  4s            sampling ==============================         3950/5000 | eta:  4s            sampling ==============================         4000/5000 | eta:  4s            sampling ===============================        4050/5000 | eta:  4s            sampling ===============================        4100/5000 | eta:  3s            sampling ================================       4150/5000 | eta:  3s            sampling ================================       4200/5000 | eta:  3s            sampling ================================       4250/5000 | eta:  3s            sampling =================================      4300/5000 | eta:  3s            sampling =================================      4350/5000 | eta:  2s            sampling =================================      4400/5000 | eta:  2s            sampling ==================================     4450/5000 | eta:  2s            sampling ==================================     4500/5000 | eta:  2s            sampling ===================================    4550/5000 | eta:  2s            sampling ===================================    4600/5000 | eta:  2s            sampling ===================================    4650/5000 | eta:  1s            sampling ====================================   4700/5000 | eta:  1s            sampling ====================================   4750/5000 | eta:  1s            sampling ====================================   4800/5000 | eta:  1s            sampling =====================================  4850/5000 | eta:  1s            sampling =====================================  4900/5000 | eta:  0s            sampling ====================================== 4950/5000 | eta:  0s            sampling ====================================== 5000/5000 | eta:  0s
iid_draws <- iid_lkj(n = n_draws)
plot(draws)

coda::gelman.diag(draws)
#> Potential scale reduction factors:
#> 
#>   Point est. Upper C.I.
#> x          1          1

stats::qqplot(
  x = as.matrix(draws),
  y = iid_draws,
  main = "lkj draws with HMC"
)

graphics::abline(0, 1)

Created on 2024-10-18 with reprex v2.1.1

Session info 
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.4.1 Patched (2024-07-08 r86915)
#>  os       macOS Sonoma 14.5
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Australia/Hobart
#>  date     2024-10-18
#>  pandoc   3.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  ! package     * version    date (UTC) lib source
#>    abind         1.4-5      2016-07-21 [1] CRAN (R 4.4.0)
#>    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)
#>    brio          1.1.5      2024-04-24 [1] CRAN (R 4.4.0)
#>    cachem        1.1.0      2024-05-16 [1] CRAN (R 4.4.0)
#>    callr         3.7.6      2024-03-25 [1] CRAN (R 4.4.0)
#>    cli           3.6.3      2024-06-21 [1] CRAN (R 4.4.0)
#>    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.1)
#>    crayon        1.5.3      2024-06-20 [1] CRAN (R 4.4.0)
#>    curl          5.2.1      2024-03-01 [1] CRAN (R 4.4.0)
#>    desc          1.4.3      2023-12-10 [1] CRAN (R 4.4.0)
#>    devtools      2.4.5      2022-10-11 [1] CRAN (R 4.4.0)
#>    digest        0.6.36     2024-06-23 [1] CRAN (R 4.4.0)
#>    ellipsis      0.3.2      2021-04-29 [1] CRAN (R 4.4.0)
#>    evaluate      0.24.0     2024-06-10 [1] CRAN (R 4.4.0)
#>    fastmap       1.2.0      2024-05-15 [1] CRAN (R 4.4.0)
#>    fs            1.6.4.9000 2024-06-26 [1] Github (r-lib/fs@714990b)
#>    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.7.0      2024-01-09 [1] CRAN (R 4.4.0)
#>  P greta       * 0.5.0      2024-10-16 [?] load_all()
#>    highr         0.11       2024-05-26 [1] CRAN (R 4.4.0)
#>    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)
#>    htmlwidgets   1.6.4      2023-12-06 [1] CRAN (R 4.4.0)
#>    httpuv        1.6.15     2024-03-26 [1] CRAN (R 4.4.0)
#>    jsonlite      1.8.8      2023-12-04 [1] CRAN (R 4.4.0)
#>    knitr         1.48       2024-07-07 [1] CRAN (R 4.4.0)
#>    later         1.3.2      2023-12-06 [1] CRAN (R 4.4.0)
#>    lattice       0.22-6     2024-03-20 [2] CRAN (R 4.4.1)
#>    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-0      2024-04-26 [2] CRAN (R 4.4.1)
#>    memoise       2.0.1.9000 2024-08-14 [1] Github (hadley/memoise@40db995)
#>    mime          0.12       2021-09-28 [1] CRAN (R 4.4.0)
#>    miniUI        0.1.1.1    2018-05-18 [1] CRAN (R 4.4.0)
#>    parallelly    1.38.0     2024-07-27 [1] CRAN (R 4.4.0)
#>    pkgbuild      1.4.4      2024-03-17 [1] CRAN (R 4.4.0)
#>    pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.4.0)
#>    pkgload       1.4.0      2024-06-28 [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.4      2024-03-16 [1] CRAN (R 4.4.0)
#>    profvis       0.3.8      2023-05-02 [1] CRAN (R 4.4.0)
#>    progress      1.2.3      2023-12-06 [1] CRAN (R 4.4.0)
#>    promises      1.3.0      2024-04-05 [1] CRAN (R 4.4.0)
#>    ps            1.8.0      2024-09-12 [1] CRAN (R 4.4.1)
#>    purrr         1.0.2      2023-08-10 [1] CRAN (R 4.4.0)
#>    R6            2.5.1      2021-08-19 [1] CRAN (R 4.4.0)
#>    Rcpp          1.0.13     2024-07-17 [1] CRAN (R 4.4.0)
#>    remotes       2.5.0      2024-03-17 [1] CRAN (R 4.4.0)
#>    reprex        2.1.1      2024-07-06 [1] CRAN (R 4.4.0)
#>    reticulate    1.38.0     2024-06-19 [1] CRAN (R 4.4.0)
#>    rlang         1.1.4      2024-06-04 [1] CRAN (R 4.4.0)
#>    rmarkdown     2.27       2024-05-17 [1] CRAN (R 4.4.0)
#>    rprojroot     2.0.4      2023-11-05 [1] CRAN (R 4.4.0)
#>    rstudioapi    0.16.0     2024-03-24 [1] CRAN (R 4.4.0)
#>    sessioninfo   1.2.2      2021-12-06 [1] CRAN (R 4.4.0)
#>    shiny         1.9.1      2024-08-01 [1] CRAN (R 4.4.0)
#>    stringi       1.8.4      2024-05-06 [1] CRAN (R 4.4.0)
#>    stringr       1.5.1      2023-11-14 [1] CRAN (R 4.4.0)
#>    tensorflow    2.16.0     2024-04-15 [1] CRAN (R 4.4.0)
#>    testthat    * 3.2.1.1    2024-04-14 [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)
#>    urlchecker    1.0.1.9000 2024-08-27 [1] Github (r-lib/urlchecker@ac38ea4)
#>    usethis       3.0.0      2024-07-29 [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.1      2024-07-31 [1] CRAN (R 4.4.0)
#>    xfun          0.46       2024-07-18 [1] CRAN (R 4.4.0)
#>    xml2          1.3.6      2023-12-04 [1] CRAN (R 4.4.0)
#>    xtable        1.8-4      2019-04-21 [1] CRAN (R 4.4.0)
#>    yaml          2.3.10     2024-07-26 [1] CRAN (R 4.4.0)
#>    yesno         0.1.3      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
#> 
#>  P ── Loaded and on-disk path mismatch.
#> 
#> ─ 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
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Related to #729
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
{greta} CRAN 0.5.0 release
Status: Weekly Todo
Milestone
No milestone
Development

No branches or pull requests
1 participant
@njtierney