copy of model_1_1_choice.stan

model_1_1_choice.stan

functions {
  real f(real mu, real s, real x){
      // separately performing numerical integration to make this computing fast
      return(normal_lpdf(x | mu, s));
      }
     
  real log_lik_Simpson(real mu, real s, real a, real b, int M) {
    // numerical integration using the Simpson's rule
    vector[M+1] lp;
    real h;
    h = (b-a)/M;
    lp[1] = f(mu, s, a);
    for (m in 1:(M/2))
      lp[2*m] = log(4) + f(mu, s, a + h*(2*m-1));
    for (m in 1:(M/2-1))
      lp[2*m+1] = log(2) + f(mu, s,  a + h*2*m);
    lp[M+1] = f(mu, s, b);
    return(log(h/3) + log_sum_exp(lp));
  }
  
  real f2(int Y, real cwl, real cwno, real ll, real lno, real C_width, real Leg_lack,
          real bias_no0,
          real bias_l_K){
    vector[3] mu;
    mu[1] = 0;
    mu[2] = bias_l_K + cwl*C_width + ll*Leg_lack;
    mu[3] = bias_no0 + cwno*C_width + lno*Leg_lack; 
    return(categorical_logit_lpmf(Y | mu));
  }
  
  
  real log_lik_Simpson_rest(int Y, real cwl, real cwno, real ll, real lno, real C_width, real Leg_lack,
                            real bias_no0,
                            real bias_l_lower, real bias_l_upper,
                            int M){
    vector[M+1] lp;
    real h;
    h = (bias_l_upper - bias_l_lower)/M;
    lp[1] = f2(Y, cwl, cwno, ll, lno, C_width, Leg_lack,
                bias_no0,
                bias_l_lower);
    for (m in 1:(M/2))
      lp[2*m] = log(4) + f2(Y, cwl, cwno, ll, lno, C_width, Leg_lack,
                            bias_no0,
                            bias_l_lower + h*(2*m-1));
    for (m in 1:(M/2-1))
      lp[2*m+1] = log(2) + f2(Y, cwl, cwno, ll, lno, C_width, Leg_lack,
                              bias_no0,
                              bias_l_lower + h*2*m);
    lp[M+1] = f2(Y, cwl, cwno, ll, lno, C_width, Leg_lack,
                  bias_no0,
                  bias_l_upper);
    return(log(h/3) + log_sum_exp(lp));
  }
}

data {
    int N; // number of actions
    int K; // number of choices
    int L; // number of individuals
    int<lower=1, upper=K> Y[N]; // Here, we have three choices
    real C_width[N]; // carapace width
    int Leg_lack[N]; // degree of leg lack
    int<lower=1, upper=L> ID[N]; // ID of the crabs
}

parameters {
  // Parameters are declared in lowercase letters
    real bias_l[L]; // local parameters incorpolated to express 'individuality' of the choices
    real bias_l0; // mean of bias_l
    real ll; 
    real cwl;
    real<lower=0> s_bias_l; // sd of bias_l
    
    real cwno;
    real bias_no0;
    real lno;
}

transformed parameters {
    matrix[N,K] mu;
    
    for (n in 1:N){
        mu[n,1] = 0; // M size choice
        mu[n,2] =  bias_l[ID[n]] + cwl*C_width[n] + ll*Leg_lack[n]; // L size choice
        mu[n,3] =  bias_no0 + cwno*C_width[n] + lno*Leg_lack[n]; // no choice (abondon the choice itself
    }
}

model {
    for (l in 1:L){
        bias_l[l] ~ normal(bias_l0, s_bias_l);
    }
    for (n in 1:N){
        Y[n] ~ categorical_logit(mu[n,]');
    }
}

generated quantities {
  vector[N] log_lik;
  real log_lik_l;
  // for making computing fast, separately compute loglik unrelated to Y[n] and other parameters
  log_lik_l = log_lik_Simpson(bias_l0, s_bias_l, 
                              bias_l0 - 5*s_bias_l, bias_l0 + 5*s_bias_l, 100); // numerical integration from -5*sigma to +5*sigma
  for (n in 1:N){
    log_lik[n] = log_lik_l + 
                 log_lik_Simpson_rest(Y[n], 
                                      cwl, cwno, ll, lno, C_width[n], Leg_lack[n],
                                      bias_no0,
                                      bias_l0 - 5*s_bias_l, bias_l0 + 5*s_bias_l,
                                      100); 
  }
}