Stan cov_exp_quad performance test

covMat_by_fun.stan

data {
    // N: num unique x values
    int<lower=1> N ;
    // x: x values
    real x[N] ;
    // y: data
    vector[N] y ;
}
transformed data{
	vector[N] zeros ;
	zeros = rep_vector(0,N) ;
}
parameters {
    // noise: sampling noise for replicates
    real<lower=0> noise ;
    // eta: GP parameter, scale of output (y)
    real<lower=0> eta ;
    // inv_rho: GP parameter, inverse lengthscale (lengthscale: x-axis scale of influence between points (wiggliness))
    real<lower=0> inv_rho ;
}
transformed parameters{
    // L: cholesky decomposition of covariance matrix
    matrix[N,N] L ;
    {
        // covMat: covariance matrix
        matrix[N,N] covMat ;
        real rho ;
        real eta_sq_plus_noise ;
        eta_sq_plus_noise = (eta^2) + noise ;
        rho = 1/inv_rho ;
        covMat = cov_exp_quad(x,eta,rho) ;
        for(n in 1:N){
        	covMat[n,n] = eta_sq_plus_noise ;
        }
        L = cholesky_decompose(covMat) ;
    }
}
model {
    // priors on eta as peaked just below 1 with nearly-linear left-tail and slow/fat left tail
    eta ~ weibull(2,1) ;
    // prior on inverse lengthscale as peaked at zero with fat tails
    inv_rho ~ student_t(4,0,1) ;
    // prior on noise (same shape as eta)
    noise ~ weibull(2,1) ;
    // assert sampling of each replicate of y as a GP
    y ~ multi_normal_cholesky(zeros,L) ;
}

covMat_by_hand.stan

data {
    // N: num unique x values
    int<lower=1> N ;
    // x: x values
    real x[N] ;
    // y: data
    vector[N] y ;
}
transformed data{
	vector[N] zeros ;
	zeros = rep_vector(0,N) ;
}
parameters {
    // noise: sampling noise for replicates
    real<lower=0> noise ;
    // eta: GP parameter, scale of output (y)
    real<lower=0> eta ;
    // inv_rho: GP parameter, inverse lengthscale (lengthscale: x-axis scale of influence between points (wiggliness))
    real<lower=0> inv_rho ;
}
transformed parameters{
    // L: cholesky decomposition of covariance matrix
    matrix[N,N] L ;
    {
        // covMat: covariance matrix
        matrix[N,N] covMat ;
        real two_times_rho_sq ;
        real eta_sq;
        real eta_sq_plus_noise ;
        two_times_rho_sq = 2*((1/inv_rho)^2) ;
        eta_sq = eta^2 ;
        eta_sq_plus_noise = eta_sq + noise ;
        for(i in 1:(N-1)){
        	covMat[i,i] = eta_sq_plus_noise ;
        	for(j in (i+1):N){
        		covMat[i,j] = exp( (-( x[i] - x[j] )^2) / two_times_rho_sq ) * eta_sq ;
        		covMat[j,i] = covMat[i,j] ;
        	}
        }
        covMat[N,N] = eta_sq_plus_noise ;
        L = cholesky_decompose(covMat) ;
    }
}
model {
    // priors on eta as peaked just below 1 with nearly-linear left-tail and slow/fat left tail
    eta ~ weibull(2,1) ;
    // prior on inverse lengthscale as peaked at zero with fat tails
    inv_rho ~ student_t(4,0,1) ;
    // prior on noise (same shape as eta)
    noise ~ weibull(2,1) ;
    // assert sampling of each replicate of y as a GP
    y ~ multi_normal_cholesky(zeros,L) ;
}

make_data_and_fit.R

library(rstan)
covMat_by_hand = stan_model(file="covMat_by_hand.stan")
covMat_by_fun = stan_model(file="covMat_by_fun.stan")

N = 200 #number of x-axis grid points
noise = 1 #amount of noise added to each replicate

#define the x-axis grid
x = seq(-10,10,length.out=N)

#define a wiggly function on x
f = scale(sin(x)*dnorm(x,5,8))[,1]

#add noise to f
y = f + rnorm(N,0,noise)

covMat_by_hand_post <- sampling(
	covMat_by_hand
	, data = list(
		x = x
		, N = length(x)
		, y = y
	)
	, iter = 2e3
	, chains = 1
	, cores = 1
	, seed = 1
	, pars = c('inv_rho','eta','noise')
)
alarm()
print(
	covMat_by_hand_post
	, pars = c('inv_rho','eta','noise')
	, digits = 2
	, probs = c(.025,.975)
)
get_elapsed_time(covMat_by_hand_post)

covMat_by_fun_post <- sampling(
	covMat_by_fun
	, data = list(
		x = x
		, N = length(x)
		, y = y
	)
	, iter = 2e3
	, chains = 1
	, cores = 1
	, seed = 1
	, pars = c('inv_rho','eta','noise')
)
alarm()
print(
	covMat_by_fun_post
	, pars = c('inv_rho','eta','noise')
	, digits = 2
	, probs = c(.025,.975)
)
get_elapsed_time(covMat_by_fun_post)