Evaluating +inf and -inf to prevent overflow and underflow in Generated Quantities?

[Max Farrell]

Hello!

I'm working on a hierarchical binomial-logit model and using the binomial_rng function in the generated quantities block for my posterior predictive checks. I'm also going to try out loo, so I included the log_lik calculation.

I've been getting a few instances of +infinity or -infinity when running my model, so to try to control for these in the generated quantities block I implemented the following:

generated quantities{

         vector[N] log_lik;

    int y_pred[N];

    real p_pred[N];

    for (n in 1:N) {

        p_pred[n] = inv_logit(beta0 + beta1 * predictor[n]);

       
        if (p_pred[n] == positive_infinity()) p_pred[n] = 1;

        if (p_pred[n] == negative_infinity()) p_pred[n] = 0;

       
        y_pred[n] = binomial_rng(cases[n], p_pred[n]);

        log_lik[n] = binomial_lpmf(deaths[n] | cases[n], p_pred[n]);

    }

}

This results in very nice posterior predictions, however, when I remove the if statements, then my predictions get much worse, and it looks like I'm underpredicting instances where the the response is zero and when the response is 100% of the number of cases.

I had the privilege to meet with Bob, who said I shouldn't need to use these control statements, but it seems like without them, instances with extremely high or extremely low p_pred are omitted from my generated quantities. 

Any thoughts? Is it okay to use these kind of controls, or should I go about this differently?

Max

[Ben Goodrich]

On Wednesday, April 12, 2017 at 1:38:09 PM UTC-4, Max Farrell wrote:

Any thoughts? Is it okay to use these kind of controls, or should I go about this differently?

I am not sure why inv_logit is overflowing, but you could call bernoulli_logit_rng() cases[n] times, which inputs the log-odds, and then sum to get a binomial realization. There is no binomial_logit_rng() although there is an issue to create one.

[Bob Carpenter]

inv_logit itself isn't very stable, as it's just this:

inv_logit(x) = 1 / (1 + exp(-x))

No way to stop that exp from overflowing or rounding under
to zero. Around x = -700 and we round to zero, and we should
round to 1 even sooner on the postive x side, probably around 25 or
30 or so.

What we may need in the RNG is something that tests for whether
the argument is 0 or 1 and just returns a string of those results.

- bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Wednesday, April 12, 2017 at 3:02:33 PM UTC-4, Bob Carpenter wrote:

inv_logit itself isn't very stable, as it's just this:

  inv_logit(x) = 1 / (1 + exp(-x))

No way to stop that exp from overflowing or rounding under
to zero.  Around x = -700 and we round to zero, and we should
round to 1 even sooner on the postive x side, probably around 25 or
30 or so.

I meant how is the OP getting values of inv_logit that overflow to positive_infinity() or negative_infinity() ?

[Bob Carpenter]

The code in the stan/math/prim/scal/fun/inv_logit.hpp is just:

inline double inv_logit(double a) {
using std::exp;
return 1.0 / (1.0 + exp(-a));
}

That should not produce an infinite result. If exp(-a) rounds
to zero, it returns 1, whereas if it goes to infinity, it should return
zero.

I'd like to see the value of the argument (beta0 + beta1 * predictor[n])
and the result p_pred[n] here if those tests succeed. So could you try
changing the generated quantities test to print like this:

p_pred[n] = inv_logit(beta0 + beta1 * predictor[n]);

if (p_pred[n] == positive_infinity()) {
print("+inf: ", "beta0 + beta1 * predictor[", n, "]=", beta0 + beta1 * predictor[n],
" p_pred = ", p_pred[n]);
p_pred[n] = 1;
}

if (p_pred[n] == negative_infinity()) {
print("-inf: ", "beta0 + beta1 * predictor[", n, "]=", beta0 + beta1 * predictor[n],
" p_pred = ", p_pred[n]);
p_pred[n] = 0;
}

and see if those tests ever get called.

- Bob

[Max Farrell]

Hi Bob and Ben, 

I ran Bob's updated GQ code and nothing got called up. I'm not sure what's going on, but it could also be due to some of the other modifications Bob suggested (i.e. vectorizing intermediate levels of the hierarchical model). 

I'm not sure if this is still the case, but I found an older post saying that as of yet intermediate hierarchical groupings cannot be vectorized (http://stackoverflow.com/questions/29379001/nested-model-in-stan). If this is the case, I might have jumped the gun on posting this error... My apologies!

However, I still do get a few errors thrown about instances of -inf for my p value in binomial_logit(cases, p), though these do not seem to cause problems in the generated quantities block.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/0rZKv0QieR4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

> On Apr 13, 2017, at 5:15 PM, Max Farrell <maxwell...@gmail.com> wrote:
>
> Hi Bob and Ben,
>
> I ran Bob's updated GQ code and nothing got called up. I'm not sure what's going on, but it could also be due to some of the other modifications Bob suggested (i.e. vectorizing intermediate levels of the hierarchical model).

Vectorization won't change the numerical results. But it can
make it easier to code things properly.

> I'm not sure if this is still the case, but I found an older post saying that as of yet intermediate hierarchical groupings cannot be vectorized (http://stackoverflow.com/questions/29379001/nested-model-in-stan). If this is the case, I might have jumped the gun on posting this error... My apologies!

You can vectorize all this now. I updated the stackoverflow
answer. The manual's where to look for the latest on what is
and isn't possible in Stan.

> However, I still do get a few errors thrown about instances of -inf for my p value in binomial_logit(cases, p), though these do not seem to cause problems in the generated quantities block.

How is p being constructed?

- Bob

[Max Farrell]

Hi Bob, 

I tried my hand at vectorizing the model once more, but again have run up against divergent transitions, and some Rhats of NaN for some of my generated quantities. 

I couldn't find a specific entry in the most recent manual about vectorizing nested hierarchical levels, so I did what I think would be the most obvious syntax.

So that you can better see how p is being constructed, attached are my original and vectorized models. The original nonvectorized model I ended up moving the intermediate grouping level "parasite" from transformed parameters into the model block, then had to use a normal_rng to get my generated quantities for the parasite effects - (This feels wrong to me as I want  to use the estimated parasite effects, not generate random ones for my ppcs, but it was the only way I could get the model to run without divergent transitions). 

I have also attached my ideal, vectorized, model that is fully vectorized and using the non-centered parameterization throughout. With this model, I get about 600 divergent transitions when running for 10k iterations, and some Rhats of NaN for my generated quantities. 

I've reverted to using my original model for now, but still wonder if / why my attempt at vectorization is causing these errors.

Max  

[Bob Carpenter]

> On Apr 18, 2017, at 3:35 PM, Max Farrell <maxwell...@gmail.com> wrote:
>
> Hi Bob,
>
> I tried my hand at vectorizing the model once more, but again have run up against divergent transitions, and some Rhats of NaN for some of my generated quantities.
>
> I couldn't find a specific entry in the most recent manual about vectorizing nested hierarchical levels, so I did what I think would be the most obvious syntax.

I don't think there's anything special about which level you're
at. The vectorization always works the same way.

> So that you can better see how p is being constructed, attached are my original and vectorized models. The original nonvectorized model I ended up moving the intermediate grouping level "parasite" from transformed parameters into the model block, then had to use a normal_rng to get my generated quantities for the parasite effects - (This feels wrong to me as I want to use the estimated parasite effects, not generate random ones for my ppcs, but it was the only way I could get the model to run without divergent transitions).

You absolutely need to use normal_rng for the posterior predictive
checks---it's the only way to draw from the proper posterior predictive
distribution.

>
> I have also attached my ideal, vectorized, model that is fully vectorized and using the non-centered parameterization throughout. With this model, I get about 600 divergent transitions when running for 10k iterations, and some Rhats of NaN for my generated quantities.
>
> I've reverted to using my original model for now, but still wonder if / why my attempt at vectorization is causing these errors.

Vectorization in and of itself won't cause divergences. It will generate
the same density. So presumably the two programs are not calculating the
same densities. So I'd try to isolate where they differ by printing out
intermediate results or breaking things out into functions to evaluate on
their own.

- Bob