Hello,
even though I know the spline-function of the mgcv-package can be used for brms, I am trying to customize my fitting function, which is why I am defining my own function. The code-example should only exemplify my current problem.
Using the brms-interface, would it be possible to define the a parameters as an k x 1 - array where I could either define the spline-parameters as random-effects (with a normal population prior with a common mu- and sigma-parameter) or giving every individual array-element its own prior?
The code for reproducibility purposes:
f_test_stan =
"
// B-Spline basis function
real f_B(int i, // index
int d, // degree
real x) { // variable
real a;
real b;
real out;
// initialize values
row_vector[15] knots = [ -4, -2, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24 ];
// knots(n = 10, d = 2, min = 0, max = 20)
if (d == 0) {
if (x >= knots[i] && x < knots[i + 1]) {
out = 1;
} else {
out = 0;
}
} else {
a = (x - knots[i]) / (knots[i + d] - knots[i]);
b = (knots[i + d + 1] - x) / (knots[i + d + 1] - knots[i + 1]);
out = a * f_B(i, d - 1, x) + b * f_B(i + 1, d - 1, x);
}
return out;
}
// Spline function
real f_spline(real x, // variable
real a1, // parameters
real a2,
real a3,
real a4,
real a5,
real a6,
real a7,
real a8,
real a9,
real a10,
real a11,
real a12) {
real out;
int d;
// initialize values
d = 2;
// There are n + d many B-splines
out = a1 * f_B(1, d, x);
out = out + a2 * f_B(2, d, x);
out = out + a3 * f_B(3, d, x);
out = out + a4 * f_B(4, d, x);
out = out + a5 * f_B(5, d, x);
out = out + a6 * f_B(6, d, x);
out = out + a7 * f_B(7, d, x);
out = out + a8 * f_B(8, d, x);
out = out + a9 * f_B(9, d, x);
out = out + a10 * f_B(10, d, x);
out = out + a11 * f_B(11, d, x);
out = out + a12 * f_B(12, d, x);
return out;
}
"
test_data = data.table(x = seq(0, 20, 0.01))
test_data[,y := sin(x) + rnorm(1, 0, 0.1), by = 'x']
prior_test = c(set_prior("normal(0, 1)", nlpar = "a1"),
set_prior("normal(0, 1)", nlpar = "a2"),
set_prior("normal(0, 1)", nlpar = "a3"),
set_prior("normal(0, 1)", nlpar = "a4"),
set_prior("normal(0, 1)", nlpar = "a5"),
set_prior("normal(0, 1)", nlpar = "a6"),
set_prior("normal(0, 1)", nlpar = "a7"),
set_prior("normal(0, 1)", nlpar = "a8"),
set_prior("normal(0, 1)", nlpar = "a9"),
set_prior("normal(0, 1)", nlpar = "a10"),
set_prior("normal(0, 1)", nlpar = "a11"),
set_prior("normal(0, 1)", nlpar = "a12"))
eq_test= bf(y ~ f_spline(x, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12),
a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10 + a11 + a12 ~ 1,
nl = TRUE)
stan_vars_test = stanvar(scode = f_test_stan, block = "functions")
mod_test = brm(eq_test,
family = gaussian(),
data = test_data,
prior = prior_test,
stanvars = stan_vars_test,
cores = 4,
iter = 1000,
control = list(adapt_delta = 0.5,
max_treedepth = 15))
- Operating System: aarch64-apple-darwin20 (darwin20)
- brms Version: 2.15.9, R-Version (4.1.0)
Plot of data:
