spectral densities of Matern Kernels of order != 3/2

[sbfnk]

I don't think the spectral densities of the Matern Kernels introduced in #742 are correct when the order is not 3/2.

From Gabriel Riutort-Mayol et al. this should be

factor = 2 * alpha * sqrt(pi) * tgamma(nu + 0.5) * pow(2 * nu, nu);
denom = tgamma(nu) * pow(rho, 2 * nu) * pow(2 * nu / rho^2 + (pi() / (2 / L) * indices)^2, nu + 0.5);

which collapses to the current specification only if nu = 1.5

factor = 2 * alpha * pow(sqrt(2 * nu) / rho, nu);
denom = (sqrt(2 * nu) / rho)^2 + pow((pi() / 2 / L) * indices, nu + 0.5);

In particular, of the three choices we offer (1/2, 3/2, 5/2), nu = 1.5 is the only one where nu + 0.5 is even and thus the square root can be done conveniently analytically. For the other two I think this has to be done numerically.

To fix this we could

1. implement the more general form
2. keep this as a specific function for the Matern 3/2 kernel and implement separate functions for orders 1/2 and 5/2, i.e. what we had before.

I suspect option (2) might be more computationally efficient and, as we're unlikely to support other orders in the near future, not too cumbersome.

[seabbs]

Linking to code:

EpiNow2/inst/stan/functions/gaussian_process.stan

Line 36 in 9cf7ad2

real factor = 2 * alpha * pow(sqrt(2 * nu) / rho, nu);

I need to look at this in a bit but on the face of it you seem correct and yes agree we should go back to having fixed solutions (i.e 2)

[sbfnk]

My cheat sheet for the general version is:

nu = 1/2: tgamma(nu) = sqrt(pi())       tgamma(nu + 0.5) = 1   pow(2 * nu, nu) = 1
nu = 3/2: tgamma(nu) = 1/2 sqrt(pi())   tgamma(nu + 0.5) = 1   pow(2 * nu, nu) = 3^(3/2)
nu = 5/2: tgamma(nu) = 3/4 sqrt(pi())   tgamma(nu + 0.5) = 2   pow(2 * nu, nu) = 5^(5/2)