I did get this to work though this is unlikely to work if you add another time series. I’m also doubtful on if this works on real data.
The issue with these models is that as you add more parameters there’s an increase in the number of combinations that would fit the data equally well. You must put fairly strong domain knowledge or other constraints onto the parameters to remove the label switching. However, at the point that these models become identifiable you almost already know what the values of the parameters are. So the utility of the model is low and this becomes especially relevant as the order of the parameter space increases.
So what did I do.
The least offensive changes were adding positive ordered constraints to A, a,phi, and f. Even with this the model was not identified due to f. You can see in the pairs plot that f could reasonably be [8.9, 9.0] or [8.0, 10.0],
where the second vector is the true parameter values. I added a repulsive prior (see Mixture Models) so that components of f near each other are “repulsed”. This is unlikely to work well if the dimension of f were > 2.
The output running on 4 chains is
> fit1
Inference for Stan model: sine_mix_forum.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
A[1] 0.52 0.00 0.03 0.45 0.50 0.52 0.54 0.59 1866 1
A[2] 5.23 0.00 0.06 5.12 5.19 5.23 5.27 5.35 2145 1
a[1] 1.05 0.00 0.04 0.95 1.03 1.05 1.08 1.12 1422 1
a[2] 1.08 0.00 0.03 1.03 1.07 1.08 1.10 1.13 2033 1
phi_raw 0.17 0.00 0.02 0.14 0.16 0.17 0.19 0.21 1479 1
dphi 0.32 0.00 0.02 0.29 0.31 0.32 0.33 0.35 1483 1
f_raw 8.01 0.00 0.00 8.00 8.00 8.01 8.01 8.01 5698 1
df 2.12 0.00 0.05 2.01 2.09 2.12 2.16 2.22 1631 1
sigma 0.15 0.00 0.01 0.12 0.14 0.14 0.15 0.17 2466 1
phi[1] 0.17 0.00 0.02 0.14 0.16 0.17 0.19 0.21 1479 1
phi[2] 0.49 0.00 0.00 0.49 0.49 0.49 0.50 0.50 4498 1
f[1] 8.01 0.00 0.00 8.00 8.00 8.01 8.01 8.01 5698 1
f[2] 10.13 0.00 0.05 10.02 10.09 10.13 10.16 10.23 1645 1
lp__ 135.03 0.06 2.31 129.69 133.65 135.42 136.73 138.44 1355 1
with a pairs plot of f now
Stan code
functions {
real repulsive_lpdf(vector mu, vector rho) {
int K = num_elements(mu);
matrix[K, K] S = diag_matrix(rep_vector(1, K));
matrix[K, K] L;
int c;
for (k1 in 1:(K - 1))
for (k2 in (k1 + 1):K){
c = K - (k2 - k1);
S[k1, k2] = exp(- squared_distance(mu[k1], mu[k2]) / (0.5 + rho[c]));
S[k2, k1] = S[k1, k2];
}
L = cholesky_decompose(S);
return 2 * sum(log(diagonal(L)));
}
}
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{
positive_ordered[k] A; // signal amplitude
positive_ordered[k] a; // signal decay constant
real<lower=0, upper=1> phi_raw;
real<lower=0, upper=1 - phi_raw> dphi;
real<lower=5, upper=15> f_raw;
real<lower=0, upper=15 - f_raw> df;
real<lower=0> sigma; // signal error
}
transformed parameters {
vector[k] phi;
vector[k] f;
// phi is an implied ordered vector with lower bound = 0 and upper bound = 1
phi[1] = phi_raw;
phi[2] = phi_raw + dphi;
// f is an implied ordered vector with lower bound = 5 and upper bound = 15
f[1] = f_raw;
f[2] = f_raw + df;
}
model {
// priors
a ~ normal(0, 2);
A ~ normal(0, 5);
sigma ~ normal(0, 0.2);
[f_raw, f_raw + df]' ~ repulsive([0.5, 0.5]');
// posterior
{
vector[N] B1 = A[1] * cos( 2 * pi() * time * f[2] + 2 * pi() * phi[1] ) .* exp( -a[1] * time );
vector[N] B2 = A[2] * cos( 2 * pi() * time * f[1] + 2 * pi() * phi[2] ) .* exp( -a[2] * time );
y ~ normal(B1 + B2, sigma);
}
}