Questions about user-defined functions

[Jonathan Gilligan]

Two questions about user-defined functions:

1. If the function takes a vector argument, is there a way to get the length of the vector so I can loop over its elements? I would like to do something like this:

vector foo(vector x, real r) {
   vector bar;

   for (i in 1:length(x)) {
       bar[i] <- x[i]^(1.-r) / (1.-r);
  }
  return(bar);
}

So as a work-around, I'm trying this to exponentiate a vector:

vector foo(vector x, real r) {
  return exp(log(x) * (1. - r)) / (1. - r);
}

2. I am having problems compiling .stan files with user-defined functions. The attached file compiles from .stan to .cpp with no complaints, but the c++ compiler gives me a long and confusing bunch of error messages, of which this seems most germane:

C:/Users/JonathanGilligan/Documents/R/win-library/3.1/rstan/include/rstan/stan_fit.hpp:1325:9:   instantiated from 'SEXPREC* rstan::stan_fit<Model, RNG_t>::grad_log_prob(SEXP, SEXP) [with Model = model22404d8f2db7_foo_namespace::model22404d8f2db7_foo, RNG_t = boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> >, SEXP = SEXPREC*]'
file22404a2774dd.cpp:648:114:   instantiated from here
C:/Users/JonathanGilligan/Documents/R/win-library/3.1/RcppEigen/include/Eigen/src/Core/Matrix.h:277:7: error: 'YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY' is not a member of 'Eigen::internal::static_assertion<false>

I am attaching an example model file that produces this error.  If I get rid of the functions block and copy and paste the body of the functions into the argument to bernoulli_logit, the model compiles and runs.

Can anyone explain what I'm doing wrong in the function definitions?

I am using rstan 2.4.0, which I'm calling from R 3.1.1

> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggmcmc_0.5.1  ggplot2_1.0.0 reshape2_1.4  plyr_1.8.1    rstan_2.4.0   inline_0.3.13 Rcpp_0.11.2  

loaded via a namespace (and not attached):
 [1] codetools_0.2-8  colorspace_1.2-4 digest_0.6.4     grid_3.1.1       gtable_0.1.2     labeling_0.2     MASS_7.3-33     
 [8] munsell_0.4.2    proto_0.3-10     scales_0.2.4     stringr_0.6.2    tools_3.1.1  

[Bob Carpenter]

On Aug 9, 2014, at 1:24 PM, Jonathan Gilligan <jonathan...@gmail.com> wrote:

> Two questions about user-defined functions:
>
> 1. If the function takes a vector argument, is there a way to get the length of the vector so I can loop over its elements? I would like to do something like this:
>
> vector foo(vector x, real r) {
> vector bar;
>
> for (i in 1:length(x)) {
> bar[i] <- x[i]^(1.-r) / (1.-r);
> }
> return(bar);
> }

You can use size(x) for arrays and vectors.

And rows(x) and cols(x) for matrices.

> So as a work-around, I'm trying this to exponentiate a vector:
> vector foo(vector x, real r) {
> return exp(log(x) * (1. - r)) / (1. - r);
> }
>
> 2. I am having problems compiling .stan files with user-defined functions. The attached file compiles from .stan to .cpp with no complaints, but the c++ compiler gives me a long and confusing bunch of error messages, of which this seems most germane:
>
> C:/Users/JonathanGilligan/Documents/R/win-library/3.1/rstan/include/rstan/stan_fit.hpp:1325:9: instantiated from 'SEXPREC* rstan::stan_fit<Model, RNG_t>::grad_log_prob(SEXP, SEXP) [with Model = model22404d8f2db7_foo_namespace::model22404d8f2db7_foo, RNG_t = boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> >, SEXP = SEXPREC*]'
> file22404a2774dd.cpp:648:114: instantiated from here
> C:/Users/JonathanGilligan/Documents/R/win-library/3.1/RcppEigen/include/Eigen/src/Core/Matrix.h:277:7: error: 'YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY' is not a member of 'Eigen::internal::static_assertion<false>


>
> I am attaching an example model file that produces this error. If I get rid of the functions block and copy and paste the body of the functions into the argument to bernoulli_logit, the model compiles and runs.
>
> Can anyone explain what I'm doing wrong in the function definitions?
>
> I am using rstan 2.4.0, which I'm calling from R 3.1.1

Thanks for the detailed info.

This looks like a bug. Let me get back to you with a full report.

- Bob

[Bob Carpenter]

On Aug 9, 2014, at 1:24 PM, Jonathan Gilligan <jonathan...@gmail.com> wrote:

> Two questions about user-defined functions:
>
> 1. If the function takes a vector argument, is there a way to get the length of the vector so I can loop over its elements? I would like to do something like this:
>
> vector foo(vector x, real r) {
> vector bar;
>
> for (i in 1:length(x)) {
> bar[i] <- x[i]^(1.-r) / (1.-r);
> }
> return(bar);
> }
>
> So as a work-around, I'm trying this to exponentiate a vector:
> vector foo(vector x, real r) {
> return exp(log(x) * (1. - r)) / (1. - r);
> }
>
> 2. I am having problems compiling .stan files with user-defined functions. The attached file compiles from .stan to .cpp with no complaints, but the c++ compiler gives me a long and confusing bunch of error messages, of which this seems most germane:

I isolated the bug and created a small reproducible example,
which I included in the issue along with an explanation:

https://github.com/stan-dev/stan/issues/849

For now, you can workaround the bug with this modified function:

vector eu_item(vector x, real r) {
if (r == 1.) {
vector[cols(x)] log_x;
log_x <- log(x);
return log_x;
// WAS JUST: return log(x);
} else {
return exp(log(x) * (1.-r)) / (1.-r); # x^(1.-r) / (1.-r)
}

The problem arises in the "return log(x)" when x is data and
r is a parameter. The return scalar type for log(x) is double (data),
but it needs to be promoted to an autodiff variable (parameter).

This should be fixed by the next release.

Having said all this, I think there's going to be an issue with
your function, because if r is real, it should never be exactly 1,
as that's a measure 0 subset of the probability mass for it. So you
could probably just drop this altogether.

Are there constraints on what r's value should be? If r > 1, that
looks like it could be problematic.

- Bob

[Jonathan Gilligan]

Thanks.

I was going to ask a followup about the "r=1" case, but wanted to test the code myself first.

This is a pretty standard utility function from economics, describing constant relative risk-aversion and plausible values span 1. The function is weird because what I really want is to know if f(x,r) > f(y,r) for some parameter r and data x: bernoulli_logit( scale * (f(x,r) - f(y,r))), where f(x,r) = x^(1-r) / (1-r) . But the function blows up near r=1 because it divides by zero. The special case for "r=1" is to avoid this singularity, and is the standard treatment in the econ literature. I probably should deal with the measure-zero issue by adjusting to use a neighborhood:

if (fabs(r=1.0) < epsilon)
  return log(x);
else
  return pow(x, 1.0 - r) / (1.0 - r);

r > 1 (x to a negative power) is not an issue because x >= 0  and I can avoid the x=0 case by adding a small number to make x strictly greater than 0.

If I could avoid r = 1, that would be nice, but the literature reports values that range from 0 to 40 or higher for different conditions and based on a lot of that literature, r=1 is a pretty reasonable mode for a prior.  A published MLE of the data set I'm working with gives a 95% CI of 0.67 to 1.02 for r.

How badly will this function screw up Stan's ability to calculate Jacobians? Do you have any suggestions for how to code something like this efficiently in Stan?

[Jonathan Gilligan]

Oops. Typo.

The code snippet should read

if (fabs(r - 1.0) < epsilon)


  return log(x);
else
  return pow(x, 1.0 - r) / (1.0 - r);


[Bob Carpenter]

[expanded subject --- Michael see below]

On Aug 10, 2014, at 5:40 PM, Jonathan Gilligan <jonathan...@gmail.com> wrote:

> Thanks.
>
> I was going to ask a followup about the "r=1" case, but wanted to test the code myself first.
>
> This is a pretty standard utility function from economics, describing constant relative risk-aversion and plausible values span 1. The function is weird because what I really want is to know if f(x,r) > f(y,r) for some parameter r and data x: bernoulli_logit( scale * (f(x,r) - f(y,r))), where f(x,r) = x^(1-r) / (1-r) . But the function blows up near r=1 because it divides by zero. The special case for "r=1" is to avoid this singularity, and is the standard treatment in the econ literature. I probably should deal with the measure-zero issue by adjusting to use a neighborhood:
>
> if (fabs(r=1.0) < epsilon)
> return log(x);
> else
> return pow(x, 1.0 - r) / (1.0 - r);
>

this has to be coded with "r == 1.0", but the compiler
will tell you that.

> r > 1 (x to a negative power) is not an issue because x >= 0 and I can avoid the x=0 case by adding a small number to make x strictly greater than 0.

If you have a constraint lower=0 on x, then x will only be 0 if you initialize to
0 (which is an error in most cases) or if there's an underflow.

> If I could avoid r = 1, that would be nice, but the literature reports values that range from 0 to 40 or higher for different conditions and based on a lot of that literature, r=1 is a pretty reasonable mode for a prior. A published MLE of the data set I'm working with gives a 95% CI of 0.67 to 1.02 for r.
>
> How badly will this function screw up Stan's ability to calculate Jacobians? Do you have any suggestions for how to code something like this efficiently in Stan?

Maybe Michael will know.

I'm not sure where you mean. The only Jacobian (determinants) we compute
are for the transforms, and that doesn't depend on anything past the
parameters block.

It will screw up gradients at r = 1 as the function's discontinuous.
I'm not sure what happens under the crossing, because of the sign change.
With a fixed x and r = 1 + epsilon, you get a really big value, whereas
with a fixed x and r = 1 - epsilon, you get a really small value. So the
function's discontinuous at r, which is a problem. It's impossible to
simulate a particle moving according to the Hamiltonian when it has to
cross a discontinuous boundary. In fact, we want the derivatives to
be continuous, as well as the function (which is why fabs() can be problematic).

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

[Tamas Papp]

Hi,

Maybe I missed something in the earlier discussion, but the
CRRA/isoelastic utility function family defined as

log(c) if r=1
(c^(1-r) - 1)/(1-r) otherwise

_is_ continuous both in c and r. The standard way of showing this is
using L'Hopital's rule. The -1 in the numerator is important, and seems
to be missing in the definition below.

The function is also continuously differentiable in c (u'(c)=c^(-r)),
and also in r --- this can be shown by calculating

\lim_{r\to1} \frac{(c^(1-r) - 1)/(1-r)-log(c)}{r-1} = -\log^2(c)/2

Best,

Tamas

[Bob Carpenter]

The function in question in the Stan program was this:

vector eu_item(vector x, real r) {
if (r == 1.) {

return log(x);
} else {
return exp(log(x) * (1.-r)) / (1.-r); # x^(1.-r) / (1.-r)
}

which is missing the -1 term in the numerator compared to Tamas's
version (below).

- Bob

[Tamas Papp]

Within the context of an (economic, not statistical) model, the CRRA
coefficient (r below) is usually fixed, and affine transformations don't
affect comparisons of utility, so it is common to drop the -1. But it is
needed for a continuous family. On top of that, economists rarely
consider differentiability in the coefficient, I had to do the proof for
continuity myself before posting, just to be sure.

I wonder if this solves the OP's problem.

Best,

Tamas

[Jonathan Gilligan]

Thanks, Tamas.

I screwed up in my original question (and my model) by leaving out the
-1 in the numerator. Thanks for catching my error.

With this, I should be able to write a function that will be fine
because the only problem is the measure-zero set (r=1) and the
probability of drawing a sample sufficiently close to r=1 to blow up the
calculation should be small enough that I can ignore it.

I really appreciate your catching my error.

Jonathan

[Bob Carpenter]

The main place to watch out for 1 values is in initialization,
because it's a round number and looks near the value you're
expecting.

- Bob

[Jonathan Gilligan]

Thanks for the pointer.

I wrote an initialization function in R that draws initial values from
my (weakly informative) priors, and then checks them for validity before
passing them to Stan. That seems to do the trick for the most part.

I still do get a few error messages from Stan in the early iterations
because something is over or underflowing and passing #IND values to
bernoulli_logit_log but once Stan has evaluated 100 or so samples, the
chains run with no further error messages, and looking at the traceplots
and density plots, it looks as though the chains mix nicely and don't
get stuck above or below 1, and with a few thousand samples, the
posterior density is smooth where it crosses 1 and Rhat for all the
parameters and lp__ are very close to 1.

Thanks for all the help and advice. Once I fixed the error Tamas pointed
out and got rid of the if statement, the stan computations sped up by a
couple of orders of magnitude (in effective samples per second).

Jonathan

--
Jonathan Gilligan
251 33rd Av N cell 615.414.5750
Nashville TN 37209 land 615.297.2960

To a cyclist, a prius is just a small hummer --- Matt Ruscigno