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