I’m getting brain freeze and I can’t figure out how to rollup the SN parameters through the density functions for likelihood and prior: here’s where I am - hopefully you can see what I’m trying to do and where I died (not very far, sorry). Thanks for your help:

stanMod <- ’
data {
int N;
vector[N] y;
real prXi;
real prOmega;
real prAlpha;
}
parameters {
real xi;
real omega;
real alpha;
//real postXi; not sure what to do here - I tried a few things and it didn’t work
//real postOmega;
//real postAlpha;

}
model {

y ~ skew_normal(xi, omega, alpha);

// really not sure what to do here - how to roll up the parameters through the distributions
// so I’m not doing anything with the prior parameters passed in yet

}
’
N = 1000
y <- rsn( N , xi=1 , omega=3 , alpha=10 )
mod <- stan_model(model_code = stanMod)
fit <- sampling(mod, data = list(y= y, prXi = 3 , prOmega = 5, prAlpha = 10, N = 1000))
print(fit)

I don’t know what “rollup” means. You do need to restrict omega to be positive by declaring it as

parameters {
real xi;
real<lower=0> omega;
real alpha;
}

I am not sure what you were intending with real postXi etc. The parameters are xi, omega, and alpha and Stan draws from their (posterior) distribution conditional on y. Similarly, having real prXi; etc. in the data block is insufficient. You need a prior distribution that describes what you believe xi, omega, and alpha to be before seeing the realizations of y. Then, in the model block you call the (log)density function for that prior distribution given whatever hyperparameters you want to specify, such as

target += normal_lpdf(xi | 0, 5);

Or better, you can pass those hyperparameters in as data

data {
...
real xi_loc;
real<lower=0> xi_scale;
}
parameters {
real xi;
real<lower=0> omega;
real alpha;
}
model {
target += skew_normal_lpdf(y | xi, omega, alpha);
target += normal_lpdf(xi | xi_loc, xi_scale);
// and similarly for omega and alpha
}