The upcoming SBC event led to me trying to implement an SBC approach to a non-linear model I’m working on.
The predictors are matrices of observed predictors over time relative to an observed outcome.
As the model is non-linear, the parameters are dependent on observed pedictor data, so I have coded the number of simulations as the number of observations. I’m curious as to whether or not SBC can/could/should be implemented to calibrating this kind of model.
Here is what I have coded as an initial attempt:
data {
int<lower=1> N; // total number of observations
vector[N] y; // response variable
int<lower=1> N_Times; // Number of times in time series
matrix[N, N_Times] X1;
matrix[N, N_Times] X2;
vector[N] X3;
}
transformed data {
real<lower=0> a_sim = normal_rng(5,5);
real<lower=0> b_sim = normal_rng(5,5);
real<lower=0> c_sim = normal_rng(5,5);
real<lower=0> kappa_sim = exponential_rng(1);
vector[N] y_pred;
vector<lower=0>[N] mu_sim;
for (n in 1:N) {
mu_sim[n] = sum(row(X1, n) .* (a_sim * exp(b_sim * row(X2, n)) + c_sim)) / X3[n];
y_pred[n] = beta_rng(inv_logit(mu_sim[n]) * kappa_sim, (1 - inv_logit(mu_sim[n)) * kappa);
}
}
parameters {
real<lower=0> a;
real<lower=0> b;
real<lower=0> c;
real<lower=0> kappa;
}
model {
vector[N] mu;
// priors including constants
a ~ normal(5,5);
b ~ normal(5,5);
c ~ normal(5,5);
kappa ~ exponential(1);
// initialize mu for beta likelihood
for (n in 1:N) {
mu[n] = sum(row(X1, n) .* (a * exp(b * row(X2, n)) + c)) / X3[n];
}
y_pred ~ beta(inv_logit(mu) * kappa, (1 - inv_logit(mu)) * kappa);
}
generated quantities {
int<lower = 0, upper = 1> lt_sim[3]
= { a < a_sim, b < b_sim, c < c_sim };
}
The code compiles, runs, and returns results and I will still need to tackle thinning and the additional steps, but I’m wondering if it’s a valid approach.
Any feedback would be greatly appreciated!