My experience is that AR unobserved parameters and normal noise can
be multi-modal in a way that nearly converges but because it’s so high-dimensional
chains can get trapped s.t. if you look at the trace plots they’ll look like
they are mixing but a lot of parameters will be slightly off. Poisson noise
doesn’t seem to have that problem. I think this version shows that problem sometimes:
data {
int T; // Number of time points.
int K; // Number of coefficients.
int M;
vector[T] y; // Measurement at each T
vector[M] kappa; // scaling.
}
parameters {
vector[K] beta;
vector<lower=0>[2] sigma;
vector[T] zeta;
}
transformed parameters {
vector[T] mu;
mu[1] = zeta[1];
for (t in 2:T) {
mu[t] = zeta[t]*sigma[1] + inv_logit(beta[1]) * mu[t-1];
// print("t: ", t, ", mu[t]: ", mu[t], ", zeta[t]: ", zeta[t],
// ", sigma[1]: ", sigma[1], ", beta[1]: ", beta[1],
// ", mu[t-1]: ", mu[t-1]
// );
}
}
model {
target += normal_lpdf(beta | 0.0, 1.0);
target += gamma_lpdf(sigma | 4.0, 10);
target += normal_lpdf(zeta | 0.0, 1);
for (t in 1:T) {
target += normal_lpdf(y[t] | mu[t], sigma[2]);
}
}
Here’s a run script for it:
T <- 400
J <- 2
x <- matrix(data=0, nrow=J, ncol=T)
y <- matrix(data=0, nrow=J, ncol=T)
beta <- runif(n=J, min=0, max=.95)
for (j in 1:J) {
for (t in 2:T) {
x[j,t] <- beta[j]*x[j,t-1] + rnorm(1,0,.7)
y[j,t] <- x[j,t] + rnorm(1,0,.3)
}
}
library(rstan)
m_ar_nc <- rstan::stan_model(file='non-centered-ar+noise.stan')
ml <- list()
for (j in 1:J) {
ml[[j]] <- rstan::sampling(object=m_ar_nc,
data=list(T=T, K=1, y=y[j,], M=1, kappa=array(1)),
control=list(adapt_delta=.90),
chains=5, iter=2000, warmup=1500)
}
o <- list()
for (j in 1:J) {
o[[j]] <- ml[[j]] %>% rstan::extract() %>% data.frame %>%
dplyr::mutate(model=as.character(j), iteration=1:nrow(.))
}
ob <- data.table::rbindlist(o)
pl <- ggplot(
data=ob,
aes(x=iteration, y=plogis(beta), colour=model)
) + geom_line() + geom_abline(intercept=beta, slope=0) +
theme_minimal() + facet_wrap( ~ model)
print(pl)
Some of these plots are useful:
sigmas <- rstan::extract(ml[[1]])$sigma %>% data.frame(check.names=FALSE) %>%
dplyr::mutate(iteration=1:2500, chain=sort(rep(1:5,500))) %>%
tidyr::gather(sigma, value, -iteration, -chain) %>%
dplyr::mutate(parameter=paste('sigma', sigma))
pl_sigmas <- ggplot(data=sigmas,
aes(x=iteration, y=value^2, colour=factor(sigma), shape=factor(chain))) +
geom_point() + theme_minimal()
betas <- rstan::extract(ml[[1]])$beta %>% data.frame(check.names=FALSE) %>%
dplyr::mutate(iteration=1:2500, chain=sort(rep(1:5,500))) %>%
tidyr::gather(beta, value, -iteration, -chain) %>%
dplyr::mutate(parameter=paste('beta', beta))
pl_betas <- ggplot(data=betas,
aes(x=iteration, y=plogis(value), colour=factor(beta), shape=factor(chain))) +
geom_point() + theme_minimal()
zetas <- rstan::extract(ml[[1]])$zeta %>%
data.frame(check.names=FALSE) %>%
dplyr::mutate(iteration=rep(1:500,5), chain=sort(rep(1:5,500))) %>%
tidyr::gather(zeta, value, -iteration, -chain)
pl_zetas <- ggplot(data=zetas,
aes(x=as.numeric(zeta), y=value, colour=factor(chain))
) + geom_point(size=0.001) + theme_minimal()
mus <- rstan::extract(ml[[1]])$mu %>%
data.frame(check.names=FALSE) %>%
dplyr::mutate(iteration=rep(1:500,5), chain=sort(rep(1:5,500))) %>%
tidyr::gather(mu, value, -iteration, -chain)
pl_mus <- ggplot(data=mus,
aes(x=as.numeric(mu), y=value, colour=factor(chain))
) + geom_point(size=0.001) + theme_minimal()