geometry? negative binomial model unstable

[Bob Carpenter]

The negative binomial malloc/free error case is more puzzling.
We hadn't ever used the negative binomial before (at least not
in BUGS examples), so it hasn't been field tested.

The following model winds up with alpha, beta going to +infinity
then segfaulting on a bad free (the seg fault needs to be fixed
as well as the model's bad behavior).

data {
int<lower=1> N;
int<lower=0> y[N];
}
parameters {
real<lower=0.0> alpha;
real<lower=0.0> beta;
}
model {
for (i in 1:N)
y[i] ~ neg_binomial(alpha, beta);
}

And here's the data:

y <- c(0, 1, 4, 0, 2, 2, 5, 0, 1)
N <- 9

The sample mean is 5/3.

This is similar to what I get from Stan if I
hard code the scale term alpha = 1 (that is, beta
around 3/5).

Even so, it's still unstable and has memory errors
with different seeds.

If I constrain alpha and beta to be above 0.1, it doesn't
crash, but adaptation takes forever and the values
of alpha and beta still tend to +infinity.

If I constrain alpha, beta to the interval (0.1,10),
everything runs fine and the estimates for alpha range
from around 1 to around 9.5 with betas around the
sample mean.

If I only constrain on (0,10), I get problems
with memory again (perhaps NaN issues on
some kind of underflow?)

Everything's also stable if I reparameterize to
probability of success and use a strong enough
location=0.5 beta prior.

Is this another geometry issue?

- Bob

[Marcus Brubaker]

I don't have time today to track this down today, but here is what I
found running it under valgrind. Running the model and data normally
produces the crash you're describing. However, when run under
valgrind it produces:

[mbrubake@dynamics stan]$ valgrind --tool=memcheck ./negbin-test
--data=negbin.data==20177== Memcheck, a memory error detector
==20177== Copyright (C) 2002-2010, and GNU GPL'd, by Julian Seward et al.
==20177== Using Valgrind-3.6.0 and LibVEX; rerun with -h for copyright info
==20177== Command: ./negbin-test --data=negbin.data
==20177==
==20177== Conditional jump or move depends on uninitialised value(s)
==20177== at 0x807B361: stan::io::dump_reader::scan_name() (dump.hpp:460)
==20177== by 0x8079D74: stan::io::dump_reader::next() (dump.hpp:769)
==20177== by 0x807974D: stan::io::dump::dump(std::istream&) (dump.hpp:833)
==20177== by 0x8051DB6: int
stan::gm::nuts_command<anon_model_namespace::anon_model>(int, char
const**) (command.hpp:229)
==20177== by 0x805105C: main (negbin-test.cpp:217)
==20177==
==20177== Conditional jump or move depends on uninitialised value(s)
==20177== at 0x807B3D1: stan::io::dump_reader::scan_name() (dump.hpp:460)
==20177== by 0x8079D74: stan::io::dump_reader::next() (dump.hpp:769)
==20177== by 0x807974D: stan::io::dump::dump(std::istream&) (dump.hpp:833)
==20177== by 0x8051DB6: int
stan::gm::nuts_command<anon_model_namespace::anon_model>(int, char
const**) (command.hpp:229)
==20177== by 0x805105C: main (negbin-test.cpp:217)
==20177==
==20177== Use of uninitialised value of size 4
==20177== at 0x2D529F: isalpha (in /lib/libc-2.12.so)
==20177== by 0x807B415: stan::io::dump_reader::scan_name() (dump.hpp:488)
==20177== by 0x8079D74: stan::io::dump_reader::next() (dump.hpp:769)
==20177== by 0x807974D: stan::io::dump::dump(std::istream&) (dump.hpp:833)
==20177== by 0x8051DB6: int
stan::gm::nuts_command<anon_model_namespace::anon_model>(int, char
const**) (command.hpp:229)
==20177== by 0x805105C: main (negbin-test.cpp:217)
==20177==

Exception: variable N not found.
Diagnostic information:
Dynamic exception type: std::runtime_error
std::exception::what: variable N not found.

==20177==
==20177== HEAP SUMMARY:
==20177== in use at exit: 0 bytes in 0 blocks
==20177== total heap usage: 29 allocs, 29 frees, 67,138 bytes allocated
==20177==
==20177== All heap blocks were freed -- no leaks are possible
==20177==
==20177== For counts of detected and suppressed errors, rerun with: -v
==20177== Use --track-origins=yes to see where uninitialised values come from
==20177== ERROR SUMMARY: 3 errors from 3 contexts (suppressed: 15 from 8)


Notice that it crashes for a different reason. It seems like there is
an uninitialized bit of memory being accessed when loading the
parameters for some reason. When it's accessed under valgrind, it is
preinitialized to zero and causes an exception to eventually be
thrown. Under normal conditions it's non-zero and results in a crash
later on for some reason.

Cheers,
Marcus

> --
>
>

[Bob Carpenter]

Thanks for the hints. I'll try to take it from
here. It'll be a good excuse to learn valgrind.

If N is getting read wrong, then it could explain
why there's no data and the model's exploding.

- Bob

[Daniel Lee]

The crash is coming from something else. With seed=1357157083, alpha
gets really small (e-21) and beta gets really big (e+24).

I've narrowed it down to lgamma crashing for a var passed in with a
value of -1.0. This happens when neg_binomial() is called with n = 0,
alpha = e-21, beta = e+24. This calls binomial_coefficient_log() with
-1, 0.

I'll try to figure out what's causing the crash.

> --
>
>

[Marcus Brubaker]

Indeed, my posted trace was a red herring. (I mistyped the data
filename when I called valgrind, d'oh.) There is an issue there when
getting to the end of the file, but it's likely benign.

When running things properly for a seed that causes the crash valgrind
reports an invalid delete():

==22641== Invalid free() / delete / delete[]
==22641== at 0x4005D5B: operator delete(void*) (vg_replace_malloc.c:387)
==22641== by 0x806C32C: stan::agrad::lgamma(stan::agrad::var
const&) (special_functions.hpp:770)
==22641== by 0x8067F60:
boost::math::tools::promote_args<stan::agrad::var, int, float, float,
float, float>::type
stan::math::binomial_coefficient_log<stan::agrad::var,
int>(stan::agrad::var, int) (special_functions.hpp:177)
...
==22761== Address 0x402b120 is 248 bytes inside a block of size 65,536 alloc'd
==22761== at 0x400682F: malloc (vg_replace_malloc.c:236)
==22761== by 0x80B4851: stan::memory::(anonymous
namespace)::eight_byte_aligned_malloc(unsigned int)
(stack_alloc.hpp:47)
==22761== by 0x80B5032:
stan::memory::stack_alloc::stack_alloc(unsigned int)
(stack_alloc.hpp:131)
==22761== by 0x80B4AF3:
stan::memory::stack_alloc::stack_alloc(unsigned int)
(stack_alloc.hpp:139)
==22761== by 0x80B4803: __cxx_global_var_init2 (in
/home/mbrubake/Projects/stan/negbin-test)

I've poked around but can't for the life of me figure out how delete()
is even being called. from there. It may be related to an exception
getting thrown, but that seems unlikely to be trying to delete memory
that was allocated as part of a vari.

Cheers,
Marcus

> --
>
>

[Bob Carpenter]

Daniel's tracked the source of the throw down to
Boost's lgamma function, which doesn't like values
very close to -1.

I'm guessing the vari memory issue is because the
crash happens at a point where we're trying to
allocate a new vari.

The dump() implementation probably deserves some
further investigation. It keeps references to the
data to avoid some copying (but then I think it copies
anyway given our use of Eigen Map).

- Bob

[Jeffrey Arnold]

Thanks for looking into them, let me know if there is any way in which I can help.

Regarding stan dumps I've been using my own function to write the dumps that adds some functionality and resolves

some issues I had when using rstan::stan.rdump.

• catches the following errors
  
• variable does not have a valid Stan name
• variable does not exist
• variable includes NA or NaN
• variable is not numeric
• add a precheck option like base::save
  
• Replace Inf and -Inf with largest and smallest floating point numbers.
  
• Ensure integers are not printed with scientific notation
  
• Written as a S4 generic function to accept either a character vector of object names or a list of objects.

I attached the code below.

[Bob Carpenter]

I see we haven't even scratched the surface
of R's dump format. I have some specific
comments below.

On 8/21/12 3:09 PM, Jeffrey Arnold wrote:
> Thanks for looking into them, let me know if there is any way in which I can help.
>
> Regarding stan dumps I've been using my own function to write the dumps that adds some functionality and resolves
> some issues I had when using rstan::stan.rdump.

You probably just want to use RStan, which won't
have most of these problems, as we read the data
directly from the R environment without dumping
it (through Rcpp).

I'm going to put many of the following on the
C++ to-do list, along with treating single values
as length-one arrays. I really would like this to
flexibly handle R output.

> * catches the following errors
> o variable does not have a valid Stan name

This will just involve periods, or is there more?

> o variable does not exist

That's on the R side.

> o variable includes NA or NaN

NA is R only, and we have no way of supporting it
in Stan. At least not without major work.

I see that R just writes the string "NaN" (without the
quotes). I should update the C++ to handle that.

> o variable is not numeric

Also on the R side.

> * add a precheck option like base::save

What does this mean?

> * Replace Inf and -Inf with largest and smallest floating point numbers.

R just writes "Inf" and "-Inf", so again I should
be able to handle those on the C++ side.

> * Ensure integers are not printed with scientific notation

I should also handle this on the C++ side.

> * Written as a S4 generic function to accept either a character vector of object names or a list of objects.

That's just on the R side.

> I attached the code below.

Thanks.

More on code contributions in a separate e-mail.

- Bob

[Jiqiang Guo]

Since this issue is brought up, I guess at least we should know that in R, objects (even an element of a vector) can have names.  So in the R dump file, there might be names, which I think we could neglect.  An example:

> a <- 1:5

> names(a) <- 1:5

> dump("a", file = 'a.rdump')

> cat(paste(readLines('a.rdump'), collapse = '\n'), '\n')

a <-

structure(1:5, .Names = c("1", "2", "3", "4", "5"))

--

Jiqiang

[Bob Carpenter]

On 8/21/12 9:05 PM, Jiqiang Guo wrote:
> Since this issue is brought up, I guess at least we should know that in R, objects (even an element of a vector) can
> have names. So in the R dump file, there might be names, which I think we could neglect. An example:
>
> > a <- 1:5
> > names(a) <- 1:5
> > dump("a", file = 'a.rdump')
> > cat(paste(readLines('a.rdump'), collapse = '\n'), '\n')
> a <-
> structure(1:5, .Names = c("1", "2", "3", "4", "5"))

RStan isn't reading in dump files, is it?

If so, it should do it using R's reader,
not ours.

- Bob

[Jiqiang Guo]

No.  Generally, it's the user's job to get data into R from whatever source.  If it is a dump file, the user can just call R's function source to read it. 

Jiqiang

[Jeffrey Arnold]

I put my responses below.  For what it is worth, if you're interested, I also have code that parses the header metadata of the csv output, gets the dimensions of the variables from their variable names, and converts table back into lists of arrays.

 

On 8/21/12 3:09 PM, Jeffrey Arnold wrote:
> Thanks for looking into them, let me know if there is any way in which I can help.
>
> Regarding stan dumps I've been using my own function to write the dumps that adds some functionality and resolves
> some issues I had when using rstan::stan.rdump.

You probably just want to use RStan, which won't
have most of these problems, as we read the data
directly from the R environment without dumping
it (through Rcpp).

I've been holding off on using rstan seriously until it stabilizes.  The lack of saving compiled objects had been a deal breaker. But now that that's implemented, I'll probably move away from the command line interface. 

 

I'm going to put many of the following on the
C++ to-do list, along with treating single values
as length-one arrays.  I really would like this to
flexibly handle R output.

>   * catches the following errors
>       o variable does not have a valid Stan name

This will just involve periods, or is there more?

Periods, not ending in __, and not being a C++ / Stan reserved name.   There's also a function that works like make.names in R, but will make valid Stan names. 

 

>       o variable does not exist

That's on the R side.

>       o variable includes NA or NaN

NA is R only, and we have no way of supporting it
in Stan.  At least not without major work.

I see that R just writes the string "NaN" (without the
quotes).  I should update the C++ to handle that.

>       o variable is not numeric

Also on the R side.

>   * add a precheck option like base::save

 

It just checks for the existence of the objects, that they are numeric, and that they have valid names, before writing to the file. I noticed that save() had that option and it seemed like a good idea.