Hi All,
I have been reading Aki’s GP demonstration with Stan (Motorcycle data) and the corresponding paper on Hilbert space approximate GPs. In the Motorcycle demonstration, Aki de-means the columns of the PHI matrix (basis functions).
Compared to non-centered basis functions, the prior with centered basis functions seems to have lower variance toward the center of the domain. As far as I can tell, non-centered basis functions are more representative of the exact GP here (constant variance). My question is: why are the columns of PHI de-meaned? Is this best practice or a quirk of this case study?
I have included Stan and R code below to draw realisations from the approximate GP prior under both cases.
gb_question.stan:
functions {
// Spectral densities.
vector diagSPD_EQ(real alpha, real rho, real L, int M) {
vector[M] j = linspaced_vector(M, 1, M);
vector[M] lambda = (j * pi() / (2*L))^2;
vector[M] s = (alpha^2) * sqrt(2*pi()) * rho * exp(-0.5 * (rho^2) * lambda);
return sqrt(s);
}
// Basis functions [N x M]
matrix PHI(int N, int M, real L, vector x, int demean) {
vector[M] j = linspaced_vector(M, 1, M);
vector[M] sqrt_lambda = j * pi() / (2 * L);
matrix[N, M] x_L = rep_matrix(x + L, M);
matrix[N, M] phi = sin(diag_post_multiply(x_L, sqrt_lambda)) / sqrt(L);
if (demean)
for (m in 1:M) phi[, m] = phi[, m] - mean(phi[, m]);
return phi;
}
}
data {
int<lower=1> N; // number of observations
vector[N] x; // univariate covariate
real<lower=0> c_f; // boundary factor
int<lower=1> M_f; // number of basis functions
real<lower=0> lengthscale_f; // lengthscale of f
real<lower=0> sigma_f; // scale of f
int<lower = 0, upper = 1> center_phi; // De-mean basis functions?
}
transformed data {
// Normalize data
real xmean = mean(x);
real xsd = sd(x);
vector[N] xn = (x - xmean)/xsd;
// Boundary value
real L_f = c_f * max(xn);
// Basis functions for f
matrix[N,M_f] PHI_f = PHI(N, M_f, L_f, xn, center_phi);
// spectral densities
vector[M_f] diagSPD_f = diagSPD_EQ(sigma_f, lengthscale_f, L_f, M_f);
}
generated quantities {
// Random draws from f
vector[N] f = PHI_f * (diagSPD_f .* to_vector(normal_rng(0, rep_vector(1, M_f))));
}
R code calling Stan:
fileq <- "./output/gb_question.stan"
# Compile model
modq <- cmdstan_model(stan_file = fileq, include_paths = ".")
# Input data
dat <- list(x = seq(0,1,length.out=100),
N = 100,
c_f = 1.5,
M_f = 40,
sigma_f = 1,
lengthscale_f = 1,
center_phi = 1 # bool for centering basis functions
)
# Centered phi
s_c <- modq$sample(data=dat, fixed_param=TRUE,
iter_warmup = 0, iter_sampling=1000)
# Non-centered phi
dat$center_phi <- 0
s_nc <- modq$sample(data=dat, fixed_param=TRUE,
iter_warmup = 0, iter_sampling=1000)
R code to visualise 90% density intervals.
# Plotting helper
f_q <- function(s, x=xn, id="x") {
f <- s$draws('f')
fmat <- matrix(0, nrow=dim(f)[1], ncol=dim(f)[3])
for (i in 1:dim(f)[1]) fmat[i, ] <- as.numeric(f[i, 1, ])
u <- apply(fmat, 2, quantile, prob=0.95)
l <- apply(fmat, 2, quantile, prob=0.05)
data.frame(x = x, u=u, l=l, phi=id)
}
# Plotting
xn <- as.numeric(scale(dat$x))
dq <- rbind(f_q(s_c, id="centered bf"), f_q(s_nc, id="non-centered bf"))
dq %>%
ggplot(aes(x=x, ymin=l, ymax=u, lty=phi, col=phi, fill=phi)) +
geom_ribbon(lwd=0.8, alpha = 0.3) +
theme_classic() +
labs(title = "90% Coverage Interval", y = "f") +
theme(legend.text = element_text(size=12),
legend.position = 'top',
legend.title = element_blank()
)
Any advice greatly appreciated. I’m guessing I have missed something simple.