Outer() converting greta array values to NA (or: how to model multiple time series in greta)

BobOHara · December 8, 2022, 3:30pm

The short version:

I want to use a variable in several parts of a model, and for reasons explained below I want to use outer() to copy the variable over. But it converts the elements of the array to NA.
x = exponential(1, dim=2)
x_arr <- ones(3,3) %o% x
I also tried sweep(), but the greta version doesn’t like 3D arrays.

Does anyone have any suggestions?

The long version:

I’m trying to model multile time series (the full version is more complicated). I got it working badly, and have been trying to follow Nick’s advice on speeding up time series (Simulations with greta? (including time-series)). I got as far as the code below, but the %o% and outer() create the NAs, so the code won’t work.

Simulate some data

NSeries <- 3
NTime <- 5
Abund <- matrix(rpois(NSeries*NTime, 10), nrow=NSeries)

Hyperparameters

lnrho = normal(0, sd=0.5, dim=NSeries)
mu0 = normal(0, sd=0.5, dim=NSeries)
sigma = exponential(1, dim=NSeries)

Matrix of time modifiers

t_seq <- seq_len(NTime)
t_mat <- outer(t_seq, t_seq, FUN = “-”)
t_mat <- pmax(t_mat, 0)
mask <- lower.tri(t_mat, diag = TRUE)

matrix of rho ^ n contributions

lnrho_arr <- outer(t_mat, lnrho, FUN="*")
mask_arr <- mask%o%ones(NSeries)

rho_arr <- exp(lnrho_arr) * mask_arr # positive rho is OK for this problem

Shocks

sigma_arr <- ones(NTime,NTime) %o% sigma
innovations = normal(0, 1, dim = c(N,N,NSeries,1))*sigma_arr

Finally to (not) model the data

loglam <- rho_arr*innovations # This is wrong (e.g. no mu0), but can be fixed later
distribution(Abund) = poisson(exp(loglam), dim=c(NSeries, NTime))

njtierney · December 12, 2022, 12:10am

Hi @BobOHara!

Thanks for posting this - I’ve just tried to reproduce your example on the dev greta branch, but it looks like the object N needs to be defined?

Here’s my reprex below:

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
# Simulate some data
NSeries <- 3
NTime <- 5
Abund <- matrix(rpois(NSeries * NTime, 10), nrow = NSeries)

# Hyperparameters
lnrho <- normal(0, sd = 0.5, dim = NSeries)
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#> 
mu0 <- normal(0, sd = 0.5, dim = NSeries)
sigma <- exponential(1, dim = NSeries)

# Matrix of time modifiers
t_seq <- seq_len(NTime)
t_mat <- outer(t_seq, t_seq, FUN = "-")
t_mat <- pmax(t_mat, 0)
mask <- lower.tri(t_mat, diag = TRUE)

# matrix of rho ^ n contributions
lnrho_arr <- outer(t_mat, lnrho, FUN = "*")
mask_arr <- mask %o% ones(NSeries)

rho_arr <- exp(lnrho_arr) * mask_arr # positive rho is OK for this problem

# Shocks
sigma_arr <- ones(NTime, NTime) %o% sigma
innovations <- normal(0, 1, dim = c(N, N, NSeries, 1)) * sigma_arr
#> Error in check_dims(mean, sd, target_dim = dim): object 'N' not found

# Finally to (not) model the data
loglam <- rho_arr * innovations # This is wrong (e.g. no mu0), but can be fixed later
#> Error in eval(expr, envir, enclos): object 'innovations' not found
distribution(Abund) <- poisson(exp(loglam), dim = c(NSeries, NTime))
#> Error in check_in_family("poisson", lambda): object 'loglam' not found

^{Created on 2022-12-12 with reprex v2.0.2}

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.2.1 (2022-06-23)
#>  os       macOS Monterey 12.3.1
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Australia/Brisbane
#>  date     2022-12-12
#>  pandoc   2.19.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date (UTC) lib source
#>  base64enc     0.1-3      2015-07-28 [1] CRAN (R 4.2.0)
#>  callr         3.7.3      2022-11-02 [1] CRAN (R 4.2.0)
#>  cli           3.4.1      2022-09-23 [1] CRAN (R 4.2.0)
#>  coda          0.19-4     2020-09-30 [1] CRAN (R 4.2.0)
#>  codetools     0.2-18     2020-11-04 [1] CRAN (R 4.2.1)
#>  crayon        1.5.2      2022-09-29 [1] CRAN (R 4.2.0)
#>  digest        0.6.30     2022-10-18 [1] CRAN (R 4.2.0)
#>  ellipsis      0.3.2      2021-04-29 [1] CRAN (R 4.2.0)
#>  evaluate      0.18       2022-11-07 [1] CRAN (R 4.2.0)
#>  fastmap       1.1.0      2021-01-25 [1] CRAN (R 4.2.0)
#>  fs            1.5.2      2021-12-08 [1] CRAN (R 4.2.0)
#>  future        1.29.0     2022-11-06 [1] CRAN (R 4.2.0)
#>  globals       0.16.2     2022-11-21 [1] CRAN (R 4.2.1)
#>  glue          1.6.2      2022-02-24 [1] CRAN (R 4.2.0)
#>  greta       * 0.4.2.9000 2022-12-12 [1] local
#>  here          1.0.1      2020-12-13 [1] CRAN (R 4.2.0)
#>  highr         0.9        2021-04-16 [1] CRAN (R 4.2.0)
#>  hms           1.1.2      2022-08-19 [1] CRAN (R 4.2.0)
#>  htmltools     0.5.3      2022-07-18 [1] CRAN (R 4.2.0)
#>  jsonlite      1.8.3      2022-10-21 [1] CRAN (R 4.2.0)
#>  knitr         1.41       2022-11-18 [1] CRAN (R 4.2.0)
#>  lattice       0.20-45    2021-09-22 [1] CRAN (R 4.2.1)
#>  lifecycle     1.0.3      2022-10-07 [1] CRAN (R 4.2.0)
#>  listenv       0.8.0      2019-12-05 [1] CRAN (R 4.2.0)
#>  magrittr      2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
#>  Matrix        1.5-3      2022-11-11 [1] CRAN (R 4.2.0)
#>  parallelly    1.32.1     2022-07-21 [1] CRAN (R 4.2.0)
#>  pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.2.0)
#>  png           0.1-8      2022-11-29 [1] CRAN (R 4.2.0)
#>  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.2.0)
#>  processx      3.8.0      2022-10-26 [1] CRAN (R 4.2.0)
#>  progress      1.2.2      2019-05-16 [1] CRAN (R 4.2.0)
#>  ps            1.7.2      2022-10-26 [1] CRAN (R 4.2.0)
#>  purrr         0.3.5      2022-10-06 [1] CRAN (R 4.2.0)
#>  R.cache       0.16.0     2022-07-21 [1] CRAN (R 4.2.0)
#>  R.methodsS3   1.8.2      2022-06-13 [1] CRAN (R 4.2.0)
#>  R.oo          1.25.0     2022-06-12 [1] CRAN (R 4.2.0)
#>  R.utils       2.12.2     2022-11-11 [1] CRAN (R 4.2.0)
#>  R6            2.5.1      2021-08-19 [1] CRAN (R 4.2.0)
#>  Rcpp          1.0.9      2022-07-08 [1] CRAN (R 4.2.0)
#>  reprex        2.0.2      2022-08-17 [1] CRAN (R 4.2.0)
#>  reticulate    1.26       2022-08-31 [1] CRAN (R 4.2.0)
#>  rlang         1.0.6      2022-09-24 [1] CRAN (R 4.2.0)
#>  rmarkdown     2.18       2022-11-09 [1] CRAN (R 4.2.0)
#>  rprojroot     2.0.3      2022-04-02 [1] CRAN (R 4.2.0)
#>  rstudioapi    0.14       2022-08-22 [1] CRAN (R 4.2.0)
#>  sessioninfo   1.2.2      2021-12-06 [1] CRAN (R 4.2.0)
#>  stringi       1.7.8      2022-07-11 [1] CRAN (R 4.2.0)
#>  stringr       1.5.0      2022-12-02 [1] CRAN (R 4.2.0)
#>  styler        1.8.1      2022-11-07 [1] CRAN (R 4.2.0)
#>  tensorflow    2.9.0      2022-05-21 [1] CRAN (R 4.2.0)
#>  tfruns        1.5.1      2022-09-05 [1] CRAN (R 4.2.0)
#>  vctrs         0.5.1      2022-11-16 [1] CRAN (R 4.2.0)
#>  whisker       0.4        2019-08-28 [1] CRAN (R 4.2.0)
#>  withr         2.5.0      2022-03-03 [1] CRAN (R 4.2.0)
#>  xfun          0.35       2022-11-16 [1] CRAN (R 4.2.0)
#>  yaml          2.3.6      2022-10-18 [1] CRAN (R 4.2.0)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.2-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.13 | packaged by conda-forge | (default, Mar 25 2022, 06:05:16)  [Clang 12.0.1 ]
#>  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
#> 
#> ──────────────────────────────────────────────────────────────────────────────

BobOHara · December 12, 2022, 2:18pm

Apologies, N should just be NTime.

njtierney · December 12, 2022, 11:37pm

All good - I should have checked the initial outer code! It returns NA values as you said, here’s the reprex:

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
NSeries <- 3
NTime <- 5
Abund <- matrix(rpois(NSeries * NTime, 10), nrow = NSeries)

# Hyperparameters
lnrho <- normal(0, sd = 0.5, dim = NSeries)
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#> 

mu0 <- normal(0, sd = 0.5, dim = NSeries)
sigma <- exponential(1, dim = NSeries)

# Matrix of time modifiers
t_seq <- seq_len(NTime)
t_mat <- outer(t_seq, t_seq, FUN = "-")
t_mat <- pmax(t_mat, 0)
mask <- lower.tri(t_mat, diag = TRUE)

# matrix of rho ^ n contributions
lnrho_arr <- outer(t_mat, lnrho, FUN = "*")

# !! outer produces NA values - subsequent calculations will be wrong as a result
lnrho_arr
#> , , 1, 1
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   NA   NA   NA   NA   NA
#> [2,]   NA   NA   NA   NA   NA
#> [3,]   NA   NA   NA   NA   NA
#> [4,]   NA   NA   NA   NA   NA
#> [5,]   NA   NA   NA   NA   NA
#> 
#> , , 2, 1
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   NA   NA   NA   NA   NA
#> [2,]   NA   NA   NA   NA   NA
#> [3,]   NA   NA   NA   NA   NA
#> [4,]   NA   NA   NA   NA   NA
#> [5,]   NA   NA   NA   NA   NA
#> 
#> , , 3, 1
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]   NA   NA   NA   NA   NA
#> [2,]   NA   NA   NA   NA   NA
#> [3,]   NA   NA   NA   NA   NA
#> [4,]   NA   NA   NA   NA   NA
#> [5,]   NA   NA   NA   NA   NA

^{Created on 2022-12-13 with reprex v2.0.2}

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.2.1 (2022-06-23)
#>  os       macOS Monterey 12.3.1
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Australia/Brisbane
#>  date     2022-12-13
#>  pandoc   2.19.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date (UTC) lib source
#>  base64enc     0.1-3      2015-07-28 [1] CRAN (R 4.2.0)
#>  callr         3.7.3      2022-11-02 [1] CRAN (R 4.2.0)
#>  cli           3.4.1      2022-09-23 [1] CRAN (R 4.2.0)
#>  coda          0.19-4     2020-09-30 [1] CRAN (R 4.2.0)
#>  codetools     0.2-18     2020-11-04 [1] CRAN (R 4.2.1)
#>  crayon        1.5.2      2022-09-29 [1] CRAN (R 4.2.0)
#>  digest        0.6.30     2022-10-18 [1] CRAN (R 4.2.0)
#>  ellipsis      0.3.2      2021-04-29 [1] CRAN (R 4.2.0)
#>  evaluate      0.18       2022-11-07 [1] CRAN (R 4.2.0)
#>  fastmap       1.1.0      2021-01-25 [1] CRAN (R 4.2.0)
#>  fs            1.5.2      2021-12-08 [1] CRAN (R 4.2.0)
#>  future        1.29.0     2022-11-06 [1] CRAN (R 4.2.0)
#>  globals       0.16.2     2022-11-21 [1] CRAN (R 4.2.1)
#>  glue          1.6.2      2022-02-24 [1] CRAN (R 4.2.0)
#>  greta       * 0.4.2.9000 2022-12-12 [1] local
#>  here          1.0.1      2020-12-13 [1] CRAN (R 4.2.0)
#>  highr         0.9        2021-04-16 [1] CRAN (R 4.2.0)
#>  hms           1.1.2      2022-08-19 [1] CRAN (R 4.2.0)
#>  htmltools     0.5.3      2022-07-18 [1] CRAN (R 4.2.0)
#>  jsonlite      1.8.3      2022-10-21 [1] CRAN (R 4.2.0)
#>  knitr         1.41       2022-11-18 [1] CRAN (R 4.2.0)
#>  lattice       0.20-45    2021-09-22 [1] CRAN (R 4.2.1)
#>  lifecycle     1.0.3      2022-10-07 [1] CRAN (R 4.2.0)
#>  listenv       0.8.0      2019-12-05 [1] CRAN (R 4.2.0)
#>  magrittr      2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
#>  Matrix        1.5-3      2022-11-11 [1] CRAN (R 4.2.0)
#>  parallelly    1.32.1     2022-07-21 [1] CRAN (R 4.2.0)
#>  pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.2.0)
#>  png           0.1-8      2022-11-29 [1] CRAN (R 4.2.0)
#>  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.2.0)
#>  processx      3.8.0      2022-10-26 [1] CRAN (R 4.2.0)
#>  progress      1.2.2      2019-05-16 [1] CRAN (R 4.2.0)
#>  ps            1.7.2      2022-10-26 [1] CRAN (R 4.2.0)
#>  purrr         0.3.5      2022-10-06 [1] CRAN (R 4.2.0)
#>  R.cache       0.16.0     2022-07-21 [1] CRAN (R 4.2.0)
#>  R.methodsS3   1.8.2      2022-06-13 [1] CRAN (R 4.2.0)
#>  R.oo          1.25.0     2022-06-12 [1] CRAN (R 4.2.0)
#>  R.utils       2.12.2     2022-11-11 [1] CRAN (R 4.2.0)
#>  R6            2.5.1      2021-08-19 [1] CRAN (R 4.2.0)
#>  Rcpp          1.0.9      2022-07-08 [1] CRAN (R 4.2.0)
#>  reprex        2.0.2      2022-08-17 [1] CRAN (R 4.2.0)
#>  reticulate    1.26       2022-08-31 [1] CRAN (R 4.2.0)
#>  rlang         1.0.6      2022-09-24 [1] CRAN (R 4.2.0)
#>  rmarkdown     2.18       2022-11-09 [1] CRAN (R 4.2.0)
#>  rprojroot     2.0.3      2022-04-02 [1] CRAN (R 4.2.0)
#>  rstudioapi    0.14       2022-08-22 [1] CRAN (R 4.2.0)
#>  sessioninfo   1.2.2      2021-12-06 [1] CRAN (R 4.2.0)
#>  stringi       1.7.8      2022-07-11 [1] CRAN (R 4.2.0)
#>  stringr       1.5.0      2022-12-02 [1] CRAN (R 4.2.0)
#>  styler        1.8.1      2022-11-07 [1] CRAN (R 4.2.0)
#>  tensorflow    2.9.0      2022-05-21 [1] CRAN (R 4.2.0)
#>  tfruns        1.5.1      2022-09-05 [1] CRAN (R 4.2.0)
#>  vctrs         0.5.1      2022-11-16 [1] CRAN (R 4.2.0)
#>  whisker       0.4        2019-08-28 [1] CRAN (R 4.2.0)
#>  withr         2.5.0      2022-03-03 [1] CRAN (R 4.2.0)
#>  xfun          0.35       2022-11-16 [1] CRAN (R 4.2.0)
#>  yaml          2.3.6      2022-10-18 [1] CRAN (R 4.2.0)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.2-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.13 | packaged by conda-forge | (default, Mar 25 2022, 06:05:16)  [Clang 12.0.1 ]
#>  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
#> 
#> ──────────────────────────────────────────────────────────────────────────────

I’ll check in with @nick about this and see what he thinks.

Cheers!

Nick

nick · December 13, 2022, 11:14am

Yeah, it looks like outer is not a function that greta overloads. So when this is run, it treats the greta array as if it were an array. Unfortunately because it’s not a function in greta we can’t issue a warning or anything (short of overloading every function in existence). So maybe the best thing to oin this instance would be to find the tensorflow equivalent and just implement the function.

As a workaround, you might be able to use kronecker() here, for which there is a greta op:

lnrho_arr <- kronecker(t_mat, lnrho)

and then reshape with the required dimensions (NOTE: I didn’t check this reshapes in the correct order)

dim(lnrho_arr) <- c(dim(t_mat), length(lnrho))

Unfortunately, when this prints it looks like you get NAs again (@njtierney something we need to fix with dim assignment, it’s not creating the unknowns object, at least in the old greta), but it actually is a greta array, unlike before, so you can do greta things with it:

calculate(lnrho_arr, nsim = 1)

njtierney · December 14, 2022, 5:47am

OK thanks for that, @nick!

I see two issues to log on github then:

We should implement outer in {greta}
dim<- doesn’t always create the unknown object

Does that sound right to you?

nick · December 14, 2022, 5:59am

Yep!

Here, have some extra characters to meet Discourse’s 20 character minimum.

njtierney · December 20, 2022, 7:01am

OK great, thanks! That is done now on github. And this response is also more than 20 characters now.

BobOHara · January 4, 2023, 10:53am

I finally got some time to look at this, and get it working (not brilliantly, but working well enough to try on the real problem). It looks like having everything as a matrix is easier to work with, as it avoids an array of problems.

Here’s the code I ended up with, it it’s helpful to anyone. I’m sure it’s not perfect.

Simulate some data

NSeries <- 3
NTime <- 100
Abund <- matrix(rpois(NSeries*NTime, 10), ncol=NSeries)
Abund_v <- c(Abund)

Hyperparameters

lnrho = normal(0, sd=0.3, dim=NSeries)
mu0 = normal(0, sd=0.5, dim=NSeries)
sigma = exponential(1, dim=NSeries)

matrix of time modifiers

t_seq <- rev(seq_len(NTime))

t_mat <- outer(t_seq, t_seq, FUN = “-”)
t_mat <- pmax(t_mat, 0)

which elements to include (don’t want upper triangular ones)

mask <- lower.tri(t_mat, diag = TRUE)

matrix of rho ^ n contributions

This code is to test dimensions

lnrho <- 10*(1:3)

lnrho_k <- kronecker(t_mat, lnrho, FUN = “+”)

Write lnrho as a matrix (this should probably be function)

lnrho_k <- kronecker(lnrho, t(t_mat))
mask_k <- kronecker(rep(1, length(lnrho)), mask)

rho_arr <- lnrho_k * mask_k

sigma_arr <- kronecker(ones(NTime,NTime), sigma)
innovations = normal(0, 1, dim = c(NTime*NSeries,NTime))*sigma_arr

Finally write the model for loglam and the likelihood

loglam_arr <- rho_arr*innovations
loglam <- apply(loglam_arr, 1, “sum”) + rep(muQ, each=NTime)
distribution(Abund_v) = poisson(exp(loglam))

Have fun

m <- model(lnrho, sigma, mu0)
mctry <- mcmc(m, chains=2, n_samples=100, warmup = 20)

njtierney · January 9, 2023, 1:53am

Hi @BobOHara!

Happy new year! That’s great news, glad to hear you got it working

Cheers,

Nick