User written Box-Cox type transformation function

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

3 Likes
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]

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

2 Likes
7

Thanks Bob!

8

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