2D interpolation and splines

[Jóhannes Guðbrandsson]

Hello

Is there a easy way to use 2D interpolation in a stan-model? If not is that on the agenda?

Same with splines. I did not find examples in the manual. But I was wondering is someone had any examples before I would dig into it myself.

Cheers

Jóhannes

[Krzysztof Sakrejda]

Splines work fine if you calculate the weights in Stan or preferably just construct the model matrix for them outside of Stan.  There are some good functions in R for making the model matrix (if that's the interface you use).  Krzysztof

[Jóhannes Guðbrandsson]

I have tried that in JAGS. I was just wondering if there was some fancy (easy) way someone had already implemented.

I am mainly interested in the interpolation part though. I haven't seen a way to get that to work. 

[Jeff]

Hi,

I was playing around with fitting splines in rstan, and using them to generate predictions. 

This might be helpful example:

https://gist.github.com/paleo13/d4d15eedd40abcbeb0a0

Cheers,

Jeff

[Jan]

What R functions for making the model matrix can you recommend?

J

[Bob Carpenter]

I'm confused --- this looks like an ordinary linear
regression.

And for prediction, you almost certainly are going to
want the noise


y_predict <- normal_rng(X_predict * beta, sigma);

otherwise, you'll underestimate prediction uncertainty.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Krzysztof Sakrejda]

On Wednesday, August 26, 2015 at 12:53:11 PM UTC-4, Bob Carpenter wrote:

I'm confused --- this looks like an ordinary linear
regression.

I imagine the spline is in the X matrix.  If you have fixed knot point locations then the spline doesn't need to be calculated in Stan. K
 

[Michael Ciere]

Hi Jóhannes,

I did some experimenting with splines in Stan -- specifically with Ramsay's M-splines (non-negative basis splines that integrate to 1) and I-splines (integrated M-splines, useful for fitting a monotone spline). The idea is the same if you use B-splines though. I found this to work very well with Stan.

As Krzysztof said, you can precompute the values of the basis splines outside of Stan and pass them through as data, and then the only model parameters are the weights of these basis splines. You probably want to use a prior for the spline weights that induces a smooth spline though, like the Bayesian P-spline specification described by Lang & Brezger (http://dx.doi.org/10.1198/1061860043010). That is, in the github example posted here you may replace 

beta ~ normal(0,10) 

by something like 

tail(beta, d-1) ~ normal(head(beta, d-1), sigma_beta) 

with some prior for sigma_beta.

If there is interest I can share some of my stuff and help construct an illustrative example -- I do think that Stan excels in this kind of thing.

Regards,

Michael

[Jeffrey Hanson]

Thanks for clarifying the normal_rng.

 

Yeah, the knots are at fixed positions. I used R’s model.matrix function and splines::ns to generate a design matrices for the training and prediction datasets. I just used stan to estimate the parameters for the pre-constructed design matrix for the training data, so the stan code I used looks very similar to that for a regression.

 

I was playing around with the code and worked out how to create interpolations beyond the range of the training dataset,

 

The updated codes here if anyone’s interested,

 

https://gist.github.com/paleo13/d4d15eedd40abcbeb0a0

[Bob Carpenter]

You have to do interpolation manually. Writing the interpolate
function from BUGS has been our to-do list. It's non-trivial to
do remotely because we don't have integer rounding functions to
compute indices.

- Bob

[Jóhannes Guðbrandsson]

Thanks everybody for your replies. 

I finally got some time to work on this again.

I wrote this function. It maybe not the most efficient but it seem to work. Any ideas from improvement are most appreciated.

  real interpol( real x_new, real y_new, vector x, vector y, vector Z){

    #Z is in column major order

    int i;

    int j;

    real deno;

    real nume;

    real large;

    real run;

    int Nx;

    int Ny;

    Nx <- num_elements(x);

    Ny <- num_elements(y);

    large <- x[Nx];

    i <- 2;

    j <- 2;

    run <- x[i];

    if (x_new > large)  i <- Nx;

    else while (x_new > run){

i <- i+1;

run <- x[i];

      }

    

    large <- y[Ny];

    run <- y[j];

    if (y_new > large) j <- Ny;

    else while (y_new > run){

j <- j+1;

run <- y[j];

      }

    

    deno <- ( x[i] - x[i-1] ) * ( y[j] - y[j-1] );

    nume <- Z[i-1+(j-1-1)*Nx] * (x[i]-x_new) * (y[j]-y_new) +

      Z[i+(j-1-1)*Nx] * (x_new-x[i-1]) * (y[j]-y_new) +

      Z[i-1+(j-1)*Nx] * (x[i]-x_new) * (y_new-y[j-1]) +

      Z[i+(j-1)*Nx] * (x_new-x[i-1]) * (y_new-y[j-1]);

    return nume/deno;

[Bob Carpenter]

That looks OK (other than indentation/spacing).

These two are the same, right:

Z[i+(j-1-1)*Nx]

Z[i-1+(j-1)*Nx]

We might try to add something built in for this. It's been
our to-do list but nobody's wanted to figure out the right interface
and implement it.

- Bob

> To post to this group, send email to stan-...@googlegroups.com.

[Jóhannes Guðbrandsson]

No. The first one is Z[ i, j-1] and the second Z[ i-1, j].