There has been some talk about dealing with either overdispersed or underdispersed count data. I have implemented the Conway-Maxwell generalization of the Poisson PMF, which is attached. It uses a brute force but numerically stable approach to calculating the normalizing constant, which will be slow if you parameterize the location parameter.
CMP.stan (2.1 KB)
3 Likes
target += sum_y * log_lambda - nu + sum_lfac - N * log_C(log_lambda, nu);
shouldn’t that be:
target += sum_y * log_lambda - nu * sum_lfac - N * log_C(log_lambda, nu);
Can nu = 0 happen? (Except by setting nu = 0)
1 Like
Thanks. I edited the OP to multiply rather than add in the target statement. In this case, the only way nu can be zero is with underflow.
Should we be putting this distribution in the math lib?
This
real sum_lfac = 0;
for (n in 1:N) sum_lfac += lgamma(y[n] + 1);
makes me realize we should really have a sum-of-lgamma-plus-1 function because it’s the log product of factorials and it’s so common for combinatoric normalization.
The above could be written as
real sum_lfac = sum(lgamma(to_vector(y) + 1));
but that’d actually be less efficient than the way Ben had it because of the extra copy in to_vector().
It could. I don’t know a clever way to do the derivative of the normalizing constant, so it would only help so much.
The CMP uses a clever variant of the taylor expension of the exp function. What could help with lgamma is to precalculate values of lgamma up to max(y) and store those values. max(y) is usually small, if large it may be worth think about approx. with some continues distribution, eg. normal.
It would make sense to precompute lgamma for discrete cases, but you have to be careful about memory paging and indexing to ensure fetch times and branching don’t dominate arithmetic times. It’d probably be a win for something intensive like lgamma that gets called a lot if there aren’t too many possible integer values.
Dear Ben
Sorry to bother but we are having problems implementing the CMP.stan code.
The problem for me is that I am trying to use that code as a custom family function for Stan via the brms package.
The relevant stan_var for the brms custom family function based on your code is at end of message along with the custom_family block for brms
The only difference being that I replace “lambda” with”mu” as that is a requirement of brms.
Also the reject() has been made an explicit string as also required. I thank Paul Buerkner for that and other important advice.
OK, the model compiles. It starts sampling but then freezes any desktop being used within a minute or two of sampling. Completely freezes my desktop.
No further sampling activity for an hour or more and locked out desktop for any task - like a kernel panic.
So two queries.
(1) have you been successful in running a CMP regression model using Stan and your CMP.stan using a real world data example?
(2) do you think that there might been something wrong with the log_C function and specifically the while statement: while(term>-745){ etc.
Regarding (2), it appears that perhaps the while loop has no appropriate stopping condition except the underflow check. Could that be the problem?
I am trying to coax Paul to make a CMP family a native brms family but he is understandably reticent to do so until we get the code working properly with a good real world data example.
The real world example I am using to test any CMP model is sourced from:
## Chanialidis C, Evers L, Neocleous T, Nobile A (2018)
## Efficient Bayesian inference for COM-Poisson regression models
## Statistics and Computing 28(3): 595–608
> require(Rchoice)
> data(Articles)
>
> # from Chanialidis et al (2018)
> # Focusing only on the students with at least one publication
> phdpublish<-subset(Articles, art>0)
> phdpublish<-transform(phdpublish, art=art-1)
>
> # Standardise all non-binary covariates
> phdpublish<-cbind(phdpublish[,c(1,2,3)],scale(phdpublish[,-c(1,2,3)],center=TRUE,scale=TRUE))
OK, so the CMP model is then: art~fem+mar+kid5+phd+ment, model likelihood=compois
So I would be most grateful for any assistance.
best regards
Milani
MacOS High Sierra 10.13.4
R 3.5.0
rstan_2.17.3
brms_2.2.2
###############
compois <- custom_family(
"CMP", dpars = c("mu", "nu"),
links = c("log", "identity"), lb = c(NA, 0),
type = "int")
###############
stan_funs <- "
real log_C(real log_mu, real nu) {
real log_C = log1p_exp(log_mu);
real lfac = 0;
real term = 0;
real k = 2;
if (nu < 0) reject(\"nu must be non-negative\");
if (nu == 0) {
if (log_mu >= 0) reject(\"log_mu must be < 0\");
return -log1m_exp(log_mu);
}
if (nu == 1) return exp(log_mu);
if (nu == positive_infinity()) return log_C;
while (term > -745) { // exp(-745) underflows
lfac += log(k);
term = k * log_mu - nu * lfac;
log_C = log_sum_exp(log_C, term);
k += 1;
}
return log_C;
}
real CMP_lpmf(int y, real mu, real nu) {
if (nu == 0) return y * log(mu) + log1m(mu);
if (nu == 1) return poisson_lpmf(y | mu);
if (nu == positive_infinity())
return bernoulli_lpmf(y | mu / (1 + mu));
{
real log_mu = log(mu);
return y * log_mu - nu * lgamma(y + 1) - log_C(log_mu, nu);
}
}
real CMP_log_lpmf(int y, real log_mu, real nu) {
if (nu == 0) return y * log_mu + log1m_exp(log_mu);
if (nu == 1) return poisson_log_lpmf(y | log_mu);
if (nu == positive_infinity()) return bernoulli_logit_lpmf(y | log_mu);
return y * log_mu - nu * lgamma(y + 1) - log_C(log_mu, nu);
}
int CMP_log_rng(real log_mu, real nu) {
int y = 0;
real log_c = log_C(log_mu, nu);
real u = uniform_rng(0,1);
real CDF = exp(-log_c);
real lfac = 0;
if (CDF >= u) return 0;
while (CDF < u) {
y += 1;
lfac += log(y);
CDF += exp(y * log_mu - nu * lfac - log_c);
}
return y - 1;
}
int CMP_rng(real mu, real nu) {
return CMP_log_rng(log(mu), nu);
}
"
[edit: escaped code]
An experimental version of the COM Poisson distribution is now available in the github version of brms. See https://github.com/paul-buerkner/brms/issues/607 for details.
2 Likes
@Bob_Carpenter any chance this might make it into a future version of Stan? Just found our way here from the goog looking for this exact thing. Thank you @bgoodri for having this here in the forums!