Ok. Here’s slightly refactored vector-based function inspired by the PR#18 in stan-docs repo. It implements 3 different weight functions, as suggested by @nhuurre . This is linear interpolation, but I also have quadratic envelope method, which is even smoother, but it sacrifices the two outer bins.
functions{
real linear_interpolation_v(real x, vector x_pred, vector y_pred){
int K = rows(x_pred);
vector[K] deltas = x - x_pred;
real ans;
real t;
real w;
int i;
if(x<x_pred[1] || x>x_pred[K]) reject("x is outside of the x_pred grid!");
if(rows(y_pred) != K) reject("x_pred and y_pred aren't of the same size");
//this is which.min()
i = sort_indices_asc(fabs(deltas))[1];
if(deltas[i]<=0) i -= 1;
ans = y_pred[i];
real x1 = x_pred[i];
real x2 = x_pred[i + 1];
real y1 = y_pred[i];
real y2 = y_pred[i + 1];
ans = y1 + (y2-y1) * (x-x1)/(x2-x1);
t = (x-x1)/(x2-x1);
w = 1-t;
//w = 1/(1 + exp(1/(1-t) - 1/t));
//w = 1 - 3*pow(t,2) + 2*pow(t,3);
ans = w*y1 + (1-w)*y2;
return ans;
}
} // end functions block
Test it with the following R code
library(tidyverse)
set.seed(42) #set the correct seed
d <- tibble(
xs = 0:5,
ys=sample(10:15, 6, replace = TRUE),
type="known")
tibble(
xs=seq(0,5, by=0.001),
ys=sapply(xs, linear_interpolation_v, d$xs, d$ys),
type="predicted"
) %>%
bind_rows(d) %>%
ggplot()+
geom_point(aes(xs, ys, color=type))+
labs(title="w = ...")
