# Fitting a custom function in rstan

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 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[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"))
``````

.

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

Great! now it works. Thank you very much.