New covariance function in rstan

[lionel68]

Summary:

Based on my understanding of https://github.com/stan-dev/stan/pull/2739, there should be new covariance function available to use via rstan.

Description:

It seems that these functions are not yet available.

Reproducible Steps:

Simple GP simulation model from the Stan manual:

data {

    int<lower=1> N;

    int<lower=1> D; //dimension of the GP
    vector[D] x[N]; //point where the GP is estimated
}

transformed data {
    matrix[N, N] L; //cholesky matrix
        vector[N] mu = rep_vector(0, N); //assume 0-mean
    matrix[N, N] K = gp_matern32_cov(x, 1.0, 1.0);

    for(n in 1:N)

        K[n, n] = K[n, n] + 0.1; // to ensure positive-definite matrix
    L = cholesky_decompose(K);

}

parameters {

    vector[N] eta;

}

model {

    eta ~ std_normal();

}
generated quantities {
    vector[N] y;
    y = mu + L * eta;
}

And in R:

N <- 100
D <- 2
x <- matrix(runif(N * D), nrow = N, ncol = D)
m <- stan("GP_mat32_sim.stan", data = list(x = x, N = 100, D = 2))

Current Output:

SYNTAX ERROR, MESSAGE(S) FROM PARSER: No matches for:

gp_matern32_cov(vector[ ], real, real)

Function gp_matern32_cov not found. error in 'model1b50107f4b65_GP_mat32_sim' at line 8, column 49

 6: transformed data {
 7:     matrix[N, N] L; //cholesky matrix
 8:     matrix[N, N] K = gp_matern32_cov(x, 1.0, 1.0);
                                                    ^
 9:     for(n in 1:N)

---

Error in stanc(file = file, model_code = model_code, model_name = model_name, : failed to parse Stan model 'GP_mat32_sim' due to the above error.

Expected Output:

NA

RStan Version:

rstan_2.19.2

R Version:

R version 3.6.1 (2019-07-05)

Operating System:

Win10

[jgabry]

Unfortunately I think this function was released in Stan v2.20 (if I'm not mistaken) and so it's not yet available in rstan (still at v2.19) until rstan 2.20.

However all hope is not lost! We've started working on a package cmdstanr, which lets you run cmdstan from R and you can use the latest Stan features without waiting for an rstan release. The package is still in its early stages and may change a lot, but I managed to get your code to run (I did have to remove the generated quantities block though because it uses a variable mu that isn't defined anywhere). Here's how you can try it out if you want (if you do try and have trouble let me know at the cmdstanr repo):

# install package from github
remotes::install_github("stan-dev/cmdstanr")

# if you don't already have cmdstan 2.20 you can run this
# (currently assumes you have C++ toolchain, but you should already if  you're able to use rstan)
install_cmdstan() 

# if you already had cmdstan 2.20 installed you can set this to that path, 
# otherwise if you used install_cmdstan() with default settings this shoud work:
set_cmdstan_path(file.path(Sys.getenv("HOME"), ".cmdstanr/cmdstan"))

# save your Stan file and pass that path to cmdstan_model()
# (i had to remove the generated quantities block of your program because 
# it uses a variable mu that hasn't been defined anywhere)
mod <- cmdstan_model("/Users/jgabry/Desktop/cmdstanr-test/gp/gp.stan")

# compile the model
mod$compile()

# prepare data and pass to mod$sample()
N <- 100
D <- 2
x <- matrix(runif(N * D), nrow = N, ncol = D)
standata <- list(x = x, N = N, D = D)

# runs 1 chain by default (like cmdstan) and can only run >1 chain
# sequentially (eventually we'll add parallelization)
fit <- mod$sample(data = standata, num_chains = 4)

# the object 'fit' doesn't do much yet (still under development) but 
# you can create a regular stanfit object from the csv output files
stanfit <- rstan::read_stan_csv(fit$output_files())

Hope that helps!

[lionel68]

Hi jgabry, Thanks for your swift reply. I just updated the code above with the missing mu line. I'll give your solution via cmdstan a try, but do you have some idea when the next rstan version will come out?

[jgabry]

Unfortunately I don’t know yet about the timing of the rstan release, but I would assume a few weeks at a minimum. @bgoodri would know better, but I’m not sure there’s a target date yet.