Repository navigation
Interpolation function(s) #58
Description
Activity
From @bob-carpenter on December 10, 2014 18:12
Related to #731 (in a good way).
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.
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), wherealphais a parameter andx(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); }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
Reacted by Deemah- changed the title
[-]Interpolation methods in Stan?[/-][+]Interpolation function(s)[/+]on Aug 18, 2016 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.
Is anyone else working on this? I might wanna steal it if not.
Reacted by Steve Bronder and John HarringtonGo!
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.,
would result in
x_tequal to a value linearly interpolated fromx[i],x[i + 1], wherei <= t <= i + 1.Copied from original issue: stan-dev/stan#1165