Fitting a sine function by specifying a circular predictor

Interfaces brms
,
1

I have some measurements (y) that follow a sine function with a known frequency. I would like to fit a non-linear model to estimate the amplitude and the phase shift. Is it possible to specify that the phase shift (ps in the below example) predictor is circular (0, 2 π)?

library(brms)

# Fake data
N <- 360 # number of measurements
x <- seq(0, 2*pi, length.out = N)
e <- rnorm(n = N, mean = 0, sd = 5) # noise
amp <- 10 # amplitude
freq <- 3 # frequency
ps <- pi/6 # phase shift
y <- amp * sin(freq * x - ps)) + e
mydata <- data.frame(x, y)

# Priors
priors <- prior(normal(0, 15), nlpar = "amp", lb = -25, ub = 25) +
  prior(cauchy(0, 2), class = 'sigma')

# Current model
mymodel <- brm(
  bf(y ~ amp * sin(3 * x - ps), amp + ps ~ 1, nl = TRUE),
  data = mydata, 
  prior = priors,
  warmup = 1000,
  iter = 4000,
  chains = 4,
  cores = 4
)
2

I have no idea about brms, but you can do this in Stan itself by declaring the variable ps to be of type unit_vector[2], then computing the angle from (x, y) coordinates. The problem with using a real<lower=0, upper= 2*pi()> variable is that the 0 and 2 * pi are on either ends of the interval, but are the same point. So if mass bunches up near 0, the posterior using a single real will be bimodal. The unit_vector[2] is overparameterized, but at least it won’t lead to a multimodal posterior.

1 Like
3

If you know the cycle length, then including a cosine term in addition to the sine term effectively estimates the phase shift. More here.

4

Just to confirm Bob’s response that I think there is currently no way to make this happy in brms directly. You may let brms generate the Stan code and then edit the code in the right places following Bob’s suggestions. Whether this is faster in your case than writing the Stan code from scratch I cannot tell.

5

Thanks for the suggestion. I have tried to implement it in Stan by declaring the variable ps to be of type unit_vector[2] (see below), but although the estimates look alright on simulated data, the fit leads to quite a large number of divergent transitions. I tried to increase adapt_delta and max_treedepth, but it didn’t help much.

Do you have any improvement to suggest to the below Stan code to avoid the divergent transitions?

Stan code:

data {
  int<lower=1> N;         // total number of observations
  vector[N] x;            // x values
  vector[N] y;            // response variable
  real<lower=0> freq;     // fixed frequency of the sine wave
}

parameters {
  real<lower=0,upper=20> A; // amplitude
  unit_vector[2] phi_vec; // phase shift represented as a unit vector
  real<lower=0> sigma;    // dispersion parameter
}

transformed parameters {
  real phi;               // phase shift
  phi = atan2(phi_vec[2], phi_vec[1]); // convert unit vector to phase in radians
}

model {
  vector[N] mu;           
  
  // priors
  A ~ normal(0, 15);   // prior for amplitude 
  phi ~ normal(0, 1);  // prior for phase shift
  sigma ~ normal(0, 2);   // prior for noise
  
  // Likelihood
  for (n in 1:N) {
    mu[n] = A * sin(freq * x[n] - phi);
  }
  y ~ normal(mu, sigma);  // Gaussian likelihood with noise

}

R code:

library(rstan)

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

N <- 400 # number of measurements
x <- seq(0, 2*pi, length.out = N)
e <- rnorm(n = N, mean = 0, sd = 5) # noise
amp <- 10 # amplitude
freq <- 3 # frequency
ps <- 1/400*(2*pi) # phase shift
y <- amp * sin(freq * x - ps) + e
mydata <- data.frame(x, y)

data_list <- list(
  N = length(mydata$x),   # Number of observations
  x = mydata$x,           # Vector of x values
  y = mydata$y,           # Vector of observed values
  freq = 3               # Fixed frequency
)

fit <- stan(
  file = "sine-model.stan",
  data = data_list,
  chains = 4,             
  warmup = 2000,          
  iter = 8000,            
  cores = 4,              
  control = list(adapt_delta = 0.999, max_treedepth = 25)
)
6

I’m no expert but I guess the divergences derive from numerical issues in the evaluation of atan2(phi_vec[2], phi_vec[1]) when phi_vec[1] passes through zero, because atan2 is based on the ratio of its two arguments. (Elsewhere I’ve noticed divergence transitions when I calculate quantities that may explode, even if these are not used to update target.)
When I replace

unit_vector[2] phi_vec; 
...
phi = atan2(phi_vec[2], phi_vec[1]);

in the above code by

simplex[2] phi_vec_sqr;
...
phi = atan2(sqrt(phi_vec_sqr[2]), sqrt(phi_vec_sqr[1]));

the divergences disappear. This way prevents the denominator in the ratio inside atan2 from passing through zero. But it only allows you to estimate values of phi in [0, pi/2], and I’ve seen it behave weirdly at the boundaries of that range. And if you can exclude some possible values of phi a priori, you could avoid the circular parametrisation in the first place by just using e.g. real<lower = 0, upper = pi()/2> phi, knowing that this won’t be multimodal due to circularity.
The innards of atan2 don’t make sense to me; if the problem is as I suspect, maybe there’s a way to make it numerically robust as its second argument’s magnitude gets dangerously small?