# User written Box-Cox type transformation function

Hi,

I have written a small function below and have a questions: can Stan do the following computation,“z = ((1 + pos).^p -1)/p”? I could not find any information in the manual on this. I read about this in 2014 and 2016 but it is not clear if this is possible yet.

``````functions {
vector g(vector y, real p) {

int N;
N = rows(y);
vector[N] pos;
vector[N] neg;
vector[N] z;
vector[N] w;
real eps;

pos = y.*(y >= 0);
neg = y.*(y <  0);

eps = 0.005;
if (p >= -eps & p <= eps) z = log(1 + pos);
else z = ((1 + pos).^p -1)/p;

if (p >= 2-eps & p <= 2+eps) w = -log(1 - neg);
else w = -((1 - neg).^(2-p) - 1)/(2-p);

return w+z;
}
}
``````

Operating System: Windows 7
Interface Version: RStan

Thanks,
Raphael

[edit: escape code]

Unless I’m missing something, I don’t think there’s an elementwise power function yet (.^), but you can do this stuff with a for loop and `pow`, sure. Looks like you got if statements and such, so make sure the derivatives with respect to your parameters are continuous across them.

Yes, Stan can calculate this, but you have to stick to Stan syntax. Looks like you’re just experimenting with R syntax, which probably won’t go well. For example, `y > 0` doesn’t exist, nor does `.^` as an operator. You can try to proceed that way, but then you have to take our error messages seriously and correct the things that don’t exist in Stan.

Also, making both the `pos` and `neg` vectors size `N` can be dicey We have `log1p(x)` which is much more stable than `log(1 + x)`.

Loops are fast in Stan.

We have a Box-Cox transform model in the code found at https://osf.io/bcr4u/

Look at the Box-Cox Normal model section.

3 Likes

Thank you all for your comments. Here is the final version of the function that works! This is based on the Yeo-Johnson transformation function which can handle both negative and positive values.

``````functions {

real g(real y, real p) {

real z;
real eps;

eps = 0.005;
if (y >= 0) {
if ( abs(p) <= eps) z = log1p(y);
else z = (pow(1 + y, p) - 1)/p;
}

else {
if (p >= 2-eps && p <= 2+eps) z = -log1m(y);
else z = -(pow(1 - y, 2-p) - 1)/(2-p);
}

return z;
}
}
``````

Yeo, I.K. and Johnson, R.A., 2000. A new family of power transformations to improve normality or symmetry. Biometrika, 87(4), pp.954-959.

Regards,
Raphael

[edit: escaped code]

1 Like

Thanks for sending a working solution. A simpler version with less nesting and fewer assignments is:

``````   real h(real a, real b) {
return (pow(a, b) - 1) / b;
}

real g(real y, real p) {
real eps = 0.005;
if (y >= 0) {
if (abs(p) <= eps)
return log1p(y);
return h(1 + y, p);
}
if (p >= 2 - eps && p <= 2 + eps)
return -log1m(y);
return -h(1 - y, 2 - p);
}
``````

It avoids recomputation of `2 - p` into the bargain.

2 Likes

Thanks Bob!

FYI, `abs` now needs to be `fabs`, else the model won’t compile.