Hi,
I am very new to Stan (Rstan). I was recently trying to model nngp using Lu’s example. I am giving here the full code as I added a weight factor while computing the covariance matrix. Initially I faced a problem on the new covariance matrix being “non-positive definite” in cholesky_decompose() after factoring in the weights. After I read few posts, I decided to force it to become pd and symmetric (which it already was). Although the problem is gone when I am running stan using NUTS (it has become slower though), I am getting an error while using the vb as:
Error in sampler$call_sampler(c(args, dotlist)) : Expecting a string vector: [type=NULL; required=STRSXP]
Can you please take a look and let me know how can I use variational bayes with advi here. I also don’t know if it is possible to make stan with nuts run faster. I am copying the R code as well. You can download the data from this link.
Thanks,
Subhomoy
/* Latent NNGP model*/
functions{
matrix make_positive_definite(matrix m){
int N = cols(m);
matrix[N,N] es= eigenvectors_sym(m);
vector[N] esv= eigenvalues_sym(m);
real delta = 2 * 1e-5;
vector[N] tau;
matrix[N,N] dm;
matrix[N,N] mdm;
for(i in 1:N){
tau[i] = fmax(0,delta-esv[i]);
}
dm = es * diag_matrix(tau) * es';
return m + dm;
}
real nngp_w1_lpdf(vector w1, real sigmasq1, real phi1, matrix NN_dist1,matrix NN_cor1,
matrix NN_distM1, matrix NN_corM1, int[,] NN_ind1, int N1, int M){
vector[N1] V;
vector[N1] I_Aw = w1;
int dim;
int h;
for (i in 2:N1) {
matrix[ i < (M + 1)? (i - 1) : M, i < (M + 1)? (i - 1): M]
iNNdistM;
matrix[ i < (M + 1)? (i - 1) : M, i < (M + 1)? (i - 1): M]
iNNCholL;
vector[ i < (M + 1)? (i - 1) : M] iNNcorr;
vector[ i < (M + 1)? (i - 1) : M] v;
row_vector[i < (M + 1)? (i - 1) : M] v2;
dim = (i < (M + 1))? (i - 1) : M;
if(dim == 1){iNNdistM[1, 1] = 1;}
else{
h = 0;
for (j in 1:(dim - 1)){
for (k in (j + 1):dim){
h = h + 1;
iNNdistM[j, k] = exp(- phi1 * NN_distM1[(i - 1), h] * NN_corM1[(i - 1), h]);
iNNdistM[k, j] = iNNdistM[j, k];
}
}
for(j in 1:dim){
iNNdistM[j, j] = 1;
}
}
iNNdistM = 0.5*(iNNdistM + iNNdistM');
iNNdistM = make_positive_definite(iNNdistM);
iNNCholL = cholesky_decompose(iNNdistM);
iNNcorr = to_vector(exp(- phi1 * (NN_dist1[(i - 1), 1:dim] .* NN_cor1[(i - 1), 1:dim]) ) );
v = mdivide_left_tri_low(iNNCholL, iNNcorr);
V[i] = 1 - dot_self(v);
v2 = mdivide_right_tri_low(v', iNNCholL);
I_Aw[i] = I_Aw[i] - v2 * w1[NN_ind1[(i - 1), 1:dim]];
}
V[1] = 1;
return - 0.5 * ( 1 / sigmasq1 * dot_product(I_Aw, (I_Aw ./ V)) +
sum(log(V)) + N1 * log(sigmasq1));
}
}
data {
int<lower=1> N1;
int<lower=1> N;
int<lower=1> M;
int<lower=1> P;
vector[N] Y;
matrix[N, P + 1] X;
int NN_ind1[N1 - 1, M];
matrix[N1 - 1, M] NN_dist1;
matrix[N1 - 1, (M * (M - 1) / 2)] NN_distM1;
matrix[N1 - 1, M] NN_cor1;
matrix[N1 - 1, (M * (M - 1) / 2)] NN_corM1;
matrix[N,N1] H1;
vector[P + 1] uB;
matrix[P + 1, P + 1] VB;
real ss1;
real st;
real ap1;
real bp1;
}
transformed data {
cholesky_factor_cov[P + 1] L_VB;
L_VB = cholesky_decompose(VB);
}
parameters{
vector[P + 1] beta;
real<lower = 0> sigma1;
real<lower = 0> tau;
real<lower = 0> phi1;
vector<lower = 0>[N1] w1;
}
transformed parameters {
real sigmasq1 = square(sigma1);
real tausq = square(tau);
}
model{
beta ~ multi_normal_cholesky(uB, L_VB);
phi1 ~ gamma(ap1, bp1);
sigma1 ~ normal(0, ss1);
tau ~ normal(0, st);
w1 ~ nngp_w1(sigmasq1, phi1, NN_dist1, NN_cor1, NN_distM1, NN_corM1, NN_ind1, N1, M);
Y ~ normal(X * beta + H1 * w1, tau);
}
# R
rm(list = ls())
library(rstan)
#library(shinystan)
#library(spNNGP)
#------------------------------ NNGP data load ---------------------------------#
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
load(file.choose())
#------------------------------ initialization vals & other vars ---------------#
myinits <- list(list(beta = c(1, 2, 2), sigma1 = 10, tau = 10, phi1 = 100))
# list(beta = c(2, 1, 3), sigma1 = 2, tau = 1, phi1 = .2),
# list(beta = c(0.5, 1.5, .8), sigma1 = 3, tau = 0.5, phi1 = .5))
parameters <- c("beta", "sigmasq1", "tausq", "phi1")
#------------------------------ run stan Variational Bayes ---------------------#
stan_model <- stan_model(file = file.choose())
# ADVI
stan_fit_vb <- vb(stan_model, data=data_nngp, output_samples=1000, seed=123,init=myinits, pars=parameters)
#------------------------------ run MCMC with NUTS -----------------------------#
samples <- stan(
file = "~/nngp_test.stan",
data = data_nngp,
init = myinits,
pars = parameters,
iter = 1000,
chains = 1,
thin = 100,
cores=1,
seed = 123
)