Speeding up the Cauchy function

[Ilya Esteban]

I am working on a regression model where I am using Cauchy priors for the coefficients. Since the built-in Cauchy function is quite slow, I am trying to refactor my code to express the Cauchy distribution as a ratio of two Normal distributions. However, the parser isnt allowing me to divide two vectors. Is the only way to do this to express the division as a an explicit loop over all the elements of the vector? Are there any better ways of speeding up sampling from the Cauchy? Thanks in advance!

Here's what the original "slow" code looks like:

data {
 int<lower=1> N; // number of observations
 int<lower=1> K; // number of features
 matrix[N,K] x; // input data
 int<lower=0,upper=1> y[N]; // outcome class
}
parameters {
 real intercept;
 vector[K] coeffs;
}
model {
 intercept ~ normal(0, 1);
 coeffs ~ cauchy(0, .01);
 y ~ bernoulli_logit(intercept + x * coeffs);
}

and here's what my non-working attempt to speed up the code looks like:

data {
 int<lower=1> N; // number of observations
 int<lower=1> K; // number of features
 matrix[N,K] x; // input data
 int<lower=0,upper=1> y[N]; // outcome class
}
parameters {
 real intercept;
 vector[K] coeffs_numerator;
 vector[K] coeffs_denominator;
}
model {
 intercept ~ normal(0, 1);
 coeffs_numerator ~ normal(0,.01);
 coeffs_denominator ~ normal(0,1);
 y ~ bernoulli_logit(intercept + x * (coeffs_numerator / coeffs_denominator));
}

[Ben Goodrich]

On Wednesday, December 26, 2012 3:51:16 PM UTC-5, Ilya Esteban wrote:

I am working on a regression model where I am using Cauchy priors for the coefficients. Since the built-in Cauchy function is quite slow, I am trying to refactor my code to express the Cauchy distribution as a ratio of two Normal distributions. However, the parser isnt allowing me to divide two vectors. Is the only way to do this to express the division as a an explicit loop over all the elements of the vector?

You need to use the "./" operator, which is documented on page 149 of the v1.1.0 manual.
 

Are there any better ways of speeding up sampling from the Cauchy? Thanks in advance!

I'm not sure that the log-posterior will evaluate much faster, but you may get more efficient sampling and hence need fewer total samples to achieve a given number of effective samples.

Ben

[Bob Carpenter]

> On Wednesday, December 26, 2012 3:51:16 PM UTC-5, Ilya Esteban wrote:
>

> ... I am trying to refactor my code to express the Cauchy distribution as a ratio of two Normal
> distributions....

>
> Are there any better ways of speeding up sampling from the Cauchy? Thanks in advance!

The other problem with the ratio is that it's not identified.
Instead, I'd recommend relying directly on the inverse CDF as
Matt's been recommending, which is quite simple for the Cauchy,
assuming the Wikipedia page is right:

http://en.wikipedia.org/wiki/Cauchy_distribution#Related_distributions

If

u ~ uniform(-pi/2, pi/2), and

v = a + b * tan(y),

then

v ~ cauchy(a,b);

To express in Stan, you'd declare:

parameters {
real<lower = -pi()/2, upper = pi()/2> u;
}
transformed parameters {
real v;
v <- a + b * tan(u);
}

where a and b are declared as type real in
the data, transformed data, parameters, or transformed
parameters block. Of course, you don't need a if it's
a constant 0 or b if it's a constant 1.

We haven't vectorized tan(y) yet, but when we do it
probably won't be any faster than just transforming
in a loop. Also, the transformed parameter's just there
for convenience to print it. You can either make it
a local variable (slightly faster as there's no printing
or range checking) or you can just use the expression
(a + b * tan(u)) where you would've used v and skip the
intermediate altogether.

- Bob

[Ilya Esteban]

That made a dramatic difference in speed. Thank you.

One quick follow-up -- if I want to extend my model and express the scale parameter for the coefficients as a half-Cauchy distribution, is this the right way to do it?: (please excuse my notation, but I am a programmer, and not a statistician)

data {
 int<lower=1> N; // number of observations
 int<lower=1> K; // number of features

 real<lower=0> S; // scale parameter

 matrix[N,K] x; // input data
 int<lower=0,upper=1> y[N]; // outcome class
}
parameters {
 real intercept;

 real<lower = -pi()/2, upper = pi()/2> u[K];  
 real<lower = -pi()/2, upper = pi()/2> scale;
}
transformed parameters {
 real<lower=0> cScale;

 cScale <- fabs(S * tan(scale));

 vector[K] coeffs;
 for (i in 1:K)
 {
  coeffs[i] <- cScale * tan(u[i]);
 }
}
model {
 intercept ~ normal(0, 2);

 y ~ bernoulli_logit(intercept + x * coeffs);
}

[Ben Goodrich]

On Thursday, December 27, 2012 9:32:09 PM UTC-5, Ilya Esteban wrote:

That made a dramatic difference in speed. Thank you.

One quick follow-up -- if I want to extend my model and express the scale parameter for the coefficients as a half-Cauchy distribution, is this the right way to do it?: (please excuse my notation, but I am a programmer, and not a statistician) 

parameters {
 real intercept;
 real<lower = -pi()/2, upper = pi()/2> u[K];  
 real<lower = -pi()/2, upper = pi()/2> scale;
}
transformed parameters {
 real<lower=0> cScale;

 cScale <- fabs(S * tan(scale));

 

I think what you wrote is not wrong, but it would be better to just restrict the support of "scale" to the (0,pi/2) interval than to take the absolute value of S * tan(scale).  Your way would imply that the posterior distribution of scale is bimodal with a reflection around zero.

Ben

[Ilya Esteban]

Just out of curiosity -- given that it takes two lines of code to transform the uniform distribution into the Cauchy, is there any particular reason why the built-in Cauchy function is implemented differently and is so much slower?

Thanks,

Ilya

On Wednesday, December 26, 2012 5:55:10 PM UTC-8, Bob Carpenter wrote:

[Bob Carpenter]

On 1/22/13 4:26 AM, Ilya Esteban wrote:
> Just out of curiosity -- given that it takes two lines of code to transform the uniform distribution into the Cauchy, is
> there any particular reason why the built-in Cauchy function is implemented differently and is so much slower?

The problem isn't that Stan's Cauchy density is inefficiently
implemented. It's the geometry, as determined (to second order)
by the Hessian of the log probability of the unconstrained
parameters. (Any parameters declared with constraints are
automatically transformed to unconstrained parameters so HMC
can operate with unconstrained support; the details of the
specific transforms are in the manual toward the back.)

By declaring a uniform variable with constraint that it be
between 0 and 1, then transforming it to a Cauchy implies
a different geometry than defining a Cauchy over an unconstrained
parameter. It's not clear this particular transform will even
help. Of course, there are then other ways to do the [0,1]
transform including inverse CDFs.

The reason he built-in Cauchy function can't do transforms
automatically is that it has no idea what the user's intent
is for the log probability function!

As we learn more about how the geometry affects sampling and which
numerical issues come up, we will hopefully be able to provide
pre-packaged solutions that aren't as clunky as defining
transforms.

Another thing we could do is perform variable transformations
under the hood. That's effectively what we're doing for
parameters declared with constraints. But I'd prefer to keep
things as direct as possible.

- Bob

> --
>
>