Monotonic splines in mcgv require “only” sampling fullfilling a constraint matrix C.
library(mgcv)
library(MASS)
library(rstan)
# x <- seq(0,1,length.out=10)
# y <- +0.2 * x + x^2*0.3 + 2
x <- 1:10
y=c(100,41,22,10,6,7,2,1,3,1)
sigma <- 0.05
y.noise <- y + rnorm(length(x), 0,sigma)
df <- data.frame(x=x, y=y.noise)
## Set up the size of the basis functions/number of knots
k <- 5
## This fits the unconstrained model but gets us smoothness parameters that
## that we will need later
unc <- gam(y ~ s(x, k = k, bs = "cr"), data = df)
## This creates the cubic spline basais functions of `x`
## It returns an object containing the penalty matrix for the spline
## among other things; see ?smooth.construct for description of each
## element in the returned object
sm <- smoothCon(s(x, k = k, bs = "cr"), df, knots = NULL)[[1]]
## This gets the constraint matrix and constraint vector that imposes
## linear constraints to enforce montonicity on a cubic regression spline
## the key thing you need to change is `up`.
## `up = TRUE` == increasing function
## `up = FALSE` == decreasing function (as per your example)
## `xp` is a vector of knot locations that we get back from smoothCon
K <- 5
df.sm=data.frame(x=x)
sm <- smoothCon(s(x, k = K, bs = "cr"), df.sm, knots = NULL
, absorb.cons=FALSE # no mean in spline. add. parameter in Stan
, diagonal.penalty = TRUE # diagonal smooth matrix
)[[1]]
Nzero <- sum(diag(sm$S[[1]]) == 0.0)
F <- mono.con(sm$xp, up = FALSE) # get constraints: up = FALSE == Decreasing constraint!
C<-F$A
The Question is now how can we sample parameters fulfilling the constraints in matrix C?
Maybe we can extend @bgoodri Constrained multinormal?