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 N; int y[N]; } parameters { real 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 N; int 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 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 https://github.com/stan-dev/stan/issues/2015#issuecomment-240409068, 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 https://github.com/stan-dev/stan/issues/2015#issuecomment-240418724, 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 https://github.com/stan-dev/stan/issues/2015#issuecomment-240418724, 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.

[bob-carpenter]

We don't have a separate list of pending tests. We just write them when we realize they're missing.

You can start an issue if there's something specific you're thinking about. Now we're off topic on this issue, so let's move it either to the stan-dev mailing list or to a new issue. Thanks.

[bob-carpenter]

I think this may be related to the incomplete gamma function. @betanalpha, mind taking a look?

[bob-carpenter]

Back to testing, I just did a sweep over our code during some refactoring and realized there are a whole lot of missing tests in the math library.

I think for the Poisson CDF we're just going to need to revert to our previous implementation because this one's not working and I don't think Michael's going to have time to fix it before our next release.

[betanalpha]

a) stan_rdump it is a wonderful function, and it makes debugging models by non RStan users much, much easier. Please provide data files and not R scripts, or (devs) assign problems to RStan users.

b) The original model works absolutely fine without truncation in CmdStan. I cannot reproduce the problems there at all.

c) Following @bob-carpenter I simplified the problem by looking at truncation on the linear scale first. There are no divergences or warnings about bounds violations, but the rate does indeed get driven to zero which appears to be an issue.

d) The problem is not in poisson_cdf but rather poisson_ccdf_log. The generated code seems to be doing the right thing, and I will look into the poisson_ccdf_log implementation to see if I can find a problem.

[betanalpha]

I rederived the `poisson_cdf_logandpoisson_ccdf_log`` derivatives and the implementations match the math.

[betanalpha]

@bob-carpenter Your suggested fix used the Poisson mass function, which works only because the cut off is exactly one. Can you get the model to run explicitly using any of the CDF functions? I rerun an old truncation model that used the CDFs directly and it worked fine, but we might as well try to verify it on this model.

[betanalpha]

Both

data {
  int<lower=1> N;
  int<lower=1> y[N];
}

parameters {
  real<lower=0> rate;
}

model {
  rate ~ normal(0.0, 1.0);
  for(n in 1:N) {
    y[n] ~ poisson(rate);
    target += - log(1 - poisson_cdf(1, rate));
  }
}

and

data {
  int<lower=1> N;
  int<lower=1> y[N];
}

parameters {
  real<lower=0> rate;
}

model {
  rate ~ normal(0.0, 1.0);
  for(n in 1:N) {
    y[n] ~ poisson(rate);
    target += - poisson_lccdf(1 | rate);
  }
}

give the same results. So either we're thinking about the model incorrectly or all of the Poisson cumulatives are broken, which would be weird given that other models using them seem to work fine. Are we missing something obvious about why this model would be fitting weird?

Also, poisson_ccdf apparently isn't exposed?

[bob-carpenter]

Looking at Michael's code, I realize my code's wrong and shouldn't have log1m(poisson_lcdf(...))---it should be log1m(exp(poisson_lcdf(...))) or more succintly and stably poisson_lccdf(...).

[bob-carpenter]

I'm not sure where that leaves us. I couldn't get the Stan program with truncation to work properly and can't tell from these messages if you got it to work.

There aren't linear ccdf functions in Stan's language --- if they're in the math lib, they'd be easy to expose with a careful day of writing test files, but that should be a separate issue.

[syclik]

@bob-carpenter, I found it! The generated code is wrong. I haven't dug into the cause, but I've definite found an instance.

This model generates the correct code. It's with the rate parameter in the natural scale:

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

The relevant part of the generated code looks like:

for (int i = 1; i <= N; ++i) {
  current_statement_begin__ = 11;
  lp_accum__.add(poisson_log<propto__>(get_base1(y,i,"y",1), rate));
  if (get_base1(y,i,"y",1) < 1) lp_accum__.add(-std::numeric_limits<double>::infinity());
  else lp_accum__.add(-poisson_ccdf_log(1, rate));

This model on the log scale fails, and it's probably due to missing cdf implementations with the parameter on the log scale:

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, ];
  }
}

The current, buggy generated code looks like:

for (int i = 1; i <= N; ++i) {
  current_statement_begin__ = 11;
  lp_accum__.add(poisson_log_log<propto__>(get_base1(y,i,"y",1), log_rate));
  if (get_base1(y,i,"y",1) < 1) lp_accum__.add(-std::numeric_limits<double>::infinity());
  else lp_accum__.add(-poisson_log(1, log_rate));
}

What it should look like is:

for (int i = 1; i <= N; ++i) {
  current_statement_begin__ = 11;
  lp_accum__.add(poisson_log_log<propto__>(get_base1(y,i,"y",1), log_rate));
  if (get_base1(y,i,"y",1) < 1) lp_accum__.add(-std::numeric_limits<double>::infinity());
  else lp_accum__.add(-poisson_ccdf_log(1, exp(log_rate)));
}

(The key difference here is: poisson_log(1, log_rate) vs poisson_ccdf_log(1, exp(log_rate))... or if we implemented more functions, poisson_log_ccdf_log(1, log_rate))

I think we need to add a poisson_log_ccdf_log(). In this case, we should not have allowed this to compile since there's no matching cdf function.

[syclik]

Bob, I just created a test on branch: feature/issue-2015-truncation.

To run the test:

./runTests.py src/test/unit/lang/parser/trunc_test.cpp