I have noisy measurements of acceleration and position and want to obtain corrected values for acceleration, velocity, and position. To do this I express acceleration using splines, and integrate it in time to obtain velocity and position.

The basic idea is the same as for a Kalman filter.

The method works well in principle, a simple example result is shown below.

The main problem I have is that Stan consistently needs very high tree depths when solving this model. In the small example I show it is 11, which I could live with, but in some instances a tree depth of 15 is not enough.

How can debug this issue in a structured way and reparametrize the model to make it more computationally tractable? I have done a lot of trial and error, including following the Gaussian Process example in the user manual, but everything Iâ€™ve tried so far shows the same problem.

Any other hints how I can improve the model would of course also be appreciated. Are there other people out there who use Stan for similar applications?

```
functions {
// b-spline code from https://github.com/milkha/Splines_in_Stan
real build_b_spline(real t, vector ext_knots, int ind, int order);
real build_b_spline(real t, vector ext_knots, int ind, int order) {
// INPUTS:
// t: the points at which the b_spline is calculated
// ext_knots: the set of extended knots
// ind: the index of the b_spline
// order: the order of the b-spline
real b_spline;
real w1 = 0;
real w2 = 0;
if (order==1)
// B-splines of order 1 are piece-wise constant
b_spline = (ext_knots[ind] <= t) && (t < ext_knots[ind+1]);
else {
if (ext_knots[ind] != ext_knots[ind+order-1])
w1 = (t - ext_knots[ind]) /
(ext_knots[ind+order-1] - ext_knots[ind]);
if (ext_knots[ind+1] != ext_knots[ind+order])
w2 = 1 - (t - ext_knots[ind+1]) /
(ext_knots[ind+order] - ext_knots[ind+1]);
// Calculating the B-spline recursively as linear interpolation of two lower-order splines
b_spline = w1 * build_b_spline(t, ext_knots, ind, order-1) +
w2 * build_b_spline(t, ext_knots, ind+1, order-1);
}
return b_spline;
}
vector sho(real t, vector y, int num_basis, vector eta, vector ext_knots, int order) {
vector[2] dydt;
dydt[1] = 0.;
for (i in 1:num_basis) dydt[1] += build_b_spline(t, ext_knots, i, order)*eta[i];
dydt[2] = y[1];
return dydt;
}
}
data {
// signal
int<lower=0> N;
vector[N] time;
vector[N] altitude;
vector[N] acceleration;
// spline parameters
int num_knots; // num of knots
vector[num_knots] knots; // the sequence of knots
int spline_degree; // the degree of spline (is equal to order - 1)
}
transformed data {
int order = spline_degree + 1;
int num_basis = num_knots + spline_degree - 1; // total number of B-splines
vector[spline_degree + num_knots] ext_knots_temp;
vector[2*spline_degree + num_knots] ext_knots; // set of extended knots
ext_knots_temp = append_row(rep_vector(knots[1], spline_degree), knots);
ext_knots = append_row(ext_knots_temp, rep_vector(knots[num_knots], spline_degree));
// Build spline matrix at measurement points
matrix[N,num_basis] B;
for (i in 1:N) {
for (j in 1:num_basis) {
B[i,j] = build_b_spline(time[i], ext_knots, j, order);
}
}
vector[N-1] dt = time[2:N] - time[1:(N-1)];
}
parameters {
real v0;
real<lower=0> h0;
vector[num_basis] eta; // for a
real<lower=0> sigma_a;
real<lower=0> sigma_h;
}
transformed parameters {
vector[2] y_init = [v0, h0]';
vector[N] a_est;
profile("splines") {
a_est = B*eta;
// for (t in 1:N) {
// a_est[t] = 0.;
// for (i in 1:num_basis) a_est[t] += build_b_spline(time[t], ext_knots, i, order)*eta[i];
// }
}
vector[N] v_est;
vector[N] h_est;
v_est[1] = y_init[1];
h_est[1] = y_init[2];
profile("time_integration") {
for (i in 1:(N-1)) {
v_est[i+1] = v_est[i] + dt[i]*a_est[i];
h_est[i+1] = h_est[i] + dt[i]*v_est[i];
}
}
}
model {
profile("priors") {
y_init[1] ~ normal(0,10);
y_init[2] ~ normal(0,1000);
eta ~ normal(0,10);
sigma_a ~ exponential(.1);
sigma_h ~ exponential(.01);
}
profile("measurements") {
acceleration ~ normal(a_est, sigma_a);
altitude ~ normal(h_est, sigma_h);
}
}
```