Inference on the rate parameter of a periodic response

Modeling
1

I’m looking to model some physiological data consisting of a periodic response induced by the heart’s pulsation. I’ve coded up a simple example here that achieves inference on the magnitude of the response and noise, but assumes a known constant pulse rate.

Here’s a version of the model that attempts to additionally do inference on the rate parameter:

functions{
  // pulse: function to generate pulse-like curve; oversimplified in shape for initial model building/testing
  vector pulse(vector time, real amplitude, real rate){
    vector[num_elements(time)] position ;
    vector[num_elements(time)] beta ;
    real period = 1/rate ;
    for(i in 1:num_elements(time)){
      position[i] = fmod(time[i],period)/period ; //this needs to be looped as fmod isn't vectorized
      beta[i] = beta_lpdf(position[i]|2,2) ; // this needs to be looped as beta_lpdf will sum it result to a real if done on a vector
    }
    return amplitude * exp(beta) ;
  }
}
data{
  int n;
  vector[n] y;
  vector[n] time;
}
parameters{
  real<lower=0> amplitude;
  real<lower=0> rate;
  real<lower=0> noise;
}
model{
  rate ~ weibull(2,2) ;
  amplitude ~ weibull(2,2) ;
  noise ~ normal(0,1) ;
  y ~ normal( pulse(time,amplitude,rate) , noise ) ;
}

However, with this version Stan can’t initialize, producing the error Rejecting initial value: Gradient evaluated at the initial value is not finite. Stan can't start sampling from this initial value. My usual solution to initialization troubles is to set init=0 in rstan::sampling(), but this wouldn’t make sense in this scenario since zero makes little sense for the amplitude and rate parameters (also, just checked and init=0 & init=1 both yield the same error message as above). Any thoughts on how to get this to initialize properly?

2

AFAIK the zero init refers to zero in unconstrained parameter space. You could also try narrowing down the range of initial values (default is -2,2) instead of setting it to zero across all chains. I think there is the init_r option for that.

3

The issue is probably beta_lpdf giving -Inf when it gets 0 or 1 as argument. For example, changing the line to the following seems to work (though this is probably not the best way to deal with it, but shows that that is probably the issue):

  position[i] = (fmod(time[i],period) + 1e-22)/(period + 1e-22) ; //this needs to be looped as fmod isn't vectorized
1 Like
4

Good catch! Unfortunately the model now samples very poorly; fast (iter=1e5 takes a couple mins) but with ridiculously low effective sample sizes and high rhats. I chose the beta as a shape for the periodic response rather arbitrarily, so I’ll try a few different shapes to see if this is an issue with the beta, or (as I pessimistically expect) with discontinuities in the gradient induced by using fmod.

5

Still working on investigating poor sampling performance, but thought I’d note that the issue with the beta can be avoided by instead using a sine:

// pulse: function to generate pulse-like curve; oversimplified in shape for initial model building/testing
functions{
  vector pulse(vector time, real amplitude, real rate){
    vector[num_elements(time)] position ;
    real period = 1/rate ;
    for(i in 1:num_elements(time)){
      position[i] = fmod(time[i],period)/period ; //this needs to be looped as fmod isn't vectorized
    }
    return amplitude*exp(sin(position*pi())) ;
  }
}