Interpolation function(s)

[syclik]

From @edbaskerville on December 9, 2014 23:33

It would be nice if Stan had some simple one-dimensional interpolation methods built in. (Not super-important, since this should be doable within a Stan model.)

One place this issue arises when implementing ODE models driven by empirical data sampled at discrete time intervals. Because the time derivative needs to be evaluated at arbitrary time points, the values of the driving variables need to be evaluated between sample times.

Even linear interpolation would be sufficient for many cases: e.g.,

data {
  int nTimepoints;
  real[nTimepoints] x;
}
// ...
x_t <- interpolate_linear(x, t);

would result in x_t equal to a value linearly interpolated from x[i], x[i + 1], where i <= t <= i + 1.

Copied from original issue: stan-dev/stan#1165

[syclik]

From @bob-carpenter on December 10, 2014 18:12

Related to #731 (in a good way).

[syclik]

From @bob-carpenter on December 10, 2014 18:13

We don't really need that in the diff eq solver. It already uses interpolation internally, with the added bonus of error bounding (in case linear interpolation isn't accurate enough).

But there are certainly other uses.

[syclik]

From @edbaskerville on December 11, 2014 0:25

Cool, sounds like interpolation based on smooth functions is on the way, glad to hear it.

For ODEs—and perhaps I'm not following—I wasn't talking about the solver per se, but rather the use of an interpolated function inside the time derivative (making it a non-autonomous system). That is dy/dt (say, infection rate) is a function of an interpolated x(t) (which is calculated from an observed, discrete time series). E.g., dy/dt = alpha * x(t), where alpha is a parameter and x(t) is an interpolated function.

In the meantime, I went ahead and implemented linear interpolation between adjacent points in an array directly in Stan. It's rather horrible, because I had to sneak around the no-real-to-int-conversion restriction by implementing an integer-floor function via bisection.

For my purposes (using an interpolated x(t) in an ODE), I think this still is OK because the partial derivatives of dy/dt with respect to the parameters should treat values of x(t) as constants (or perhaps I'm missing something?)—and the function, though not differentiable, is continuous, and thus can be integrated.

Here it is, in its mostly untested, frightening glory:

int intFloor(int leftStart, int rightStart, real iReal)
{
  // This is absurd. Use bisection algorithm to find int floor.
  int left;
  int right;

  left <- leftStart;
  right <- rightStart;

  while((left + 1) < right) {
    int mid;
    // print("left, right, mid, i, ", left, ", ", right, ", ", mid, ", ", iReal);
    mid <- left + (right - left) / 2;
    if(iReal < mid) {
      right <- mid;
    }
    else {
      left <- mid;
    }
  }
  return left;
}

// Interpolate arr using a non-integral index i
// Note: 1 <= i <= length(arr)
real interpolateLinear(real[] arr, real i)
{
  int iLeft;
  real valLeft;
  int iRight;
  real valRight;

  // print("interpolating ", i);

  // Get i, value at left. If exact time match, then return value.
  iLeft <- intFloor(1, size(arr), i);
  valLeft <- arr[iLeft];
  if(iLeft == i) {
    return valLeft;
  }

  // Get i, value at right.
  iRight <- iLeft + 1;
  valRight <- arr[iRight];

  // Linearly interpolate between values at left and right.
  return valLeft + (valRight - valLeft) * (i - iLeft);
}

[bob-carpenter]

Copied from #64:

From @bob-carpenter on January 12, 2015 19:20

Jan Scholz (via e-mail) suggested implementing the general form of interpolation used in SciPy.

https://github.com/scipy/scipy/blob/v0.15.0/scipy/ndimage/interpolation.py#L230

It should be able to handle 2D and 3D structures. Jan says the linear form (order 1) should be enough for his purposes.

For background, BUGS provides an interp.lin function for the 1D case:

http://www.mrc-bsu.cam.ac.uk/wp-content/uploads/manual14.pdf

The goal would be to get the derivatives working properly in a way that encapsulates the index-fiddling relative to the user and takes care of error conditions.

This would be a nice standalone project.

• example for model-examples
• documentation using example
• C++ function with tests
• plumbing into Stan's language
• tests in Stan programs to ensure compilability

Copied from original issue: stan-dev/stan#1214

[bob-carpenter]

We need a concrete plan for what one of these functions should look like to be put into the issue itself with some description of how the derivatives work.

[bbbales2]

Is anyone else working on this? I might wanna steal it if not.

[wds15]

Go!