Hi @Kartik!
Just wanted to follow up say that we have fixed this in the latest release of greta.dynamics (https://cran.r-project.org/web/packages/greta.dynamics/), thanks for bringing this to our attention!
Here’s a reproducible example:
library(greta.dynamics)
#> Loading required package: 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
greta_sitrep()
#> ℹ checking if python available
#> ✔ python (v3.10) available
#>
#> ℹ checking if TensorFlow available
#> ✔ TensorFlow (v2.15.0) available
#>
#> ℹ checking if TensorFlow Probability available
#> ✔ TensorFlow Probability (v0.23.0) available
#>
#> ℹ checking if greta conda environment available
#> ✔ greta conda environment available
#>
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#>
#> ℹ greta is ready to use!
library(deSolve)
LVmod <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
Ingestion <- rIng * Prey * Predator
GrowthPrey <- rGrow * Prey * (1 - Prey / K)
MortPredator <- rMort * Predator
dPrey <- GrowthPrey - Ingestion
dPredator <- Ingestion * assEff - MortPredator
return(list(c(dPrey, dPredator)))
})
}
pars <- c(
rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2, # /day, mortality rate of predator
assEff = 0.5, # -, assimilation efficiency
K = 10
) # mmol/m3, carrying capacity
yini <- c(Prey = 1, Predator = 2)
times <- seq(0, 30, by = 1)
out <- ode(yini, times, LVmod, pars)
# simulate observations
jitter <- rnorm(2 * length(times), 0, 0.1)
y_obs <- out[, -1] + matrix(jitter, ncol = 2)
# fit a greta model to infer the parameters from this simulated data
# greta version of the function
lotka_volterra <- function(y, t, rIng, rGrow, rMort, assEff, K) {
Prey <- y[1, 1]
Predator <- y[1, 2]
Ingestion <- rIng * Prey * Predator
GrowthPrey <- rGrow * Prey * (1 - Prey / K)
MortPredator <- rMort * Predator
dPrey <- GrowthPrey - Ingestion
dPredator <- Ingestion * assEff - MortPredator
cbind(dPrey, dPredator)
}
# priors for the parameters
rIng <- uniform(0, 2) # /day, rate of ingestion
rGrow <- uniform(0, 3) # /day, growth rate of prey
rMort <- uniform(0, 1) # /day, mortality rate of predator
assEff <- uniform(0, 1) # -, assimilation efficiency
K <- uniform(0, 30) # mmol/m3, carrying capacity
# initial values and observation error
y0 <- uniform(0, 5, dim = c(1, 2))
obs_sd <- uniform(0, 1)
# solution to the ODE
y <- ode_solve(lotka_volterra, y0, times, rIng, rGrow, rMort, assEff, K)
# sampling statement/observation model
distribution(y_obs) <- normal(y, obs_sd)
# we can use greta to solve directly, for a fixed set of parameters (the true
# ones in this case)
values <- c(
list(y0 = t(1:2)),
as.list(pars)
)
vals <- calculate(y, values = values)[[1]]
vals
#> [,1] [,2]
#> [1,] 1.000000 2.000000
#> [2,] 1.626868 1.863281
#> [3,] 2.488485 1.871156
#> [4,] 3.408719 2.058952
#> [5,] 4.060276 2.457740
#> [6,] 4.200195 3.055239
#> [7,] 3.852383 3.753852
#> [8,] 3.242341 4.386750
#> [9,] 2.612406 4.810123
#> [10,] 2.101931 4.978732
#> [11,] 1.746161 4.934839
#> [12,] 1.528342 4.754134
#> [13,] 1.420247 4.507028
#> [14,] 1.398437 4.245808
#> [15,] 1.446404 4.005420
#> [16,] 1.551344 3.808030
#> [17,] 1.699583 3.666968
#> [18,] 1.872662 3.588978
#> [19,] 2.046175 3.574837
#> [20,] 2.192831 3.619023
#> [21,] 2.289637 3.709183
#> [22,] 2.325193 3.826968
#> [23,] 2.302773 3.950750
#> [24,] 2.237499 4.060039
#> [25,] 2.149706 4.139833
#> [26,] 2.058567 4.183017
#> [27,] 1.978655 4.190254
#> [28,] 1.918413 4.167923
#> [29,] 1.880715 4.125443
#> [30,] 1.867847 4.073026
#> [31,] 1.876775 4.020200
Created on 2024-11-26 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.1
#> system aarch64, darwin20
#> ui X11
#> language (EN)
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Australia/Hobart
#> date 2024-11-26
#> 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-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.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.2)
#> crayon 1.5.3 2024-06-20 [1] CRAN (R 4.4.0)
#> deSolve * 1.40 2023-11-27 [1] CRAN (R 4.4.0)
#> digest 0.6.37 2024-08-19 [1] CRAN (R 4.4.1)
#> evaluate 1.0.1 2024-10-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 2024-11-06 [1] local
#> greta.dynamics * 0.2.2 2024-11-13 [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.9 2024-09-20 [1] CRAN (R 4.4.1)
#> knitr 1.49 2024-11-08 [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-1 2024-10-18 [2] CRAN (R 4.4.2)
#> parallelly 1.39.0 2024-11-07 [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.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.8.1 2024-10-28 [1] CRAN (R 4.4.1)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.4.0)
#> Rcpp 1.0.13-1 2024-11-02 [1] CRAN (R 4.4.1)
#> reprex 2.1.1 2024-07-06 [1] CRAN (R 4.4.0)
#> reticulate 1.40.0 2024-11-15 [1] CRAN (R 4.4.1)
#> rlang 1.1.4 2024-06-04 [1] CRAN (R 4.4.0)
#> 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.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.2 2024-10-28 [1] CRAN (R 4.4.1)
#> xfun 0.49 2024-10-31 [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
#>
#> ─ 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
#>
#> ──────────────────────────────────────────────────────────────────────────────
Let us know how you go 
