# Sampling a discrete parameter from a bernoulli prior

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (stan) with clear spacing and indentation. For example, use


//
data {
int<lower=0> N; //number of observations
int T; //number of timepoint
int<lower=0> Ngroups; //number of unique subjects
int subjects[N]; //60
vector[N] treatment;
int familyRichness[N];
real IgA[N];
matrix[N, T] X_t; //intercept per time point
matrix[N, T] X_trt; // treatment effect per time point
matrix[N, T] X_fam; //family richness per timepoint

}

// The parameters accepted by the model. Our model
parameters {
vector[T] mux; //intercept per time point for poisson model
vector[T] muy; //intercept per time point for normal model
vector[T] beta; // treatment effect per time point for normal
vector[T] alpha; // treatment effect for poisson model
vector phi;
vector[T] delta;
real<lower=0> p4;
real<lower=0> sigma;
real<lower=0> sigma_iga; // random intercept variance for iga
real<lower=0> sigma_rich; // random intercept variance for richness
cholesky_factor_corr L_iga_rich;
matrix[2, Ngroups] z_iga_rich; //random effect from N(0, 1) to post multiply

}

transformed parameters {

real<lower=0> lambda[N];
vector[T] phis;
real mu[N];
matrix[2, Ngroups] ran_effs;
vector L_sigma;
vector[T] gamma;

L_sigma = sigma_iga;
L_sigma = sigma_rich;

phis = 0;
phis = 0;
phis = phi;
phis = phi;

for(i in 1:4)
gamma[i] = phi[i]*delta[i];

ran_effs = diag_pre_multiply(L_sigma, L_iga_rich) * z_iga_rich;

for(i in 1:N) {
lambda[i] = exp(X_t[i, ] * mux + X_trt[i, ] * alpha + ran_effs[2, subjects[i]]);//poisson
mu[i] =  X_t[i, ] * muy + X_trt[i] * beta + X_fam[i, ]* gamma + ran_effs[1, subjects[i]]; //normal
}

}

model {

int p1;
int p2;

L_iga_rich ~ lkj_corr_cholesky(2.0);
to_vector(z_iga_rich) ~ normal(0, 1);
sigma ~ cauchy(0, 5);
sigma_iga ~ cauchy(0, 5);
sigma_rich ~ cauchy(0, 5);
mux ~ cauchy(0, 10);
muy ~ cauchy(0, 10);
alpha ~ normal(0, 1);
beta ~ normal(0, 1);
delta ~ normal(0, 1);

phi ~ bernoulli(0.5);
phi ~ bernoulli(p4);
p4 ~ uniform(0,1);

for(i in 1:N) {
familyRichness[i] ~ poisson(lambda[i]);
IgA[i] ~ normal(mu[i], sigma);
}
}

generated quantities {

real rho;
matrix[2, 2] Sigma;

//compute sigma
Sigma = diag_pre_multiply(L_sigma, L_iga_rich * L_iga_rich');

//compute rho
rho = Sigma[1, 2] / (pow(Sigma[1, 1], 0.5) * pow(Sigma[2, 2], 0.5));
}
`

X_{ij} \sim dpois(\lambda_{ij}), log(\lambda_{ij}) = \mu_X[j] + b1_i + \alpha[j] Z_i and
Y_{ij} \sim dnorm(\mu_{ij},\sigma), \mu_{ij} = \mu_Y[j] + b2_i + \beta[j] Z_i + \phi[j]*\delta[j] X_{ij}.
The \phi[j]'s are the challenge here. I fixed its first two components to be 0, while the last two components are to be simulated from the Bernoulli distribution. I understand that stan does not allow for declaration of discrete parameters and one has to marginalize over discrete parameters, but how do I do this in my case? Thanks in anticipation.

Hi Olajumoke! Welcome to the Stan forum! :)

Sorry it took us so long to respond to your question.The model that you posted is not so trivial it seems. Maybe you could give a bit more context, for example about what those \phi 's are doing. Also, it would be great if you could post (a subset) of your data or data that could look like your real data for us to be able to run your model.

I think then it will become a bit easier to talk about marginalizing those Bernoulli r.v.

Cheers!
Max

2 Likes