Thanks for the detailed response!
I was able to get the mixture to work without the overinflation by changing the code a bit. You’ll see I changed the specification of the repulsive distribution. I’m only adding the pairwise comparisons to lp. That means I have to build the full matrix for S. I’ve let rho vary by each pairwise comparison but it seems to like 0.5 for all the pairs.
But this breaksdown if I run 3 chains. I’m guessing this is the sensitivity to initial values.
Full code and results for 1 chain looks good:
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 = 0;
for (k1 in 1:(K - 1))
for (k2 in (k1 + 1):K){
c += 1;
S[k1, k2] = log(1 - exp(- squared_distance(mu[k1], mu[k2]) / rho[c]));
S[k2, k1] = S[k1, k2];
}
L = cholesky_decompose(S);
return 2 * sum(log(diagonal(L)));
}
}
data {
int<lower=1> K;
int<lower=1> N;
real y[N];
}
parameters {
ordered[K] mu;
vector<lower=0>[choose(K, 2)] rho;
real<lower=0> sigma[K];
simplex[K] lambda;
}
model {
// Prior model
mu ~ normal(0, 5);
sigma ~ std_normal();
lambda ~ dirichlet(rep_vector(3, K));
rho ~ gamma(1, 2);
mu ~ repulsive(rho);
// Observational model
for (n in 1:N) {
real comp_lpdf[K];
for (k in 1:K) {
comp_lpdf[k] = log(lambda[k]) + normal_lpdf(y[n] | mu[k], sigma[k]);
}
target += log_sum_exp(comp_lpdf);
}
}
Inference for Stan model: beta_alpha_sean.
1 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=1000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
mu[1] -5.05 0.00 0.07 -5.17 -5.09 -5.05 -5.00 -4.91 643 1
mu[2] -0.98 0.00 0.03 -1.04 -1.00 -0.98 -0.97 -0.93 1489 1
mu[3] 4.04 0.00 0.01 4.01 4.03 4.04 4.05 4.06 999 1
mu[4] 6.99 0.00 0.03 6.93 6.97 6.99 7.01 7.05 1197 1
rho[1] 0.49 0.01 0.49 0.01 0.13 0.34 0.71 1.81 1449 1
rho[2] 0.49 0.01 0.46 0.01 0.15 0.35 0.68 1.66 1519 1
rho[3] 0.51 0.01 0.51 0.01 0.14 0.36 0.73 1.94 1437 1
rho[4] 0.48 0.01 0.48 0.01 0.13 0.31 0.67 1.77 1778 1
rho[5] 0.49 0.01 0.49 0.02 0.14 0.33 0.71 1.71 1711 1
rho[6] 0.51 0.01 0.50 0.01 0.15 0.36 0.70 1.94 1387 1
sigma[1] 0.95 0.00 0.05 0.85 0.91 0.94 0.98 1.05 1588 1
sigma[2] 0.55 0.00 0.02 0.51 0.53 0.54 0.56 0.59 1125 1
sigma[3] 0.24 0.00 0.01 0.22 0.24 0.24 0.25 0.26 1380 1
sigma[4] 0.32 0.00 0.02 0.27 0.30 0.32 0.33 0.37 1394 1
lambda[1] 0.18 0.00 0.01 0.16 0.18 0.18 0.19 0.21 1313 1
lambda[2] 0.40 0.00 0.02 0.37 0.38 0.40 0.41 0.43 1524 1
lambda[3] 0.32 0.00 0.01 0.29 0.31 0.32 0.33 0.35 1653 1
lambda[4] 0.10 0.00 0.01 0.08 0.09 0.10 0.11 0.12 1520 1
lp__ -1899.03 0.17 3.04 -1905.97 -1900.81 -1898.75 -1896.88 -1893.98 314 1
Samples were drawn using NUTS(diag_e) at Mon Sep 7 08:00:32 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
singular_fit <- stan(file='beta_alpha_sean.stan', data=input_data,
chains = 1, iter = 2000, warmup = 1000,
seed = 483892929, refresh = 500)
singular_fit
Here’s the plot
params1 <- as.data.frame(extract(singular_fit, permuted=FALSE)[,1,])
c_light_trans <- c("#DCBCBCBF")
c_light_highlight_trans <- c("#C79999BF")
c_mid_trans <- c("#B97C7CBF")
c_mid_highlight_trans <- c("#A25050BF")
c_dark_trans <- c("#8F2727BF")
c_dark_highlight_trans <- c("#7C0000BF")
par(mar = c(4, 4, 0.5, 0.5))
plot(params1$"mu[1]", params1$"mu[2]", col=c_dark_highlight_trans, pch=16, cex=0.8,
xlab="pair_1", xlim=c(-6, 8), ylab="pair_2", ylim=c(-2, 8))
points(params1$"mu[1]", params1$"mu[3]", col="blue", pch=16, cex=0.8)
points(params1$"mu[1]", params1$"mu[4]", col="green", pch=16, cex=0.8)
points(params1$"mu[2]", params1$"mu[3]", col="orange", pch=16, cex=0.8)
points(params1$"mu[2]", params1$"mu[4]", col="purple", pch=16, cex=0.8)
points(params1$"mu[3]", params1$"mu[4]", col=c_light_trans, pch=16, cex=0.8)
Adding 3 chains and you’ll see the separation break down in chains 2 and 3
