1-dimensional Interpolation Functions

Functions
array[] vector	interp_1d_cubic (array[] vector y, data array[] real x, data array[] real x_out, array[] int a0, vector theta)
array[] vector	interp_1d_linear (array[] vector y, data array[] real x, data array[] real x_out)

Detailed Description

These functions perform interpolation one dimension at a time.

 
array[] vector interp_1d_cubic(array[] vector y, data array[] real x,
                               data array[] real x_out, array[] int a0,
                               vector theta) {
  int left = 1;
  int right = 1;
  real h = 0.0;
  real w = 0.0;
  int N_in = size(x);
  int N_out = size(x_out);
  int D = size(y[1]);
  vector[D] f_left;
  vector[D] f_right;
  real h00;
  real h10;
  real h01;
  real h11;
  array[N_out] vector[D] y_out;
  for (j in 1 : N_out) {
    // Find left and right point indices
    while (x[right] < x_out[j]) {
      right = right + 1;
    }
    while (x[left + 1] < x_out[j]) {
      left = left + 1;
    }
    
    // Evaluate derivatives
    f_left = derivative_fun(0.0, y[left], a0, theta);
    f_right = derivative_fun(0.0, y[right], a0, theta);
    
    // Hermite basis functions
    h = x[right] - x[left];
    w = (x_out[j] - x[left]) / h;
    h00 = 2.0 * w ^ 3 - 3.0 * w ^ 2 + 1.0;
    h10 = w ^ 3 - 2.0 * w ^ 2 + w;
    h01 = -2.0 * w ^ 3 + 3.0 * w ^ 2;
    h11 = w ^ 3 - w ^ 2;
    
    // Compute interpolation
    y_out[j] = h00 * y[left] + h10 * h * f_left + h01 * y[right]
               + h11 * h * f_right;
  }
  return y_out;
}
 

 
array[] vector interp_1d_linear(array[] vector y, data array[] real x,
                                data array[] real x_out) {
  int left = 1;
  int right = 1;
  real w = 1.0;
  int N_in = size(x);
  int N_out = size(x_out);
  int D = size(y[1]);
  array[N_out] vector[D] y_out;
  for (j in 1 : N_out) {
    while (x[right] < x_out[j]) {
      right = right + 1;
    }
    while (x[left + 1] < x_out[j]) {
      left = left + 1;
    }
    w = (x[right] - x_out[j]) / (x[right] - x[left]);
    y_out[j] = w * y[left] + (1 - w) * y[right];
  }
  return y_out;
}

Function Documentation

interp_1d_cubic()

array[] vector interp_1d_cubic	(	array[]vector	y,
		data array[]real	x,
		data array[]real	x_out,
		array[]int	a0,
		vector	theta
	)

1d interpolation using cubic Hermite splines

The interpolated value \(y = g(x)\) is given by a cubic polynomial \(g\), defined uniquely by four values \(y_0 = g(x_0)\), \(y_1 = g(x_1)\), \(f_0 = g'(x_0)\), and \(f_1 = g'(x_1)\), where \(x_0 < x \leq x_1\). Assumes that derivative_fun(), which evaluates \(g'(x)\), is defined in th functions block before this function. It must have signature

vector derivative_fun(real t, vector y, int[] a0, vector a1)

i.e. same as what can be used with the Stan ODE solvers but it shouldn't actually use the t argument. It will be always called with t = 0.0.

Info: https://en.wikipedia.org/wiki/Cubic_Hermite_spline

Author

Juho Timonen

Parameters

y	array of D-vectors, length N_in
x	incresing array of reals, length N_in
x_out	increasing array of reals, length N_out, values must be in \((\min(x), \max(x)]\)
a0	array of integer inputs given to derivative_fun()
theta	parameter vector given to derivative_fun()

Returns

array of D-vectors, length N_out, corresponding to interpolated values y_out

interp_1d_linear()

array[] vector interp_1d_linear	(	array[]vector	y,
		data array[]real	x,
		data array[]real	x_out
	)

Linear 1d interpolation

Author

Juho Timonen

Parameters

y	array of D-vectors, length N_in
x	incresing array of reals, length N_in
x_out	increasing array of reals, length N_out, values must be in \((\min(x), \max(x)]\)

Returns

array of D-vectors, length N_out, corresponding to interpolated values y_out