gamma_cdf and gamma_cdf_log

[Christopher Steele]

Hi all,

I am running a model that requires the cdf of a gamma distribution:

gamma_cdf_log(threshold,dose[j,h]/scale,1/scale);

When I run with gamma_cdf I get an overflow error, which is not suprising as I had the same problem when I previously implemented the model in C++. However, I also get the same error when I run with gamma_log_cdf:

Error in function boost::math::tgamma<e>(e): Result of tgamma is too large to represent.

I was under the impression that gamma_log_cdf should use lgamma instead of tgamma, or at least that is what my C++ code does. 

Any help would be appreciated.

Chris Steele

> sessionInfo()

R version 3.1.3 (2015-03-09)

Platform: x86_64-unknown-linux-gnu (64-bit)

Running under: CentOS release 6.4 (Final)

locale:

 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              

 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    

 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   

 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 

 [9] LC_ADDRESS=C               LC_TELEPHONE=C            

[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] likeLTD_6.0.0   rtf_0.4-11      gtools_3.1.0    ggplot2_0.9.3.1

[5] DEoptim_2.2-2   rstan_2.6.0     inline_0.3.13   Rcpp_0.11.5    

loaded via a namespace (and not attached):

 [1] codetools_0.2-10   colorspace_1.2-4   dichromat_2.0-0    digest_0.6.3      

 [5] gdata_2.13.2       grid_3.1.3         gtable_0.1.2       labeling_0.2      

 [9] MASS_7.3-39        munsell_0.4.2      plyr_1.8           proto_0.3-10      

[13] RColorBrewer_1.0-5 reshape2_1.2.2     R.methodsS3_1.6.1  R.oo_1.18.0       

[17] scales_0.2.3       stats4_3.1.3       stringr_0.6.2      tcltk_3.1.3       

[21] tools_3.1.3       

[Christopher Steele]

Note, it is definitely this line causing the problem, because if I comment out the line the model runs just fine (using gamma densities elsewhere in the model).

[Bob Carpenter]

Our log CDFs and CCDFs are not that robust yet.
They are just implemented as logs of the CDF and
log(1 - CDF). There's already an issue to fix these,
but it's not on our short-term to-do list.

One thing you can do is try to constrain whichever parameter
is causing a problem.

- 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.
> For more options, visit https://groups.google.com/d/optout.

[Christopher Steele]

Thanks for the information, I think I will try to implement a true log cdf if further constraining does not work. 

Chris

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

[Bob Carpenter]

You might want to check the gamma_cdf_log.hpp in

stan/src/stan/math/prim/scal/prob/

This one actually looks reasonable in that it's not just a naive
log of the gamma_cdf.

The problem you might be having is that you often get overflow or
underflow of derivatives when you have log CDFs, because in the
tails you wind up dividing by the (non-log) CDF.

- Bob

[Christopher Steele]

Thanks for the tip, how would I go about calling the function from the .hpp file?

Chris

[Bob Carpenter]

You need to compile using the same includes as our makefiles,
then you should just be able to pass in double values if you just
want to test behavior.

It's a little trickier if you want to also evaluate gradients.
Then you either need to use our functionals and create an appropriately
functor or use the library directly.

You can also see derivative behavior by using finite differences on
an implementation you trust.

- Bob

[Christopher Steele]

Sorry for being thick, but how exactly would I do that? I have been looking around the internet and through the stan manual but cannot find anything.

[Bob Carpenter]

Everything's header only in Stan's C++, so it's just a matter
of getting the right includes. If STAN_HOME is where Stan
is installed, then the libraries you need are here

EIGEN_HOME = $(STAN_HOME)/lib/eigen_3.2.4
BOOST_HOME = $(STAN_HOME)/lib/boost_1.55.0

and the paths to include are

-I. -I$(STAN_HOME)/src -I$(BOOST_HOME) -I$(EIGEN_HOME)

We don't have any user-facing doc for the C++ API. It sort
of starts from the makefiles and the unit-test files. Running
through the unit tests would be anotehr alternative. They run
from a top-level python script which includes its own doc.

- Bob

[Christopher Steele]

I had a go at this but it doesn't seem to have worked. I downloaded the stan source files from your github, and linked to the eigen and boost libraries when compiling my model using the Rstan function stan(), but I still get the tgamma errors when I run my model. 

So instead I have been looking into coding up a log gamma cdf function in C++. I read on another mailing list post that any functions you use in C code has to be understood by autodiff. Does this mean I cannot use binary arithmetic, break statements and the like? Or are there functions that autodiff understands, but are not included in the stan manual?

Thanks for ll the help!

[Daniel Lee]

Hi Christopher,

On Tuesday, April 21, 2015, Christopher Steele <steele...@googlemail.com> wrote:

I had a go at this but it doesn't seem to have worked. I downloaded the stan source files from your github, and linked to the eigen and boost libraries when compiling my model using the Rstan function stan(), but I still get the tgamma errors when I run my model. 

Unless you've rebuilt RStan from scratch, I would expect the interface to change just based on downloading the Stan source code from GitHub. The instructions that Bob posted were for compiling C++ directly from the command line. If you're trying to write custom C++, you'll have to start with the command line before trying to get it integrated through RStan. It can be done, it just requires more config and more steps. 

 

So instead I have been looking into coding up a log gamma cdf function in C++. I read on another mailing list post that any functions you use in C code has to be understood by autodiff. Does this mean I cannot use binary arithmetic, break statements and the like? Or are there functions that autodiff understands, but are not included in the stan manual?

The Stan C++ api has binary arithmetic, you can use break statements, and we've also reimplemented all of the standard math library, so that's all available, the documentation is a little light, but what we have is written in doxygen and if you don't have doxygen installed, grab a tagged version of the repository. It'll have the generated html webpages under doc/api.

There's also a wki page that describes the use of the autodiff library. 

[Bob Carpenter]

> On Apr 21, 2015, at 9:22 PM, Christopher Steele <steele...@googlemail.com> wrote:
>
> I had a go at this but it doesn't seem to have worked. I downloaded the stan source files from your github, and linked to the eigen and boost libraries when compiling my model using the Rstan function stan(), but I still get the tgamma errors when I run my model.

I'm not sure what you were trying. That function's just what gets called
in a Stan program.

>
> So instead I have been looking into coding up a log gamma cdf function in C++. I read on another mailing list post that any functions you use in C code has to be understood by autodiff. Does this mean I cannot use binary arithmetic, break statements and the like?

Break statements are fine. By binary arithmetic, you mean bit fiddling like the
binary | and & or << and <<< shift operations?

The trick is that any functions or operations that are involved in the result
have to be auto-diffable. So if you bit fiddle to do something like choose which
iterative algorithm to run, that'll probably be OK, but you'll have to be careful
in how you pull the bits out.

We're happy to advise if you have something specific in mind, like a stable
C++ algorithm that works on double values that you just want to translate to
something that can be autodiffed.

And it's no holds barred if you implement custom derivatives for the code. I think
you'll find the problem is that there's not a very big range through which CDFs can
be autodiffed due to floating point representations themselves (that is, they don't
have enough range to represent the derivatives in many cases).

> Or are there functions that autodiff understands, but are not included in the stan manual?

Nope --- they're all in the manual. But control's separate --- you can use
the regular loops, while, conditionals, etc., the ternary operator, etc. You just
can't import a third-party library and use its functions.

- Bob

[Christopher Steele]

Thanks for the help, I am going to try reparameterizing the model to use the exponential cdf for now, but I may come back to this in the future.