Hello everyone I’m quite new to the stan programing and I’m at the moment attempting to fit time series data to the sum of two sinusoids (possible more in the future). My main problem is label switching as running one chain of the model works fine, but as soon as 4 chains are running the model breaks down.
My model has the following structure:
y(t)=\sum_{j=1}^{k}(A_ j\cdot cos(2\pi f_jt+2\pi\phi_j )\cdot \exp(-a_jt) ).
I assume the following priors for each parameter (a,A,\phi,f) :
a\sim N(0,2.2)
A\sim N(0,\Sigma)
\phi\sim U(0,1)
f\sim U(0,15)
My likelihood I presume to be a normal distribution with a single variance parameter and a sum of two mean values:
y|(a,A,\phi,f)=N(\mu_1+\mu_2,\sigma^2)
where \mu_1=A_ 1\cdot cos(2\pi f_1t+2\pi\phi_1 )\cdot \exp(-a_1t) and \mu_2=A_ 2\cdot cos(2\pi f_2t+2\pi\phi_2 )\cdot \exp(-a_2t)
and finally \sigma^2\sim InvGamma(2,10)
I simulate my data as the following: (utilizing rstan)
# two sinusoid models #
# data and model params
N=100 # points
dt=0.1 # dwell time
t_max=0.99
# time points
j_idx=seq(0,t_max,length.out = N)
time=j_idx*1
iu=as.complex(1i) # complex imaguinary unit
# component 1 parms
A_1=1
phi_1=0.2
f_1=10
alpha_1=3
# component 2 parms
A_2=5
phi_2=0.5
f_2=8
alpha_2=1
# component 1 and 2
y1=cos(2*pi*time*f_1+2*pi*phi_1)*exp(-alpha_1*time)*A_1
y2=cos(2*pi*time*f_2+2*pi*phi_2)*exp(-alpha_2*time)*A_2
# full data
y=y1+y2
# add data noise
y_noise=y+rnorm(length(y),0,0.1) # data with noise
library(loo)
library(StanHeaders)
library(rstan)
library(RcppEigen)
library(bayesplot)
library(gridExtra)
library(rstanarm)
# data set up for rstan
dat=list(N = length(y_noise),
time=c(t(time)),
k=2,
y = c(t(Re(y_noise))))
# fitting with rstan - note set your own working directory or paste in the stan model
fit1 <- stan(file = 'testing stan_2_15_5M.stan',cores=1,chains = 1, data = dat,iter=2000,control = list(adapt_delta=0.8,max_treedepth=10)) # using 1 chain fit is fine
fit2<- stan(file = 'testing stan_2_15_5M.stan',cores=4,chains = 4, data = dat,iter=2000,control = list(adapt_delta=0.8,max_treedepth=10)) # using 4 chains label switching
The code of the stan script is as of the following
data{
int k; // no of signals
int N;// total number of points
vector[N] time; // observed time points
vector[N] y; // observed signal values
}
parameters{
vector <lower=0> [k] A; // signal amplitude
vector <lower=0> [k] a; // signal decay constant
vector <lower=0> [k] phi; // phase of signals
vector <lower=0> [k] f; // frequency of signals
real<lower=0> sigma; // signal error
}
model {
vector[N] B1;
vector[N] B2;
for (i in 1:N){
B1[i]=A[1]*cos(2*pi()*time[i]*f[1]+2*pi()*phi[1] )*exp(-a[1]*time[i]);
B2[i]=A[2]*cos(2*pi()*time[i]*f[2]+2*pi()*phi[2] )*exp(-a[2]*time[i]);
}
// priors
a~normal(0,2.2);
A[1]~normal(0,sqrt(dot_product(B1',B1)));
A[2]~normal(0,sqrt(dot_product(B2',B2)));
phi~uniform(0,1);
f~uniform(0,15);
// posterior
y~normal(B1+B2,sigma);
}
The output of the the 1 chain fit is stating close to the correct values of the parameters, whilst the 4 chain is much off. Utilizing some diagnostics from the 4 chain run I found the following:
none of the parameter chains are mixing - I tried to use ordered "insert parameter" [k] and also postive_ordered [k]
but that did not work. Does anyone have a suggestion as to what might be done about the problem? :)