fix kernel spectral densities

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

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

vignettes/gaussian_process_implementation_details.Rmd

@@ -104,28 +104,28 @@ 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}
 
 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 \omega^2 \right)
 \end{equation}
 
 The functions $\phi_{j}(x)$ are the eigenfunctions of the Laplace operator,