add poisson_log_lcdf and poisson_log_lccdf functions

[davharris]

Summary:

[From @bob-carpenter: I edited the title to highlight the root cause. You can skip all the discussion; just missing the underlying functions at this point.]

I'm fitting a model with the poisson_log function (i.e. a Poisson distribution parameterized by log(lambda) instead of lambda).

I'm running into two problems:

• Large number of divergent transitions with a simple model (well over 50% of the post-warmup iterations)
• If I add truncation, I start getting errors about negative lambda values and my log(rate) parameter is squashed up against 0.

My best guess is that truncation flips a switch somewhere that leads to stan::math::poisson_log_log when it should go to to stan::math::poisson_log. I can't explain the divergent transitions, but I think it also might be related. Maybe an incorrect nonlinear transform?

Reproducible Steps:

compiled = stan_model(model_code =
"data {
  int<lower=1> N;
  int<lower=1> y[N];
}
parameters {
  real log_rate;
}
model {
  log_rate ~ normal(0.0, 1.0);
  for(i in 1:N){
    y[i] ~ poisson_log(log_rate) T[1, ];
  }
}
"
)

y = rpois(1000, 1)
y = y[y>0]
data = list(y = y, N = length(y))

m = sampling(compiled, data = data)

As noted above, the non-truncated version (i.e. the above code with T[1, ] removed on line 12) also has problems.

Current Output:

Sampling without truncation:

• There were 2768 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
• samples of log_rate approximately Gaussian and can span 0, depending on the distribution of y in data.

Sampling with truncation:

• If the model initializes, I get many iterations of Exception thrown at line 12: poisson_log: Rate parameter is -2.24632e-13, but must be >= 0! 1
• With an earlier (more complicated) version of the model that didn't initialize, I instead got [596] "Exception thrown at line 12: stan::math::poisson_log: Rate parameter is -0.989166, but must be >= 0!" Note the reference to stan::math::poisson_log instead of stan::math::poisson_log_log
• There were 1065 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help

Expected Output:

MCMC samples

Current Version:

rstan 2.11.1

[davharris]

I think this may actually be two problems, one with poisson_log and one with truncated Poissons.

I get absurdly negative values of log_rate if I swap in the non-log poisson distribution as in y[i] ~ poisson(exp(log_rate)) T[1, ]; (log_rate < -250). So I don't think it's just a problem with poisson_log.

Either that or I'm deeply confused, in which case I apologize for taking up your time with this!

[rtrangucci]

I can't run this example right now, but the parameter sigma will cause problems if it isn't integrated into your model or if it doesn't have a proper prior on it. Try removing sigma as a parameter and see if that fixes things.

I agree that the error messages are odd, I'll look into those as well when I'm back in front of a computer.

Rob

On Aug 16, 2016, at 8:09 PM, David J. Harris notifications@github.com wrote:

I think this may actually be two problems, one with poisson_log and one with truncated Poissons.

I get absurdly negative values of log_rate if I swap out in the non-log poisson distribution as in y[i] ~ poisson(exp(log_rate)) T[1, ]; (<-250). So I don't think it's just a problem with poisson_log.

Either that or I'm deeply confused, in which case I apologize for taking up your time with this!

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub, or mute the thread.

[rtrangucci]

Whoops, looks like you edited the code. I can confirm that using the truncation with the stan::math::poisson_log_log does seem to lead to sampling problems. My guess is that there isn't a stan::math::poisson_cdf_log_log, but that T[1,] for some reason gets expanded to an expression involving stan::math::poisson_cdf_log. @davharris Thank you for catching this!

[bob-carpenter]

@rtrangucci Thanks for taking a look, but if that was a problem, the C++ for the translated Stan program wouldn't compile.

@davharris Does the equivalent program run if you use the basic Poisson?

[bob-carpenter]

I verified there's definitely a bug with Poisson truncation, even without the log scale. No idea why our unit tests didn't catch this.

data {
  int<lower=1> N;
  int<lower=1> y[N];
}
parameters {
  real<lower=0> rate;
}
model {
  rate ~ normal(0.0, 1.0);
  for(i in 1:N) {
    y[i] ~ poisson(rate) T[1, ];
  }
}

Fit I get (with original N = 100) is:

> fit <- stan("poistrunc.stan", chains=1, iter=500, data=data)
...
> print(fit, digits=10, "rate", probs=c())
Inference for Stan model: poistrunc.
1 chains, each with iter=500; warmup=250; thin=1; 
post-warmup draws per chain=250, total post-warmup draws=250.

     mean se_mean sd n_eff Rhat
rate    0       0  0   250  NaN

[bob-carpenter]

For a workaround, you can use the version recommended in the manual for truncating Poisson at 1:

data {
  int<lower=1> N;
  int<lower=1> y[N];
}
parameters {
  real<lower=0> rate;
}
model {
  rate ~ normal(0.0, 1.0);
  for(i in 1:N) {
    y[i] ~ poisson(rate);
    target += -log1m(poisson_lpmf(0 | rate));
  }
}

which gives me the fit

     mean se_mean   sd n_eff Rhat
rate 1.07       0 0.04   191    1

This also works with a log-scale parameterization (not quite same prior, but it's dominated by data here anyway):

data {
  int<lower=1> N;
  int<lower=1> y[N];
}
parameters {
  real log_rate;
}
model {
  log_rate ~ normal(0.0, 1.0);
  for(i in 1:N) {
    y[i] ~ poisson_log(log_rate);
    target += -log1m(poisson_log_lpmf(0 | log_rate));
  }
}
generated quantities {
  real<lower=0> rate;
  rate = exp(log_rate);
}

the fit for which is

         mean se_mean   sd n_eff Rhat
rate     1.07       0 0.03   250    1
log_rate 0.07       0 0.03   250    1

[davharris]

Thanks! Also, let me know when your testing infrastructure is ready. I'd be happy to contribute a few tests.

Dave

On Aug 17, 2016, at 6:53 AM, Bob Carpenter notifications@github.com wrote:

For a workaround, you can use the version recommended in the manual for truncating Poisson at 1:

data {
int<lower=1> N;
int<lower=1> y[N];
}
parameters {
real<lower=0> rate;
}
model {
rate ~ normal(0.0, 1.0);
for(i in 1:N) {
y[i] ~ poisson(rate);
target += -log1m(poisson_lpmf(0 | rate));
}
}
which gives me the fit

 mean se_mean   sd n_eff Rhat

rate 1.07 0 0.04 191 1
This also works with a log-scale parameterization (not quite same prior, but it's dominated by data here anyway):

data {
int<lower=1> N;
int<lower=1> y[N];
}
parameters {
real log_rate;
}
model {
log_rate ~ normal(0.0, 1.0);
for(i in 1:N) {
y[i] ~ poisson_log(log_rate);
target += -log1m(poisson_log_lpmf(0 | log_rate));
}
}
generated quantities {
real<lower=0> rate;
rate = exp(log_rate);
}
the fit for which is

     mean se_mean   sd n_eff Rhat

rate 1.07 0 0.03 250 1
log_rate 0.07 0 0.03 250 1
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

[bob-carpenter]

Not sure what you mean by testing infrastructure. All of our repos have pretty extensive (though obviously not complete!) unit tests in place.

[davharris]

Hey Bob, sorry—I didn’t mean to imply that the project didn’t have any tests!

I was referring to the relative difficulty of writing additional integration tests until the “infamous refactor”: https://twitter.com/betanalpha/status/757741262468448256 https://twitter.com/betanalpha/status/757741262468448256

On Aug 17, 2016, at 9:23 AM, Bob Carpenter notifications@github.com wrote:

Not sure what you mean by testing infrastructure. All of our
repos have pretty extensive (though obviously not complete!)
unit tests in place.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub #2015 (comment), or mute the thread https://github.com/notifications/unsubscribe-auth/AAzdCfiIm1cRqwqERNrTT-YFgpJncYgnks5qgwtFgaJpZM4Jl9nP.

[syclik]

I don't think any of us read that as you implying that we don't have tests.
I had the same confusion as Bob, though -- not sure what you mean by
testing infrastructure.

Regarding Michael's tweet -- he's referring to a set of statistical tests
that are waiting for changes to the Stan infrastructure. This doesn't need
to wait for anything. Almost everything in the math library can be tested
with deterministic unit tests. We have plenty of unit tests to use as
examples if you're wanting to add more tests.

Daniel

On Wed, Aug 17, 2016 at 9:56 AM, David J. Harris notifications@github.com
wrote:

Hey Bob, sorry—I didn’t mean to imply that the project didn’t have any
tests!

I was referring to the relative difficulty of writing additional
integration tests until the “infamous refactor”: https://twitter.com/
betanalpha/status/757741262468448256 <https://twitter.com/
betanalpha/status/757741262468448256>

On Aug 17, 2016, at 9:23 AM, Bob Carpenter notifications@github.com
wrote:

Not sure what you mean by testing infrastructure. All of our
repos have pretty extensive (though obviously not complete!)
unit tests in place.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub <
https://github.com/stan-dev/stan/issues/2015#issuecomment-240409068>, or
mute the thread <https://github.com/notifications/unsubscribe-
auth/AAzdCfiIm1cRqwqERNrTT-YFgpJncYgnks5qgwtFgaJpZM4Jl9nP>.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
#2015 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AAZ_F42NYDJ3zILkgIxpwKSC37KEI0odks5qgxL9gaJpZM4Jl9nP
.

[davharris]

Oh, okay. Sorry for the confusion. Is there a to-do list for tests that need to be written? If so, I could start there.

Dave

On Aug 17, 2016, at 10:10 AM, Daniel Lee notifications@github.com wrote:

I don't think any of us read that as you implying that we don't have tests.
I had the same confusion as Bob, though -- not sure what you mean by
testing infrastructure.

Regarding Michael's tweet -- he's referring to a set of statistical tests
that are waiting for changes to the Stan infrastructure. This doesn't need
to wait for anything. Almost everything in the math library can be tested
with deterministic unit tests. We have plenty of unit tests to use as
examples if you're wanting to add more tests.

Daniel

On Wed, Aug 17, 2016 at 9:56 AM, David J. Harris notifications@github.com
wrote:

Hey Bob, sorry—I didn’t mean to imply that the project didn’t have any
tests!

I was referring to the relative difficulty of writing additional
integration tests until the “infamous refactor”: https://twitter.com/
betanalpha/status/757741262468448256 <https://twitter.com/
betanalpha/status/757741262468448256>

On Aug 17, 2016, at 9:23 AM, Bob Carpenter notifications@github.com
wrote:

Not sure what you mean by testing infrastructure. All of our
repos have pretty extensive (though obviously not complete!)
unit tests in place.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub <
https://github.com/stan-dev/stan/issues/2015#issuecomment-240409068>, or
mute the thread <https://github.com/notifications/unsubscribe-
auth/AAzdCfiIm1cRqwqERNrTT-YFgpJncYgnks5qgwtFgaJpZM4Jl9nP>.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
#2015 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AAZ_F42NYDJ3zILkgIxpwKSC37KEI0odks5qgxL9gaJpZM4Jl9nP
.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.