How to define a array of parameter values

Interfaces brms
,
1

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:
data

2

Hi,
sorry for not getting to you earlier. I don’t think you can easily define an array of variables in brms. There might be some options if you decide to completely hijack what brms does (as discussed e.g. at Brms & pk models? - #16 by martinmodrak). But I think you are running into limits of what brms can do for you and it might make more sense to just move to plain Stan.

Best of luck with your model!