search

alashworth / test-issue-import

0 stars 0 forks source link

add poisson_log_lcdf and poisson_log_lccdf functions #107

Open alashworth opened 5 years ago

alashworth commented 5 years ago

Issue by davharris Tuesday Aug 16, 2016 at 23:48 GMT Originally opened as https://github.com/stan-dev/stan/issues/2015

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:

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:

Sampling with truncation:

MCMC samples

Current Version:

rstan 2.11.1

alashworth commented 5 years ago

Comment by davharris Wednesday Aug 17, 2016 at 00:09 GMT

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!

alashworth commented 5 years ago

Comment by rtrangucci Wednesday Aug 17, 2016 at 02:43 GMT

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.

alashworth commented 5 years ago

Comment by rtrangucci Wednesday Aug 17, 2016 at 03:12 GMT

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!

alashworth commented 5 years ago

Comment by bob-carpenter Wednesday Aug 17, 2016 at 10:10 GMT

@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?

alashworth commented 5 years ago

Comment by bob-carpenter Wednesday Aug 17, 2016 at 10:44 GMT

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
alashworth commented 5 years ago

Comment by bob-carpenter Wednesday Aug 17, 2016 at 10:53 GMT

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
alashworth commented 5 years ago

Comment by davharris Wednesday Aug 17, 2016 at 12:44 GMT

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.

alashworth commented 5 years ago

Comment by bob-carpenter Wednesday Aug 17, 2016 at 13:23 GMT

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

alashworth commented 5 years ago

Comment by davharris Wednesday Aug 17, 2016 at 13:56 GMT

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.

alashworth commented 5 years ago

Comment by syclik Wednesday Aug 17, 2016 at 14:10 GMT

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 .

alashworth commented 5 years ago

Comment by davharris Wednesday Aug 17, 2016 at 14:29 GMT

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.

alashworth commented 5 years ago

Comment by bob-carpenter Wednesday Aug 17, 2016 at 14:48 GMT

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.

alashworth commented 5 years ago

Comment by bob-carpenter Friday Aug 19, 2016 at 12:04 GMT

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

alashworth commented 5 years ago

Comment by bob-carpenter Monday Aug 22, 2016 at 15:05 GMT

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.

alashworth commented 5 years ago

Comment by betanalpha Tuesday Aug 23, 2016 at 00:15 GMT

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.

alashworth commented 5 years ago

Comment by betanalpha Tuesday Aug 23, 2016 at 00:44 GMT

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

alashworth commented 5 years ago

Comment by betanalpha Tuesday Aug 23, 2016 at 00:52 GMT

@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.

alashworth commented 5 years ago

Comment by betanalpha Tuesday Aug 23, 2016 at 02:06 GMT

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?

alashworth commented 5 years ago

Comment by bob-carpenter Tuesday Aug 23, 2016 at 12:15 GMT

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(...).

alashworth commented 5 years ago

Comment by bob-carpenter Tuesday Aug 23, 2016 at 12:21 GMT

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.

alashworth commented 5 years ago

Comment by syclik Friday Sep 02, 2016 at 18:52 GMT

@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.

alashworth commented 5 years ago

Comment by syclik Saturday Sep 03, 2016 at 01:34 GMT

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