I want to fit a variation of skellam which is diinflated at 0 and bounded to {-3,-2,-1,1,2,3}. For this reason in second function skellam_without_lpmf I give value -700 in log probability of these not allowed values in order to have these values a probability almost equal to 0. The probability mass function is converted to $\log(f(x)/f(-3)+f(-2)+f(-1)+f(1)+f(2)+f(3)). The problem is that when sampling finishes a message for divergent transitions is coming. In 3000 postwarmum iterations 2987 are the divergent transitions. Furthermore, if adapt delta is increased at 0.99, no improvement happens. Which can be the problem with this model? To mention that a corresponding model only "âdi-inflated at 0â, has not problem in convergence. The code is similar with this. Below is presented my stan code:
functions {
real skellam_lpmf(int k, real mu1, real mu2) {
real total;
real log_prob;
total = (- mu1 - mu2) + (log(mu1) - log(mu2)) * k / 2;
log_prob = total + log(modified_bessel_first_kind(k, 2 * sqrt(mu1*mu2)));
return log_prob;
}
real skellam_without_lpmf(int k, real mu1, real mu2) {
real log_prob_new;
real f_one;
real f_two;
real f_three;
real f_neg_one;
real f_neg_two;
real f_neg_three;
f_one=exp(skellam_lpmf(1|mu1,mu2));
f_two=exp(skellam_lpmf(2|mu1,mu2));
f_three=exp(skellam_lpmf(3|mu1,mu2));
f_neg_one=exp(skellam_lpmf(-1|mu1,mu2));
f_neg_two=exp(skellam_lpmf(-2|mu1,mu2));
f_neg_three=exp(skellam_lpmf(-3|mu1,mu2));
if (k > 3 ) {
log_prob_new=-700;
return log_prob_new;
} else if (k==0) {
log_prob_new=-700;
return log_prob_new;
} else if (k <(-3)){
log_prob_new=-700;
return log_prob_new;
} else {
log_prob_new=skellam_lpmf(k|mu1,mu2)-log(f_one+f_two+f_three+f_neg_one+f_neg_two+f_neg_three);
return log_prob_new;
}
}
}
data {
int <lower=1> n_games; //number of games 132
int <lower=1> n_teams; //number of teams 12
int <lower=1,upper=n_teams> home_team[n_games];
int <lower=1,upper=n_teams> away_team[n_games];
int <lower=0,upper=3> home_sets[n_games];//0-3 sets can have each team
int <lower=0,upper=3> away_sets[n_games];//0-3 sets can have each team
}
parameters {
real mu;
real home;
real attack_raw[n_teams - 1];
real defense_raw[n_teams - 1];
}
transformed parameters {
// Enforce sum-to-zero constraints
vector[n_teams] attack;
vector[n_teams] defense;
vector[n_games] lambda1;
vector[n_games] lambda2;
for (t in 1:(n_teams-1)) {
attack[t] = attack_raw[t];
defense[t] = defense_raw[t];
}
attack[n_teams] = -sum(attack_raw);
defense[n_teams] = -sum(defense_raw);
//likelihood-systematic component
for (g in 1:n_games) {
lambda1[g]=exp(mu+home+attack[home_team[g]]+defense[away_team[g]]);
lambda2[g]=exp(mu+attack[away_team[g]]+defense[home_team[g]]);
}
}
model {
int sets_diff[n_games];
//Priors
target+=normal_lpdf(mu|0,100);
target+=normal_lpdf(home|0,100);
target+=normal_lpdf(attack|0,100);
target+=normal_lpdf(defense|0,100);
//likelihood-systematic component
for (g in 1:n_games) {
sets_diff[g]=home_sets[g]-away_sets[g];
target+=skellam_without_lpmf(sets_diff[g]|lambda1[g],lambda2[g]);
}
}