I’ve been struggling with modelling simple periodic functions, discovering that the likelihood topography is massively multimodal in very unintuitive ways. But recently I had the idea to reframe things with a hidden markov model applied to the slope of the function derivative. I’ve had some seeming success with the model below, but I’m having trouble working out how to take the successful model of the derivative and integrate to get back to the original function.

In fact, ultimately I’ll be looking to do a whole lot more modelling of the observations, including things like the observations being the additive combination of a known number of periodic functions with well-separated frequencies, a covariate that additionally influences the observations, etc. So I’m thinking that maybe I’m going to have to abandon use of `hmm_marginal`

anyway, and maybe that’s what I need to do here to do inference on the integrated function?

Here’s an R script to play. It’ll produce data that looks like (dots are data, line is true function):

And it’ll show that the model does a great job of classifying the slopes:

But my hamfisted attempts to work back to the original function are clearly not working:

And here’s the model itself:

```
// derivative_hmm.stan
data{
int n ;
vector[n] dy ; // derivative (lag-1 difference) of y
real y0 ; // zeroth y observation (for attempting to align f)
}
parameters{
real<lower=0> noise ; //measurement noise
real<lower=0> pos_slope ; // slope for positive slope state
real<upper=0> neg_slope ; // slope for negative slope state
simplex[2] p0 ; // initial state probabilities
simplex[2] tp_pn ; // transition probability for positive to negative
simplex[2] tp_np ; // transition probability for negative to positive
}
transformed parameters {
matrix[2, 2] tps = transpose(append_col(tp_pn,tp_np));
matrix[2, n] state_ll ;
// Compute the log likelihoods in each possible state
for(i_n in 1:n) {
state_ll[1, i_n] = normal_lpdf(dy[i_n] | pos_slope, noise);
state_ll[2, i_n] = normal_lpdf(dy[i_n] | neg_slope, noise);
}
}
model{
noise ~ std_normal() ;
pos_slope ~ std_normal() ;
neg_slope ~ std_normal() ;
target += hmm_marginal(state_ll, tps, p0);
}
generated quantities{
// get the probability of a positive slope at each timepoint
row_vector[n] p_pos = hmm_hidden_state_prob(state_ll, tps, p0)[1,];
// compute the most probable slope at each timepoint
// vector[n] df = p_pos > .5 ? pos_slope : neg_slope ; // doesn't work; conditionals aren't vectorized
vector[n] df ;
for(i_n in 1:n){
if(p_pos[i_n]>.5){
df[i_n] = pos_slope ;
}else{
df[i_n] = neg_slope ;
}
}
// compute the forward and backward integrals
vector[n] fi ;
vector[n] bi ;
for(i_n in 1:n){
fi[i_n] = sum(df[1:i_n]);
bi[n-i_n+1] = sum(-df[(n-i_n+1):n]);
}
// average
vector[n] f = (fi+bi)/2 + y0 ;
}
```