Well, 2.5 years later…I was interested in monotone splines and did the polar decomposition. No divergences or treedepth issues with the simulated code. ESS is way higher too.
The polar decomposition can be accomplished with SVD. What’s really cool is that the rc of Stan now has SVD in it! I pass the P and U matrices to Stan.
library(splines2)
library(cmdstanr)
set.seed(42)
tmin <- 0
tmax <- 10
Npoints <- 100
ninterior_knots <- 5
sigma_noise <- .1
times <- seq(tmin,tmax,length.out = Npoints)
knots <- seq(tmin, tmax, length.out = ninterior_knots+2)[2:ninterior_knots+1]
nknots <- length(knots)
mspline_degree<- 3
i_spline_basis_evals <- iSpline(times, knots=knots, degree=mspline_degree+1,
intercept=TRUE)
m_spline_basis_evals <- deriv(i_spline_basis_evals)
nbasis <- dim(m_spline_basis_evals)[2]
gammas <- abs(rnorm(nbasis, 0, 1))
y_mean <- as.vector(i_spline_basis_evals %*% as.matrix(gammas))
y = rnorm(Npoints, y_mean, sigma_noise)
spline_svd <- svd(i_spline_basis_evals )
spline_P <- spline_svd$v %*% diag(spline_svd$d) %*% t(spline_svd$v)
spline_U <- spline_svd$u %*% t(spline_svd$v)
stan_data <- list(
Npoints = Npoints,
m = nbasis,
spline_U = spline_U,
spline_P = spline_P,
y = y
)
I have no idea if the thing I do in transformed data is generalizable to other problems. I inspected P and P \gamma to see if any constraint makes sense from Ben’s comment above. But it seems to sample fine without any extra constraints. However, gamma is not quite identified, if you’re really interested in that.
data {
int<lower=1> Npoints;
int<lower=1> m;
matrix[Npoints, m] spline_U;
matrix[m, m] spline_P;
vector[Npoints] y;
}
parameters {
vector[m] theta;
real<lower=0> sigma; // error standard deviation
}
model {
sigma ~ std_normal();
y ~ normal(spline_U * theta, sigma);
}
generated quantities {
vector[m] gamma = mdivide_left_spd(spline_P, reverse(theta));
}
Here’s the summary of gamma from the model
> fit$summary("gamma")
# A tibble: 9 x 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 gamma[1] 1.26 1.26 0.119 0.116 1.07 1.46 1.00 2670. 2339.
2 gamma[2] 0.802 0.807 0.298 0.296 0.306 1.28 1.00 2458. 2470.
3 gamma[3] 0.00609 -0.000957 0.402 0.394 -0.657 0.681 1.00 2444. 2359.
4 gamma[4] 1.07 1.07 0.420 0.417 0.391 1.76 1.00 2450. 2459.
5 gamma[5] 0.0849 0.0886 0.382 0.373 -0.522 0.701 1.00 2354. 2458.
6 gamma[6] 0.347 0.342 0.360 0.339 -0.240 0.939 1.00 2292. 2345.
7 gamma[7] 1.29 1.30 0.353 0.344 0.704 1.86 1.00 2281. 2154.
8 gamma[8] 0.261 0.259 0.297 0.290 -0.203 0.745 1.00 2256. 2396.
9 gamma[9] 1.92 1.92 0.200 0.199 1.60 2.24 1.00 2196. 2561.
vs true gamma
> gammas
[1] 1.37095845 0.56469817 0.36312841 0.63286260 0.40426832 0.10612452 1.51152200 0.09465904 2.01842371
Pairs plot
