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]