Prior predictive check doesn't make sense for exp_mod_normal_rng

Helo rstan experts,
I am trying to fit my response time data to an exGaussian model.
As I kept getting the error without being able to sample,

exp_mod_normal_lpdf: Inv_scale parameter is -2.36578, but must be > 0!
exp_mod_normal_lpdf: Scale parameter is -2.08391, but must be > 0!

I decided to check the prior first.

(I have 2 factors (a 2 by 3 design) with 5 contrasts, and then I set 5 contrasts respectively for mu, sigma, and lambda. Therefore I have 15 betas in total but I just center all betas on 0)

data {
int<lower = 1> N_trials;
int<lower = -1,upper = 1> cond[N_trials];
int<lower = -1,upper = 1> type1[N_trials];
int<lower = -1,upper = 1> type2[N_trials];
}

generated quantities {
real y_sim[N_trials];
real mu[N_trials];
real sigma[N_trials];
real lambda[N_trials];
real alpha1 = normal_rng (210,50);
real alpha2 = normal_rng (40,5);
real alpha3 = normal_rng (0.005,0.0005);
real beta1 = normal_rng (0, 10);
real beta2 = normal_rng (0, 1);
real beta3 = normal_rng (0, 0.00025);
real beta4 = normal_rng (0, 10);
real beta5 = normal_rng (0, 1);
real beta6 = normal_rng (0, 0.000125);
real beta7 = normal_rng (0, 10);
real beta8 = normal_rng (0, 1);
real beta9 = normal_rng (0, 0.000125);
real beta10 = normal_rng (0, 10);
real beta11 = normal_rng (0, 1);
real beta12 = normal_rng (0, 0.0000625);
real beta13 = normal_rng (0, 10);
real beta14 = normal_rng (0, 1);
real beta15 = normal_rng (0, 0.0000625);
for(i in 1:N_trials){
mu[i] = alpha1 + beta1 * cond[i] + beta4 * type1[i] + beta7 * type2[i] + beta10 * cond[i] * type1[i] + beta13 * cond[i] * type2[i];
sigma[i] = alpha2 + beta2 * cond[i] + beta5 * type1[i] + beta8 * type2[i] + beta11 * cond[i] * type1[i] + beta14 * cond[i] * type2[i];
lambda[i] = alpha3 + beta3 * cond[i] + beta6 * type1[i] + beta9 * type2[i] + beta12 * cond[i] * type1[i] + beta15 * cond[i] * type2[i];
y_sim[i] = exp_mod_normal_rng (mu[i] , sigma[i] , lambda[i]);}
}

This prior predictive check gives me results of mu, sigma, lambda being what I want (mu around 210, sigma around 40, and lambda around 1/200). However, my y_sim’s distribution is so restricted (373.0940, 450.9805), which is not exGaussian at all.

If I simply do on R,
range(rnorm(50000,210,40)+rexp(50000,1/200)),
I would get something like
[1] 60.69044 2546.52426 which has a long tail

I wonder what I have done wrong. At first I thought this is because the mu in rstan might stand for grand mean of the distribution so I also tried real alpha1 = normal_rng (410,50), but it still gave me an obviously non-exGaussian distribution.

I also attach here my stan script for the real sampling (where I also got the “scale parameter must >0” error), if that helps with solving my ultimate issue preventing me from fitting the RTs.
fitting_exGaussian.stan (1.6 KB)
df_forstan.csv (898.4 KB)

df ← read.csv(“df_forstan.csv”)
list_forstan ← list(N_trials=nrow(df),
cond=df$cond,
type1=df$type1,
type2=df$type2,
rt=df$rt,
I=length(unique(df$subject)),
J=length(unique(df$item)),
subject=df$subject,
item=df$item)
full_model ← stan_model(file = “fitting_exGaussian.stan”)

Thanks a lot.

best,
Kuan-Jung