From 6ad0e8f75a51c3ad66e0d4724e41c107aaa7f11a Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Thu, 24 Oct 2024 12:20:28 +0100 Subject: [PATCH 1/2] fix kernel spectral densities --- inst/stan/functions/gaussian_process.stan | 4 ++-- tests/testthat/test-stan-guassian-process.R | 4 ++-- vignettes/gaussian_process_implementation_details.Rmd | 8 ++++---- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/inst/stan/functions/gaussian_process.stan b/inst/stan/functions/gaussian_process.stan index 94f488398..ab5ee5eb7 100644 --- a/inst/stan/functions/gaussian_process.stan +++ b/inst/stan/functions/gaussian_process.stan @@ -33,7 +33,7 @@ vector diagSPD_EQ(real alpha, real rho, real L, int M) { vector diagSPD_Matern12(real alpha, real rho, real L, int M) { vector[M] indices = linspaced_vector(M, 1, M); real factor = 2; - vector[M] denom = rho * ((1 / rho)^2 + (pi() / 2 / L) * indices); + vector[M] denom = rho * ((1 / rho)^2 + pow(pi() / 2 / L * indices, 2)); return alpha * sqrt(factor * inv(denom)); } @@ -65,7 +65,7 @@ vector diagSPD_Matern32(real alpha, real rho, real L, int M) { vector diagSPD_Matern52(real alpha, real rho, real L, int M) { vector[M] indices = linspaced_vector(M, 1, M); real factor = 3 * pow(sqrt(5) / rho, 5); - vector[M] denom = 2 * (sqrt(5) / rho)^2 + pow((pi() / 2 / L) * indices, 3); + vector[M] denom = 2 * pow((sqrt(5) / rho)^2 + pow((pi() / 2 / L) * indices, 2), 3); return alpha * sqrt(factor * inv(denom)); } diff --git a/tests/testthat/test-stan-guassian-process.R b/tests/testthat/test-stan-guassian-process.R index 61be6f479..0a4717065 100644 --- a/tests/testthat/test-stan-guassian-process.R +++ b/tests/testthat/test-stan-guassian-process.R @@ -36,7 +36,7 @@ test_that("diagSPD_Matern functions return correct dimensions and values", { # Check specific values for known inputs indices <- linspaced_vector(M, 1, M) factor12 <- 2 - denom12 <- rho * ((1 / rho)^2 + (pi / 2 / L) * indices) + denom12 <- rho * ((1 / rho)^2 + (pi / 2 / L * indices)^2) expected_result12 <- alpha * sqrt(factor12 / denom12) expect_equal(result12, expected_result12, tolerance = 1e-8) @@ -58,7 +58,7 @@ test_that("diagSPD_Matern functions return correct dimensions and values", { # Check specific values for known inputs indices <- linspaced_vector(M, 1, M) factor52 <- 3 * (sqrt(5) / rho)^5 - denom52 <- 2 * (sqrt(5) / rho)^2 + ((pi / 2 / L) * indices)^3 + denom52 <- 2 * ((sqrt(5) / rho)^2 + (pi / 2 / L * indices)^2)^3 expected_result52 <- alpha * sqrt(factor52 / denom52) expect_equal(result52, expected_result52, tolerance = 1e-8) }) diff --git a/vignettes/gaussian_process_implementation_details.Rmd b/vignettes/gaussian_process_implementation_details.Rmd index 7ea2e7a29..69d575bb4 100644 --- a/vignettes/gaussian_process_implementation_details.Rmd +++ b/vignettes/gaussian_process_implementation_details.Rmd @@ -104,7 +104,7 @@ S_k(x) = \alpha^2 \frac{2 \sqrt{\pi} \Gamma(\nu + 1/2) (2\nu)^\nu}{\Gamma(\nu) \ For $\nu = 3 / 2$ (the default in `EpiNow2`) this simplifies to \begin{equation} - S_k(x) = + S_k(\omega) = \alpha^2 \frac{4 \left(\sqrt{3} / \rho \right)^3}{\left(\left(\sqrt{3} / \rho\right)^2 + \omega^2\right)^2} = \left(\frac{2 \alpha \left(\sqrt{3} / \rho \right)^{\frac{3}{2}}}{\left(\sqrt{3} / \rho\right)^2 + \omega^2} \right)^2 \end{equation} @@ -112,20 +112,20 @@ For $\nu = 3 / 2$ (the default in `EpiNow2`) this simplifies to For $\nu = 1 / 2$ it is \begin{equation} -S_k(x) = \alpha^2 \frac{2}{\rho \left(1 / \rho^2 + \omega^2\right)} +S_k(\omega) = \alpha^2 \frac{2}{\rho \left(1 / \rho^2 + \omega^2\right)} \end{equation} and for $\nu = 5 / 2$ it is \begin{equation} -S_k(x) = +S_k(\omega) = \alpha^2 \frac{3 \left(\sqrt{5} / \rho \right)^5}{2 \left(\left(\sqrt{5} / \rho\right)^2 + \omega^2\right)^3} \end{equation} In the case of a squared exponential the spectral kernel density is given by \begin{equation} -S_k(x) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 w^2 \right) +S_k(\omega) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 w^2 \right) \end{equation} The functions $\phi_{j}(x)$ are the eigenfunctions of the Laplace operator, From 328d6c739241a01adad30e79c8b8e20b62e3345d Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Thu, 24 Oct 2024 12:22:33 +0100 Subject: [PATCH 2/2] sneak in typo fix --- vignettes/gaussian_process_implementation_details.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/gaussian_process_implementation_details.Rmd b/vignettes/gaussian_process_implementation_details.Rmd index 69d575bb4..a92091162 100644 --- a/vignettes/gaussian_process_implementation_details.Rmd +++ b/vignettes/gaussian_process_implementation_details.Rmd @@ -125,7 +125,7 @@ S_k(\omega) = In the case of a squared exponential the spectral kernel density is given by \begin{equation} -S_k(\omega) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 w^2 \right) +S_k(\omega) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 \omega^2 \right) \end{equation} The functions $\phi_{j}(x)$ are the eigenfunctions of the Laplace operator,