Fitting a custom function in rstan

Modeling
,
1

I’m trying to fit a model: Y ~ TestFunc(CBr, aBr, oBr). In the following, I provide the test data and my stan model. The model gives me this error: Stan model does not contain samples, making me believe I formulate the model incorrectly. I’m new to the Bayesian models and very confused, so I appreciate any suggestions to put me on the right track.

stanmodel3 <- "
functions{
  real TestFunc(real yearlingBR,
                    real adultBR, 
                    real oldBR){ 
    
    matrix[13, 13] Lm = rep_matrix(0, 13, 13);
    Lm[1,1] = 0;
    Lm[1, 2] = 0;
    Lm[1, 3] = 0.4 * (0.92 - adultBR);
    Lm[1, 4] = 0.4 * (0.92 - adultBR);
    Lm[1, 5] = 0.4 * (0.95 - adultBR);
    Lm[1, 6] = 0.4 * (0.95 - adultBR);
    Lm[1, 7] = 0.4 * (0.95 - adultBR);
    Lm[1, 8] = 0.4 * (0.9 - adultBR);
    Lm[1, 9] = 0.4 * (0.9 - adultBR);
    Lm[1, 10] = 0.4 * (0.8 - adultBR);
    Lm[1, 11] = 0.4 * (0.67 - oldBR);
    Lm[1, 12] = 0;
    Lm[1, 13] = 0;

    Lm[2,1] = 0.8 - yearlingBR;
    Lm[3,2] = 0.92 - adultBR;
    Lm[4,3] = 0.92 - adultBR;
    Lm[5,4] = 0.95 - adultBR;
    Lm[6,5] = 0.95 - adultBR;
    Lm[7,6] = 0.95 - adultBR;
    Lm[8,7] = 0.9 - adultBR;
    Lm[9,8] = 0.9 - adultBR;
    Lm[10,9] = 0.8 - adultBR;
    Lm[11,10] = 0.67 - oldBR;
    Lm[12,11] = 0.51 - oldBR;
    Lm[13,12] = 0.1 - oldBR;
    vector[13] eigVals = eigenvalues_sym(Lm);
    real gR = max(eigVals);
    return gR ;
  }
}
data {
    int<lower=1> N;   // total number of observations
    vector[N] obsY;      
    vector<lower=0>[N] cBr; 
    vector<lower=0>[N] aBr; 
    vector<lower=0>[N] oBr;
} 
parameters {
    real<lower=0> tau;
} 
transformed parameters {
    real sigma = 1 / sqrt(tau);
    vector[N] mu;
    for(i in 1:N){
        mu[i] = TestFunc(cBr[i], aBr[i], oBr[i]);
    }
} 
model {
    // Priors
    cBr ~ normal(0.4, 0.2);
    aBr ~ normal(0.4, 0.2);
    oBr ~ normal(0.05, 0.025);
    tau ~ gamma(.0001, .0001);
    // Likelihood
    for(i in 1:N){ 
        obsY[i] ~ normal(mu[i], sigma);
    }
       
}
generated quantities{
    vector[N] Y_mean;
    vector[N] Y_pred;
    for(i in 1:N){
        // Posterior parameter distribution of the mean
        Y_mean[i] = TestFunc(cBr[i], aBr[i], oBr[i]);
        // Posterior predictive distribution
        Y_pred[i] = normal_rng(Y_mean[i], sigma);
    }
}
"
testDf <- data.frame(CBR = c(0.2145573, 0.2145573, 0.2145573, 0.2145573, 0.1593126, 0.1593126, 0.1593126, 0.1593126), 
                     ABR = c(0.3983057, 0.3983057, 0.4777045, 0.4777045, 0.3983057, 0.3983057, 0.4777045, 0.4777045), 
                     OBR = c(0.1337283, 0.0581019, 0.1337283, 0.0581019, 0.1337283, 0.0581019, 0.1337283, 0.0581019), 
                     obsY = c(0.6574670, 0.6585795, 0.5765599, 0.5776401, 0.6683592, 0.6694080, 0.5863081, 0.5873235))
dat <- list(
  "N" = nrow(testDf), 
  "obsY" = testDf$obsY,
    "cBr" = testDf$CBR,
    "aBr" = testDf$ABR,
    "oBr" = testDf$OBR)


rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit <- stan(model_code = stanmodel3, 
            model_name = "TestModel", 
            data = dat)
print(fit, c("cBr" "aBr", "oBr"))

.

2

Running that model CmdStan complains:

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: Exception: eigenvalues_sym: m is not symmetric. m[1,2] = 0, but m[2,1] = 0.585443

The model initializes correctly after symmetrizing the matrix

    Lm[13,12] = 0.1 - oldBR;
+   for (i in 1:13) {
+     for (j in 1:i) {
+       Lm[i,j] = Lm[j,i];
+     }
+   }
    vector[13] eigVals = eigenvalues_sym(Lm);
    real gR = max(eigVals);
    return gR ;

Although it looks like the matrix is deliberately non-symmetric. In that case you should use eigenvalues instead of eigenvalues_sym. The eigenvalues of a general matrix can be complex numbers so you need to deal with that (probably by taking the real part only?).

    Lm[13,12] = 0.1 - oldBR;
    real gR = max(get_real(eigenvalues(Lm)));
    return gR ;
1 Like
3

Great! now it works. Thank you very much.