greta-dev / greta

simple and scalable statistical modelling in R
https://greta-stats.org
Other
527 stars 63 forks source link

Cholesky factor should be upper triangular to match R's behaviour #692

Open njtierney opened 4 weeks ago

njtierney commented 4 weeks ago
# In R a cholesky factor is upper triangular (zeros on the lower triangle)
x <- rWishart(1, 4, diag(3))[, , 1]
chol(x)
#>         [,1]     [,2]       [,3]
#> [1,] 2.15256 1.348369 -0.5904764
#> [2,] 0.00000 1.486747 -0.5220772
#> [3,] 0.00000 0.000000  2.0428290

# in greta, we want to mimic this:
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
xg <- wishart(df = 4, Sigma = diag(3))
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#> 
# zeroes are already there?
xg
#> greta array <variable following a wishart distribution>
#>      [,1] [,2] [,3]
#> [1,]  ?    ?    ?  
#> [2,] 0     ?    ?  
#> [3,] 0    0     ?
#> 
chol_xg <- chol(xg)
# cholesky factor is lower triangle -  the zeroes are in the upper triangle
calculate(chol_xg, nsim = 1)
#> $chol_xg
#> , , 1
#> 
#>          [,1] [,2] [,3]
#> [1,] 1.848067    0    0
#> 
#> , , 2
#> 
#>           [,1]     [,2] [,3]
#> [1,] 0.6420776 1.114142    0
#> 
#> , , 3
#> 
#>          [,1]       [,2]     [,3]
#> [1,] 3.021485 -0.1259585 1.464119

Created on 2024-08-13 with reprex v2.1.1

Session info ``` r 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-08-13 #> pandoc 3.2.1 @ /opt/homebrew/bin/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> 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.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) #> digest 0.6.36 2024-06-23 [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) #> greta * 0.4.5.9000 2024-08-01 [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.8.8 2023-12-04 [1] CRAN (R 4.4.0) #> knitr 1.48 2024-07-07 [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) #> parallelly 1.38.0 2024-07-27 [1] CRAN (R 4.4.0) #> 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.4 2024-03-16 [1] CRAN (R 4.4.0) #> progress 1.2.3 2023-12-06 [1] CRAN (R 4.4.0) #> ps 1.7.7 2024-07-02 [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) #> 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) #> 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) #> 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.1 2024-07-31 [1] CRAN (R 4.4.0) #> xfun 0.46 2024-07-18 [1] CRAN (R 4.4.0) #> 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 #> #> ─ 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 #> #> ────────────────────────────────────────────────────────────────────────────── ```