Numeric overflow

[Fil Krynicki]

Hi everyone.

I've been wrestling with this issue for a couple days now and I'm pretty much out of ideas. I'm using Pystan to fit a model (model string and dat attached). Python version is 3.3.3. 64bit Linux, Kernel 3.12.3-1-ARCH, 8GB memory. gcc is version 4.8.2 and I believe that's what Stan is using (I've seen gcc errors before).

I'm getting a numeric overflow exception when trying to fit a model with 3 chains and 10000 saved steps. The model is fairly simple and the data seems valid, so I'm not really sure what's going on. The model is being fit in a big loop - I'm fitting on the order of 100 models in a long batch job. The exception occurs non-deterministically in this loop, no doubt due to the Monte Carlo nature of the sampler.

Is it reasonable to call a model valid if I continually re-fit until the exception doesn't occur, or should I be changing the model so that it does not numerically overflow 100% of the time. Is there a clear reason why I'm seeing this behaviour?

Thanks for any help.

[Michael Betancourt]

Are you just seeing a warning or is exception causing the code to die? Numerical overflow
should lead to rejected samples, which is fine statistically if inefficient.

For best performance, however, try to make your model more well-posed. Try tighter priors
or reparameterizing so that the bulk of the posterior mass is within a reasonable range.

> --
> 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.
> For more options, visit https://groups.google.com/groups/opt_out.
> <exception_model_string.txt><exception_dat.json>

[Bob Carpenter]

I'm surprised you're getting overflow given the fairly dispersed
data you have. Even some mild priors on alpha and beta should
solve the problem, though --- you're using an improper prior that's
uniform on (0,inf). Gelman et al. discuss priors for beta distributions
in chapter 5 and reparameterize in terms of mean alpha/(alpha + beta)
and total count (alpha + beta) --- we use that as an example in the
reparameterization chapter of the manul.

Unfortunately for this model choice, your data don't look like they come from
a beta distribution --- the data looks multimodal concentrated around .99
and .36 (with all the values around 0.36 apparently identical).
So a simple beta distribution won't be very good for this
data, even when you get the model to fit. So I'm guessing the model runs
off to fit one of the modes really well.

I'd also suggest a constraint on the data variables N and y, because your
model requires N >= 0 and 0 < y[n] < 1. It won't change behavior,
only help with error checking.

- Bob

[Fil Krynicki]

The exception is causing the code to die, which is the problem. The exact exception is:

File "/usr/lib/python3.3/site-packages/pystan/api.py", line 366, in stan

    n_jobs=n_jobs, **kwargs)

  File "/usr/lib/python3.3/site-packages/pystan/model.py", line 679, in sampling

    ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)

  File "/usr/lib/python3.3/site-packages/pystan/model.py", line 81, in _map_parallel

    return list(mymap(function, args))

  File "/usr/lib/python3.3/multiprocessing/pool.py", line 228, in map

    return self._map_async(func, iterable, mapstar, chunksize).get()

  File "/usr/lib/python3.3/multiprocessing/pool.py", line 564, in get

    raise self._value

OverflowError: Error in function boost::math::digamma<d>(d): numeric overflow

Thanks Bob. I'm constrained by a beta model at the moment, but I'll take a run at your suggestions when I return from New Years festivities. Happy Holidays. :)

[Bob Carpenter]

That's a bug in PyStan. I'm not sure if it originates in
Stan or not. I'd run your models and see, but I don't know
how to handle the JSON-formatted data.

What should happen is that the error should be caught and the sample
rejected. You still might not be able to fit the model, but
an overflow in digamma() shouldn't cause the sampler to die.

However we fix Stan, you need a better model. The
computational problem is a symptom of a problem in the model ---
what Andrew calls the "folk theorem". In this case, there are
two issues with the model: (a) no priors on the parameters to
the beta distribution, and (b) the beta distribution doesn't fit
the data. With a prior, you should be able to get a (bad)
fit computationally.

- Bob

On 1/1/14, 12:53 AM, Fil Krynicki wrote:
> The exception is causing the code to die, which is the problem. The exact exception is:
>
> File "/usr/lib/python3.3/site-packages/pystan/api.py", line 366, in stan
> n_jobs=n_jobs, **kwargs)
> File "/usr/lib/python3.3/site-packages/pystan/model.py", line 679, in sampling
> ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
> File "/usr/lib/python3.3/site-packages/pystan/model.py", line 81, in _map_parallel
> return list(mymap(function, args))
> File "/usr/lib/python3.3/multiprocessing/pool.py", line 228, in map
> return self._map_async(func, iterable, mapstar, chunksize).get()
> File "/usr/lib/python3.3/multiprocessing/pool.py", line 564, in get
> raise self._value
> OverflowError: Error in function boost::math::digamma<d>(d): numeric overflow
>
>

> On Tuesday, December 31, 2013 2:40:49 PM UTC-5, Michael Betancourt wrote:
>

> Are you just seeing a warning or is exception causing the code to die? Numerical overflow
> should lead to rejected samples, which is fine statistically if inefficient.
>
> For best performance, however, try to make your model more well-posed. Try tighter priors
> or reparameterizing so that the bulk of the posterior mass is within a reasonable range.
>

> On Dec 31, 2013, at 1:43 PM, Fil Krynicki <filipk...@gmail.com <javascript:>> wrote:
>
> > Hi everyone.
> >
> > I've been wrestling with this issue for a couple days now and I'm pretty much out of ideas. I'm using Pystan to
> fit a model (model string and dat attached). Python version is 3.3.3. 64bit Linux, Kernel 3.12.3-1-ARCH, 8GB memory.
> gcc is version 4.8.2 and I believe that's what Stan is using (I've seen gcc errors before).
> >
> > I'm getting a numeric overflow exception when trying to fit a model with 3 chains and 10000 saved steps. The
> model is fairly simple and the data seems valid, so I'm not really sure what's going on. The model is being fit in a
> big loop - I'm fitting on the order of 100 models in a long batch job. The exception occurs non-deterministically in
> this loop, no doubt due to the Monte Carlo nature of the sampler.
> >
> > Is it reasonable to call a model valid if I continually re-fit until the exception doesn't occur, or should I be
> changing the model so that it does not numerically overflow 100% of the time. Is there a clear reason why I'm seeing
> this behaviour?
> >
> > Thanks for any help.
> >
> > --
> > 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 <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.
> > <exception_model_string.txt><exception_dat.json>

[Bob Carpenter]

On 1/1/14, 12:53 AM, Fil Krynicki wrote:> The exception is causing the code to die, which is the problem. The exact
exception is:
>
> File "/usr/lib/python3.3/site-packages/pystan/api.py", line 366, in stan
> n_jobs=n_jobs, **kwargs)
> File "/usr/lib/python3.3/site-packages/pystan/model.py", line 679, in sampling
> ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
> File "/usr/lib/python3.3/site-packages/pystan/model.py", line 81, in _map_parallel
> return list(mymap(function, args))
> File "/usr/lib/python3.3/multiprocessing/pool.py", line 228, in map
> return self._map_async(func, iterable, mapstar, chunksize).get()
> File "/usr/lib/python3.3/multiprocessing/pool.py", line 564, in get
> raise self._value
> OverflowError: Error in function boost::math::digamma<d>(d): numeric overflow

I just verified that this is a PyStan-specific problem
and opened this issue to track its being fixed:

https://github.com/stan-dev/pystan/issues/53

Running in CmdStan, there are warning messages about the
overflow causing rejection, then sampling settles down.

Here's the fit I get using CmdStan (optimization level 0, which
is why it's so slow):


|
Inference for Stan model: exception_model
2 chains: each with iter=(1000,1000); warmup=(0,0); thin=(1,1); 2000 iterations saved.

Warmup took (14, 13) seconds, 28 seconds total
Sampling took (14, 15) seconds, 29 seconds total

Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
lp__ 784 4.5e-02 1.1 781 784 785 597 21 1.0e+00
accept_stat__ 0.92 2.2e-03 0.10 0.71 0.96 1.0 2000 69 1.0e+00
stepsize__ 0.70 2.8e-02 0.028 0.67 0.73 0.73 1.0 0.034 4.1e+13
treedepth__ 0.81 1.3e-02 0.57 0.00 1.0 2.0 1981 68 1.0e+00
n_divergent__ 0.00 0.0e+00 0.00 0.00 0.00 0.00 2000 69 nan
alpha 0.84 1.2e-03 0.035 0.79 0.84 0.90 788 27 1.0e+00
beta 0.38 4.7e-04 0.013 0.36 0.38 0.40 766 26 1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at
convergence, R_hat=1).
|


The values of alpha and beta are so low because it has to
accomodate the two modes in the data. If the model truly fit the
data, you'd expect (alpha + beta) to be closer to the data size.

Here's the config for the second chain:


|
method = sample (Default)
sample
num_samples = 1000 (Default)
num_warmup = 1000 (Default)
save_warmup = 0 (Default)
thin = 1 (Default)
adapt
engaged = 1 (Default)
gamma = 0.050000000000000003 (Default)
delta = 0.80000000000000004 (Default)
kappa = 0.75 (Default)
t0 = 10 (Default)
init_buffer = 75 (Default)
term_buffer = 50 (Default)
window = 25 (Default)
algorithm = hmc (Default)
hmc
engine = nuts (Default)
nuts
max_depth = 10 (Default)
metric = diag_e (Default)
stepsize = 1 (Default)
stepsize_jitter = 0 (Default)
id = 0 (Default)
data
file = exception.data.R
init = 2 (Default)
random
seed = 4294967295 (Default)
output
file = output2.csv
diagnostic_file = (Default)
refresh = 100 (Default)
|


- Bob


>

[Filip Krynicki]

Thanks for all the help, the response has been well beyond what I could have hoped. I have a very naive question. How would I specify a prior on the alpha and beta parameters? My understanding is that the lower and upper keywords specify hard constraints (or as close as is allowed given numerical precision issues). Do I specify a prior by sampling the alpha and beta parameters from a distribution defined by a hyperparameter?

I'm looking through the manual for any mention of priors and coming up with nothing. My thinking would be (using a normal prior for both parameters for simplicity's sake):

data {

   ...

   real prior_mu;

   real prior_sigma;

}

parameters {

   real alpha;

   real beta;

}

model {

   alpha ~ normal(prior_mu, prior_sigma)

   beta ~ normal(prior_mu, prior_sigma)

   for (i in 1:N) {

      y[i] ~ beta(alpha, beta);

   }

}

This model compiles and runs but I have little confidence that it's doing what I think it is from a theoretical standpoint.

--
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/ia-k26aucy0/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

On 1/1/14, 9:31 PM, Filip Krynicki wrote:
> Thanks for all the help, the response has been well beyond what I could have hoped.

We take bugs seriously, and diagnosing them often involves
us having to dive into the particulars of the models.

We're also interested in bringing more people up to speed in
Bayesian modeling and want to understand where new users have
problems.

> I have a very naive question. How
> would I specify a prior on the alpha and beta parameters? My understanding is that the lower and upper keywords specify
> hard constraints (or as close as is allowed given numerical precision issues).

They restrict alpha and beta to the constrained ranges (both > 0 in this case,
modulo arithmetic underflow, which can make them 0). Then you get a uniform
prior over that range, which is improper in this case because

INTEGRAL_0^infinity const dx = infinity

> Do I specify a prior by sampling the
> alpha and beta parameters from a distribution defined by a hyperparameter?

Yes, but the sampling statements don't actually do sampling. So
you don't actually sample alpha and beta from their priors, per se,
though that's the generative story of your model. In Stan, sampling
statements just add terms to the log probability function.

> I'm looking through the manual for any mention of priors and coming up with nothing.

They're all over the place in the manual. In the
very first example there's a beta prior on theta.

We're assuming a fair degree of background in Bayesian
stats, though there's a chapter near the end of the manual
on Bayesian data analysis that recapitulates the initial
chapter(s) of Gelman et al.'s book of the same name.

> My thinking would be (using a
> normal prior for both parameters for simplicity's sake):
>
> data {
> ...
> real prior_mu;
> real prior_sigma;
> }
>
> parameters {
> real alpha;
> real beta;
> }
>
> model {
> alpha ~ normal(prior_mu, prior_sigma)
> beta ~ normal(prior_mu, prior_sigma)
> for (i in 1:N) {
> y[i] ~ beta(alpha, beta);
> }
> }
>
> This model compiles and runs but I have little confidence that it's doing what I think it is from a theoretical standpoint.

This is OK, theoretically at least. The prior is
technically a normal truncated to the constrained values.

A more typical prior might be an exponential or lognormal, if you
put it direclty on alpha and beta. But it's more typical to do
what we discuss in section 17.1 of the manual, which specifically discusses
priors for the beta distribution in terms of a reparameterization
in terms of the mean (alpha / (alpha + beta)) and
total prior count (alpha + beta). There we use a Pareto prior following
Gelman et al.'s Bayesian Data Analysis, but a lognormal or
exponential (or more geneally gamma) would also make sense.

- Bob

>
>
> On Wed, Jan 1, 2014 at 2:48 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
> On 1/1/14, 12:53 AM, Fil Krynicki wrote:> The exception is causing the code to die, which is the problem. The exact
> exception is:
>
> >

> > File "/usr/lib/python3.3/site-__packages/pystan/api.py", line 366, in stan
> > n_jobs=n_jobs, **kwargs)
> > File "/usr/lib/python3.3/site-__packages/pystan/model.py", line 679, in sampling
> > ret_and_samples = _map_parallel(call_sampler___star, call_sampler_args, n_jobs)
> > File "/usr/lib/python3.3/site-__packages/pystan/model.py", line 81, in _map_parallel
> > return list(mymap(function, args))
> > File "/usr/lib/python3.3/__multiprocessing/pool.py", line 228, in map

> > return self._map_async(func, iterable, mapstar, chunksize).get()

> > File "/usr/lib/python3.3/__multiprocessing/pool.py", line 564, in get

> > raise self._value
> > OverflowError: Error in function boost::math::digamma<d>(d): numeric overflow
>
>
> I just verified that this is a PyStan-specific problem
> and opened this issue to track its being fixed:
>

> https://github.com/stan-dev/__pystan/issues/53 <https://github.com/stan-dev/pystan/issues/53>

> To unsubscribe from this topic, visit https://groups.google.com/d/__topic/stan-users/ia-k26aucy0/__unsubscribe
> <https://groups.google.com/d/topic/stan-users/ia-k26aucy0/unsubscribe>.
> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@__googlegroups.com
> <mailto:stan-users%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.
>
>
> --

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