[EXAMPLES] Add example for tensor splines

[schalkdaniel]

This would be pretty nice, calculate a bivariate function based on the tensor spline.

[schalkdaniel]

nsim = 5000L

x1 = runif(nsim, 0, 5)
x2 = runif(nsim, 0, 5)

y = sin(x1) + cos(x2) + rnorm(nsim) 

# Grid for plotting:
x1pred = seq(min(x1), max(x1), length.out = 100)
x2pred = seq(min(x2), max(x2), length.out = 100)
df_plt = expand.grid(x = x1pred, y = x2pred)

df_plt$z = sin(df_plt$x) + cos(df_plt$y)# + rnorm(nrow(df_plt))

library(ggplot2)

gg1 = ggplot() +
  geom_contour_filled(data = df_plt, aes(x = x, y = y, z = z), bins = 15) +
  ggtitle("sin(x) + cos(y)")
#  geom_rug(data = dat0, aes(x = age, y = capital.gain))


predictTensor = function(coef, x1, x2, knots1, knots2, degree = 3) {
  basis1 = cpsp::createSplineBasis(values = x1, degree = degree, knots = knots1)
  basis2 = cpsp::createSplineBasis(values = x2, degree = degree, knots = knots2)
  
  tensor = cpsp::rowWiseTensor(basis1, basis2)
  return(tensor %*% coef)
}
myEstimator = function(X, y, penmat = 0, xtx = NULL, xty = NULL) {
  if (! (missing(X) || missing(y))) {
    xtx = crossprod(X)
    xty = crossprod(X, y)
  }
  L = chol(xtx + penmat)
  z = backsolve(L, xty, transpose = TRUE)
  return(as.vector(backsolve(L, z)))
}

# Create spline bases + tensor product:
k1 = cpsp::createKnots(x1, 10, 3)
k2 = cpsp::createKnots(x2, 10, 3)
b1 = cpsp::createSplineBasis(x1, 3, k1)
b2 = cpsp::createSplineBasis(x2, 3, k2)

bt = cpsp::rowWiseTensor(b1, b2)
xtxt = t(bt) %*% bt

p1 = cpsp::penaltyMat(ncol(b1), 3)
p2 = cpsp::penaltyMat(ncol(b2), 3)

pt = kronecker(p1, p2)
lambda = cpsp::demmlerReinsch(xtxt, pt, 50)

betat = myEstimator(bt, y, penmat = lambda * pt)

df_plt$pred = predictTensor(betat, df_plt$x, df_plt$y, k1, k2)

gg2 = ggplot() +
  geom_contour_filled(data = df_plt, aes(x = x, y = y, z = pred), bins = 15) +
  geom_rug(data = data.frame(x = x1, y = x2), aes(x = x, y = y)) +
  ggtitle("Bivariate spline prediction")

gridExtra::grid.arrange(gg1, gg2, ncol = 2L)