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

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="*")

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))

1 Like

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 = "*")

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
distribution(Abund) <- poisson(exp(loglam), dim = c(NSeries, NTime))
``````

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
#>
#> ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
``````

Apologies, N should just be NTime.

1 Like

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
#>
#> ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
``````

Cheers!

Nick

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)
``````
1 Like

OK thanks for that, @nick!

I see two issues to log on github then:

1. We should implement `outer` in {greta}
2. `dim<-` doesnβt always create the unknown object

Does that sound right to you?

Yep!

Here, have some extra characters to meet Discourseβs 20 character minimum.

1 Like

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

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)

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

lnrho_k <- kronecker(lnrho, t(t_mat))

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)

1 Like

Hi @BobOHara!

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

Cheers,

Nick