Modelling a periodic response

Modeling
1

I have some photoplethysmography data where the heart beat causes a quasi-regular pattern in the data (image). There are a variety of features of a given pulsation wave form (ex. amplitude, width, etc) that I’m interested in modelling (eventually as varying auto-regressively through time, but for this post let’s focus on simply assuming they’re constant for a given data set). If I had a single pulsation, I’d define a function that models a single pulse with arguments for the parameters I’m interested in, yielding the height of the wave for a given time:

functions{
  real get_pulse(real x, real amplitude, real width){
    ...
  }
}
data{
    int n;
    vector[n] y;
    vector[n] x;
}
parameters{
    real<lower=0> amplitude;
    real<lower=0> width;
    real lower=0> noise;
}
model{
    .... \\ priors here
  y ~ normal(get_pulse(x,amplitude,width),noise);
}

But I’m not sure how I’d go about modelling a series of pulses. Presumably, as a periodic phenomenon, I should be using the unit_normal variable type somehow to reflect that at a certain value for x we’re wrapping back around to the beginning of a new pulse, but I’m not sure how to achieve this. Help?

2

Would a Gaussian process with a kernel that imposes periodicity work? e.g., eq. 20 in http://www.robots.ox.ac.uk/~sjrob/Pubs/philTransA_2012.pdf

2 Likes
3

Thanks for the idea, but I feel like a GP with a periodic kernel isn’t a good fit for this scenario as I have very good expectation for the precise shape that is going to be periodic and I want to do direct inference on the parameters of this shape. With enough data a GP with periodic kernel would figure out the periodic shape by itself, but do so only with a large compute demand and wouldn’t permit direct inference on the shape parameters.

4

@Bob_Carpenter Any thoughts on this? I feel like the solution is going to be something simple that I’m missing.

5

Ah! I think the solution involves fmod where I’d compute the remainder from dividing the time (x in the code above) by the period length (width in the code above), then dividing said remainder by the period length to get the position in the period (on a scale from zero to 1). Off to implement, will report back…

(n.b. as per the warning at the beginning of 41.7, I’m expecting poor sampling efficiency by using fmod, but hopefully not too bad…)

6

Yup, the fmod approach seemed to do the trick, and even sampled seemingly fairly efficiently despite the warning in section 41.7.

Here’s a gist of an example: https://gist.github.com/mike-lawrence/67d1f7cbf9fafc56a20339af80ef13d3

2 Likes
7

I just saw this now. I’m the wrong person to be asking about GPs!