# User written Box-Cox type transformation function

#1

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]

#2

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.

#3

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.

#4

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

Look at the Box-Cox Normal model section.

Problem recovering generated predictions for box-cox model
#5

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]

#6

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.

Reversing Yao-Johnson normality transform?
#7

Thanks Bob!

#8

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