I’m trying to fit a mixture of Generalized Poisson’s to simulated data. I was able to define a LPMF for a GP and fit a non-mixture dataset fine:
functions {
real genpoiss_lpmf(int y, real lambda1, real lambda2) {
return log(lambda1) +
(y - 1) * log(lambda1 + y * lambda2) -
lgamma(y + 1) -
y * lambda2 - lambda1;
}
}
data {
// data
int N;
int y[N];
vector[N] count;
// priors
real<lower=0> mu0;
real<lower=0> sigma0;
real<lower=0> phi0;
}
parameters {
real<lower=0> mu;
real<lower=0> sigma2;
}
transformed parameters {
real<lower=0> lambda1;
real lambda2;
lambda1 = mu * sqrt(mu / sigma2);
lambda2 = 1 - sqrt(mu / sigma2);
}
model {
real contributions;
mu ~ normal(mu0, sigma0);
sigma2 ~ normal(0, phi0);
for (n in 1:N) {
contributions = genpoiss_lpmf(width[n] | lambda1, lambda2);
target += contributions * count[n];
}
}
However, even an “easy” mixture model (large separation between means, only two components) does not recover the parameters well - most of the transitions after warmup are divergent despite increase the adapt_delta. Additionally, initializing stan with the true parameter values throws
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can't start sampling from this initial value.
until failing after 100 attempts, so I feel like I must have an error in my stan model.
Below is the code to simulate the data in R as well as the stan code to fit the model.
functions {
real genpoiss_lpmf(int y, real lambda1, real lambda2) {
return log(lambda1) +
(y - 1) * log(lambda1 + y * lambda2) -
lgamma(y + 1) -
y * lambda2 - lambda1;
}
}
data {
// data
int N; // number observations
int K; // number clusters
int y[N]; // fragment size
vector[N] count; // count of fragments of size n
// priors
vector<lower=0>[K] mu0;
real<lower=0> sigma0;
real<lower=0> phi0;
real<lower=0> alpha;
}
parameters {
vector<lower=0>[K] mu;
vector<lower=0>[K] sigma2;
simplex[K] theta;
}
transformed parameters {
vector<lower=0>[K] lambda1;
vector[K] lambda2;
lambda1 = mu .* sqrt(mu ./ sigma2);
lambda2 = 1 - sqrt(mu ./ sigma2);
}
model {
real contributions[K];
// prior
for (k in 1:K) {
mu ~ normal(mu0[k], sigma0);
}
sigma2 ~ normal(0, phi0);
theta ~ dirichlet(rep_vector(alpha / K, K));
for (n in 1:N) {
for (k in 1:K) {
contributions[k] = log(theta[k]) + genpoiss_lpmf(y[n] | lambda1[k], lambda2[k]);
}
target += log_sum_exp(contributions) * count[n];
}
}
library(rstan)
library(HMMpa)
# simulate data
mu <- c(20, 50)
sigma2 <- c(3, 4)
theta <- c(0.4, 0.6)
lambda1 <- mu * sqrt(mu / sigma2)
lambda2 <- 1 - sqrt(mu / sigma2)
set.seed(1)
n <- 1000
z <- sample(1:length(theta), n, replace=TRUE, prob=theta)
y <- rep(0, n)
for (i in 1:n) {
y[i] <- rgenpois(1, lambda1[z[i]], lambda2[z[i]])
}
tmp <- table(y)
y.value <- as.numeric(names(tmp))
y.count <- as.numeric(tmp)
plot(y.value, y.count)
# fit model in stan
model <- stan_model("stan/gpd-mixture.stan", model_name="gpd-mixture")
stan_data <- list(N=length(y.value),
K=length(mu),
y=y.value,
count=y.count,
mu0=mu,
sigma0=1e-2,
phi0=2,
alpha=1)
initf <- function() {
list(mu=mu, sigma2=sigma2, theta=theta)
}
results <- sampling(model, data=stan_data,
chains=4, cores=4,
init=initf,
#control=list(adapt_delta=0.999, max_treedepth=15),
seed=428, iter=2000)