Using the skew-t to model error

[Sam Dickson]

The paraphrased version of my question is this: I'm trying to use the skew-t distribution to fit the error in a model and I think I might be doing it wrong. I've fit the error with a t-distribution and a skew-normal, but when I try to do the custom distribution I get an error.  I think the problem must lie in my implementation of increment_log_prob() function, which is shown in the following code:

 

  for (n in 2:N){
    x[n] ~ student_t(nu,mu[n],sigma);
    increment_log_prob(-log(sigma)
                       + student_t_cdf_log(skew*z[n]*sqrt((nu+1)/(pow(z[n],2)+nu)),nu+1,0,1));
  }

What is the correct way to use the skew-t distribution to model the error of the model?

 

 

The long format of my question is below.

 

I'm trying to fit a multiple components cosinor model to telemetry data (like body temperature or heart rate).  I have 5 days of baseline data measured every half hour.  I first use a Lomb-Scargle periodogram to look for periodicity at 24 hours, 12 hours, and 8 hours, then attempt to fit the model accordingly.  I have an AR(1) parameter to account for autocorrelation.  I started of by fitting the error with a t distribution in order to allow for occasional outliers in the data, so if there was 24- and 12-hour periodicity, my stan code would look like this:

 

data {
  int<lower=0> N;
  real x[N];
  real time[N];
  real periodc;
  real npi;
  real ppi;
  real<lower=0> sigmaLow;
  real<lower=0> sigmaHigh;
  real smean;
  real<lower=0> nuParm;
  real<lower=0> rgx;
}
parameters {
  real M;
  real<lower=0> A1;
  real<lower=0> A2;
  real<lower=npi,upper=ppi> phi1;
  real<lower=npi,upper=ppi> phi2;
  real<lower=-1,upper=1> ar1;
  real err;
  real<lower=sigmaLow,upper=sigmaHigh> sigma;
  real<lower=0> nuMinusOne;
}
transformed parameters {
  real<lower=1> nu;
  real mu[N];
  nu <- nuMinusOne + 1;
  for (n in 2:N)
    mu[n] <- M + A1*cos(periodc*time[n]-phi1)
               + A2*cos(2*periodc*time[n]-phi2/2)
               + ar1*(x[n-1]-(M + A1*cos(periodc*time[n-1]-phi1)
                                + A2*cos(2*periodc*time[n-1]-phi2/2)));
}
model {
  M ~ normal(smean, sigmaHigh);
  sigma ~ uniform(sigmaLow, sigmaHigh);
  nuMinusOne ~ exponential(nuParm);
  ar1 ~ uniform(-1,1);
  A1 ~ uniform(0,rgx);
  A2 ~ uniform(0,rgx);
  phi1 ~ von_mises(0,0.0001);
  phi2 ~ von_mises(0,0.0001);
  err ~ student_t(nu,0,sigma);
  for (n in 2:N)
    x[n] ~ student_t(nu,mu[n],sigma);
}

Here, M is the intercept, A1 and A2 are the amplitudes for the 24- and 12-hour effects, respectively, phi1 and phi2 are to shift the cosinor component, ar1 is the AR(1) term, sigma is the variance of the error, and nu is the normality parameter of the error.

 

Most of the declared variables that are not data are for convenience--like 'periodc' is the periodic constant 2*pi because the time variable goes from 0 to 1, 'npi' and 'ppi' are -pi and pi, but I might want to shift the upper and lower bounds of the periodicity in case one of the 'phi's was close to the boundary.  The 'nuParm' variable is the parameter for the exponential prior on nu.  When I pass the data from R, it looks like this:

 

data_list = list(x = subdatb[,var],
                 time = subdatb$Daily,
                 periodc=2*pi,
                 npi=-pi,
                 ppi=pi,
                 sigmaLow=sd(subdatb[,var],na.rm=T)/10000,
                 sigmaHigh=10000*sd(subdatb[,var],na.rm=T),
                 N = nrow(subdatb),
                 smean=mean(subdatb[,var]),
                 nuParm=1/29,
                 rgx=diff(range(subdatb[,var])))

 

This is my stan call:

 

stan.fit <- stan(model_code = model_string, data = data_list,
                 iter = 5000, chains = 4)

 

This code worked fine, but for some subjects and some parameters there was enough skewness in the data that the t-distributed error model didn't appear to be robust enough for the data, so I fit it using the skew-normal distribution:

 

data {
  int<lower=0> N;
  real x[N];
  real time[N];
  real periodc;
  real npi;
  real ppi;
  real<lower=0> sigmaLow;
  real<lower=0> sigmaHigh;
  real smean;
  real<lower=0> rgx;
}
parameters {
  real M;
  real<lower=0> A1;
  real<lower=0> A2;
  real<lower=npi,upper=ppi> phi1;
  real<lower=npi,upper=ppi> phi2;
  real<lower=-1,upper=1> ar1;
  real err;
  real<lower=sigmaLow,upper=sigmaHigh> sigma;
  real skew;
}
transformed parameters {
  real mu[N];
  for (n in 2:N)
    mu[n] <- M + A1*cos(periodc*time[n]-phi1)
               + A2*cos(2*periodc*time[n]-phi2/2)
               + ar1*(x[n-1]-(M + A1*cos(periodc*time[n-1]-phi1)
                                + A2*cos(2*periodc*time[n-1]-phi2/2)));
}
model {
  M ~ normal(smean, sigmaHigh);
  sigma ~ uniform(sigmaLow, sigmaHigh);
  skew ~ normal(0,10000);
  ar1 ~ uniform(-1,1);
  A1 ~ uniform(0,rgx);
  A2 ~ uniform(0,rgx);
  phi1 ~ von_mises(0,0.0001);
  phi2 ~ von_mises(0,0.0001);
  err ~ skew_normal(0,sigma,skew);
  for (n in 2:N)
    x[n] ~ skew_normal(mu[n],sigma,skew);
}

Again, this code works fine but the model wasn't quite right for some of the data.

 

Up to this point, I've done all this using both rstan and rjags.  In rjags, I was able to fit a model that uses a skew-t error (using the zeros trick), and it usually fixed most of the problems I had with the fit using the other two error models.  When the fit still wasn't right, there were some mixing issues and I suspect we may have been hitting the wall on what jags could handle with the zeros trick, so I was hoping to be able to fit the model in stan to give the Hamiltonian Mone Carlo a chance at it, but I haven't been able to get the skew-t to work in stan.

 

In my stan version of the skew-t, I have everything the same from the models I've included above, except I'm combining the nu parameter from the t model and the skew parameter from the skew-normal model.  Essentially, I start with the t model and the increment the log likelihood to add the skewness.  I've also added a transformed variable, z, to increase readability:

 

data {
  int<lower=0> N;
  real x[N];
  real time[N];
  real periodc;
  real npi;
  real ppi;
  real<lower=0> sigmaLow;
  real<lower=0> sigmaHigh;
  real smean;
  real<lower=0> nuParm;
  real<lower=0> rgx;
}
parameters {
  real M;
  real<lower=0> A1;
  real<lower=0> A2;
  real<lower=npi,upper=ppi> phi1;
  real<lower=npi,upper=ppi> phi2;
  real<lower=-1,upper=1> ar1;
  real<lower=sigmaLow,upper=sigmaHigh> sigma;
  real skew;
  real<lower=0> nuMinusOne;
}
transformed parameters {
  real<lower=1> nu;
  real mu[N];
  real z[N];
  nu <- nuMinusOne + 1;
  for (n in 2:N) {
    mu[n] <- M + A1*cos(periodc*time[n]-phi1)
               + A2*cos(2*periodc*time[n]-phi2/2)
               + ar1*(x[n-1]-(M + A1*cos(periodc*time[n-1]-phi1)
                                + A2*cos(2*periodc*time[n-1]-phi2/2)));
    z[n] <- (x[n]-mu[n])/sigma;
  }
}
model {
  M ~ normal(smean, sigmaHigh);
  sigma ~ uniform(sigmaLow, sigmaHigh);
  skew ~ normal(0,10000);
  nuMinusOne ~ exponential(nuParm);
  ar1 ~ uniform(-1,1);
  A1 ~ uniform(0,rgx);
  A2 ~ uniform(0,rgx);
  phi1 ~ von_mises(0,0.0001);
  phi2 ~ von_mises(0,0.0001);
  for (n in 2:N){
    x[n] ~ student_t(nu,mu[n],sigma);
    increment_log_prob(-log(sigma)
                       + student_t_cdf_log(skew*z[n]*sqrt((nu+1)/(pow(z[n],2)+nu)),nu+1,0,1));
  }
}

This code does not work.  I get the following error:

 

Error : Error in function ibeta_derivative<e>(e,e,e): Overflow Error
error occurred during calling the sampler; sampling not done

The main difference between the working code and the non-working code is the increment_log_prob() function, so I'm hoping that I've implemented it wrong and that someone might have a simple fix for this.

 

The other question I have has to do with my 'err' parameter, which I'm using to put confidence bands on the model fit.  I don't know how to sample from the skew-t in stan for this parameter as I did with the others.  I can use the 'rst()' function from the 'sn' package in R, but it would be smoother if I could just do this from within stan as I did with the first two models.

 

[Michael Betancourt]

This issue refers to the fact that you're really pushing the numerical stability
of the Student CDF. Remember that Stan requires derivatives, and although
we have working derivatives for the CDFs they are not the most numerically
stable implementations, which are an incredibly difficult problem.

[Andrew Gelman]

Would it be better to model the t as a mixture, or would that create other problems?

> --
> 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/groups/opt_out.

[Michael Betancourt]

Anything not involving the CDF would be more better numerically, but may be worse in terms of sampling.

[Sam Dickson]

I'll see if implementing a mixture will fit the data better better than just the t or the skew-normal.  Now I'm wondering if I should attempt a mixture of t's, a mixture of skew-normals, or a t with a skew-normal?  Sounds like a can of worms, but at least I have somewhere to go.

 

Thanks.

[Sam Dickson]

Do you think this problem would persist if I were able to create a skew-t module based on the skew-normal?

On Thursday, December 19, 2013 11:02:46 AM UTC-5, Michael Betancourt wrote:

[Michael Betancourt]

The normal CDF is better-behaved numerically, but you still have to be careful with
extremely large values.

[Andrew Gelman]

Or even something based on a Cauchy?  At least this has a simple cdf.

[Bob Carpenter]

Working on cdfs, complementary cdfs, and their log versions
would be a great project for someone who wants to dive
deep into numerical processing. The hard part is getting
derivatives w.r.t. the parameters and stable complementary cdfs.

There are intrinsic problems with floating point for representing
cdfs near 1 --- the closest you can get is about 1e-14. That's
why you get both erf() and erfc() in numerical libs, and both log()
and log1m().

- Bob