// // This Stan program implements a poisson binomial // regression where the probability vector is // parameterized as a probit. The object is to // recover the probit coefficients // I use a poisson_binomial_lpmf function taken from the // Stan forums, see: https://discourse.mc-stan.org/t/poisson-binomial-distribution-any-existing-stan-implementation/ functions { /** * Return the Poisson-binomial log probability mass for the specified * count y and vector of probabilities theta. The result is the log * probability of y successes in N = size(theta) independent * trials with probabilities of success theta[1], ..., theta[N]. * * See: https://en.wikipedia.org/wiki/Poisson_binomial_distribution * * @param y number of successes * @param theta vector of probabilities * @return Poisson-binomial log probability mass */ real poisson_binomial_lpmf(int y, vector theta) { int N = rows(theta); matrix[N + 1, N + 1] alpha; vector[N] log_theta = log(theta); vector[N] log1m_theta = log1m(theta); if (y < 0 || y > N) reject("poisson binomial variate y out of range, y = ", y, " N = ", N); for (n in 1:N) if (theta[n] < 0 || theta[n] > 1) reject("poisson binomial parameter out of range,", " theta[", n, "] =", theta[n]); if (N == 0) return y == 0 ? 0 : negative_infinity(); // dynamic programming with cells // alpha[n + 1, tot + 1] = log prob of tot successes in first n trials alpha[1, 1] = 0; for (n in 1:N) { // tot = 0 alpha[n + 1, 1] = alpha[n, 1] + log1m_theta[n]; // 0 < tot < n for (tot in 1:min(y, n - 1)) alpha[n + 1, tot + 1] = log_sum_exp(alpha[n, tot] + log_theta[n], alpha[n, tot + 1] + log1m_theta[n]); // tot = n if (n > y) continue; alpha[n + 1, n + 1] = alpha[n, n] + log_theta[n]; } return alpha[N + 1, y + 1]; } } data { int N; // N obs int K; // N covariates int G; // number of groups int y[G]; // matrix[N, K] X; // covariates int gsize[G]; // vector of group sizes } // The parameters accepted by the model. We // estimate an intercept and the covariates delta parameters { real mu; vector[K] delta; } // Make the probability param transformed parameters { vector[N] p; // probability p = Phi(mu + X * delta); } // The model to be estimated. We model the output // 'y' to be poisson-binomially distributed with // individual probabilities p that are a function of // the parameters delta and mu model { int pos; pos = 1; mu ~ std_normal(); delta ~ normal(0, 5); for (g in 1:G) { target += poisson_binomial_lpmf(y[g] | segment(p, pos, gsize[g])); pos = pos + gsize[g]; } }