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)
)