Yet another case of (Rejecting initial value: Gradient evaluated at the initial value is not finite.)

Hi, I have a sophisticated Decision Theory model to estimate based on an experiment. The model treats subjects’ responses as a function of parameters and Gaussian errors (sigma_r, sigma_a, sigma_e- one error for each task type in the experiment). The model values are computed using a numerical integration (to give some context: consider an urn with 100 blue and red balls, whose composition is unknown. The model assumes that subjects considers all compositions possible: 1 blue + 99 red, 2 blue + 98 red, …, 98 blue + 2 red, 99 blue + 1 red. Such “beliefs” are modelled as beta distributions and there’s a summation over 99 possible beliefs using a for loop from 1 to 99). I have used print statements in different locations to check that the model values are being computed as intended. Everything seems fine but there’s no sampling that takes place. I’ve looked at different threads and tried many ideas- making my priors stronger, constraining parameters to narrower ranges, using lpdf instead of densities, initializing close to expected values, etc. Despite all this, there’s no resolution to the error. The code and data are attached below. I use inc_beta_inverse suggested by @spinkney (removed in code below for brevity) and a few functions that I defined to carry out the numerical integration. In a few places, there’s underflow/ overflow due to beta distribution which I’ve rectified using is_nan. I also print the parameter initial values, model values for this set of initial parameters and the log-likelihoods below and surprisingly, they’re all finite and completely reasonable. So I don’t understand why the gradient is not finite. Is the model too complicated for gradient computation? I don’t want to switch back to JAGS (I just migrated to Stan and I highly appreciate the capability to build our own functions- in JAGS, I had to copy-paste my code 99 times instead of using a for loop!) but I could at least run the chain in JAGS.


functions{
  real wpkt(real a, real p){return (p^a)/((p^a+(1-p)^a)^(1/a));}
  real phi(real theta, real rho, real x){return((1-exp(-theta*(x^rho)))/(1-exp(-theta)));}
  real phi_inv(real theta, real rho, real y){return(((-log(1-(y*(1-exp(-theta)))))/theta)^(1/rho));}
  matrix beta_xy(real x, real y){
    matrix[101,5] mat;
    for (i in 1:101){
    mat[i,1] = (i-1)*0.01;
    mat[i,2] = exp(beta_lpdf(mat[i,1] | x, y))*0.01;
    mat[i,3] = exp(beta_lpdf(mat[i,1] | y, x))*0.01;
    if(i==1)mat[i,4] = mat[i,2];
    if(i==1)mat[i,5] = mat[i,3];
    if(i>1) mat[i,4] = mat[i-1,4] + mat[i,2];
    if(i>1) mat[i,5] = mat[i-1,5] + mat[i,3];
    }
    return mat;
  }
  vector model_ambi(vector params){
    vector[32] ambiguous;
    vector[32] sum_amb;
    real u03agg = (0^params[1])/(200^params[1]);
    real uf3agg = (200^params[1])/(200^params[1]);
    real u13agg = (0^params[1])/(200^params[1]);
    real ue3agg = (150^params[1])/(200^params[1]);
    real u23agg = (0^params[1])/(200^params[1]);
    real ud3agg = (100^params[1])/(200^params[1]);
    real u33agg = (0^params[1])/(200^params[1]);
    real uc3agg = (50^params[1])/(200^params[1]);
    real u43agg = (50^params[1])/(200^params[1]);
    real ub3agg = (200^params[1])/(200^params[1]);
    real u53agg = (50^params[1])/(200^params[1]);
    real ua3agg = (150^params[1])/(200^params[1]);
    real u63agg = (25^params[1])/(200^params[1]);
    real u93agg = (125^params[1])/(200^params[1]);
    real u73agg = (25^params[1])/(200^params[1]);
    real u83agg = (75^params[1])/(200^params[1]);
    matrix[101,5] BRagg;
    matrix[101,5] GPagg;
    BRagg = beta_xy(params[6],params[7]);
    GPagg = beta_xy(params[8],params[9]);
    sum_amb = rep_vector(0,32);
   
    for (k in 1:99) {
    sum_amb[01] = sum_amb[01] + phi(params[3],params[4],uf3agg*wpkt(params[2],(100-k)*0.01)+u03agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[02] = sum_amb[02] + phi(params[3],params[4],ue3agg*wpkt(params[2],(100-k)*0.01)+u13agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[03] = sum_amb[03] + phi(params[3],params[4],ud3agg*wpkt(params[2],(100-k)*0.01)+u23agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[04] = sum_amb[04] + phi(params[3],params[4],uc3agg*wpkt(params[2],(100-k)*0.01)+u33agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[05] = sum_amb[05] + phi(params[3],params[4],ub3agg*wpkt(params[2],(100-k)*0.01)+u43agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[06] = sum_amb[06] + phi(params[3],params[4],ua3agg*wpkt(params[2],(100-k)*0.01)+u53agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[07] = sum_amb[07] + phi(params[3],params[4],u93agg*wpkt(params[2],(100-k)*0.01)+u63agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[08] = sum_amb[08] + phi(params[3],params[4],u83agg*wpkt(params[2],(100-k)*0.01)+u73agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])))?(0):(wpkt(params[5],BRagg[k+1,5])-wpkt(params[5],BRagg[k,5])));  
    sum_amb[09] = sum_amb[09] + phi(params[3],params[4],uf3agg*wpkt(params[2],(100-k)*0.01)+u03agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[10] = sum_amb[10] + phi(params[3],params[4],ue3agg*wpkt(params[2],(100-k)*0.01)+u13agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[11] = sum_amb[11] + phi(params[3],params[4],ud3agg*wpkt(params[2],(100-k)*0.01)+u23agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[12] = sum_amb[12] + phi(params[3],params[4],uc3agg*wpkt(params[2],(100-k)*0.01)+u33agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[13] = sum_amb[13] + phi(params[3],params[4],ub3agg*wpkt(params[2],(100-k)*0.01)+u43agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[14] = sum_amb[14] + phi(params[3],params[4],ua3agg*wpkt(params[2],(100-k)*0.01)+u53agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[15] = sum_amb[15] + phi(params[3],params[4],u93agg*wpkt(params[2],(100-k)*0.01)+u63agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[16] = sum_amb[16] + phi(params[3],params[4],u83agg*wpkt(params[2],(100-k)*0.01)+u73agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])))?(0):(wpkt(params[5],BRagg[k+1,4])-wpkt(params[5],BRagg[k,4])));  
    sum_amb[17] = sum_amb[17] + phi(params[3],params[4],uf3agg*wpkt(params[2],(100-k)*0.01)+u03agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[18] = sum_amb[18] + phi(params[3],params[4],ue3agg*wpkt(params[2],(100-k)*0.01)+u13agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[19] = sum_amb[19] + phi(params[3],params[4],ud3agg*wpkt(params[2],(100-k)*0.01)+u23agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[20] = sum_amb[20] + phi(params[3],params[4],uc3agg*wpkt(params[2],(100-k)*0.01)+u33agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[21] = sum_amb[21] + phi(params[3],params[4],ub3agg*wpkt(params[2],(100-k)*0.01)+u43agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[22] = sum_amb[22] + phi(params[3],params[4],ua3agg*wpkt(params[2],(100-k)*0.01)+u53agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[23] = sum_amb[23] + phi(params[3],params[4],u93agg*wpkt(params[2],(100-k)*0.01)+u63agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[24] = sum_amb[24] + phi(params[3],params[4],u83agg*wpkt(params[2],(100-k)*0.01)+u73agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])))?(0):(wpkt(params[5],GPagg[k+1,5])-wpkt(params[5],GPagg[k,5])));  
    sum_amb[25] = sum_amb[25] + phi(params[3],params[4],uf3agg*wpkt(params[2],(100-k)*0.01)+u03agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[26] = sum_amb[26] + phi(params[3],params[4],ue3agg*wpkt(params[2],(100-k)*0.01)+u13agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[27] = sum_amb[27] + phi(params[3],params[4],ud3agg*wpkt(params[2],(100-k)*0.01)+u23agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[28] = sum_amb[28] + phi(params[3],params[4],uc3agg*wpkt(params[2],(100-k)*0.01)+u33agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[29] = sum_amb[29] + phi(params[3],params[4],ub3agg*wpkt(params[2],(100-k)*0.01)+u43agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[30] = sum_amb[30] + phi(params[3],params[4],ua3agg*wpkt(params[2],(100-k)*0.01)+u53agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[31] = sum_amb[31] + phi(params[3],params[4],u93agg*wpkt(params[2],(100-k)*0.01)+u63agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    sum_amb[32] = sum_amb[32] + phi(params[3],params[4],u83agg*wpkt(params[2],(100-k)*0.01)+u73agg*(1-wpkt(params[2],(100-k)*0.01)))*(is_nan((wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])))?(0):(wpkt(params[5],GPagg[k+1,4])-wpkt(params[5],GPagg[k,4])));  
    }
    ambiguous[01] = 200*(phi_inv(params[3],params[4],sum_amb[01])^(1/params[1]));
    ambiguous[02] = 200*(phi_inv(params[3],params[4],sum_amb[02])^(1/params[1]));
    ambiguous[03] = 200*(phi_inv(params[3],params[4],sum_amb[03])^(1/params[1]));
    ambiguous[04] = 200*(phi_inv(params[3],params[4],sum_amb[04])^(1/params[1]));
    ambiguous[05] = 200*(phi_inv(params[3],params[4],sum_amb[05])^(1/params[1]));
    ambiguous[06] = 200*(phi_inv(params[3],params[4],sum_amb[06])^(1/params[1]));
    ambiguous[07] = 200*(phi_inv(params[3],params[4],sum_amb[07])^(1/params[1]));
    ambiguous[08] = 200*(phi_inv(params[3],params[4],sum_amb[08])^(1/params[1]));
    ambiguous[09] = 200*(phi_inv(params[3],params[4],sum_amb[09])^(1/params[1]));
    ambiguous[10] = 200*(phi_inv(params[3],params[4],sum_amb[10])^(1/params[1]));
    ambiguous[11] = 200*(phi_inv(params[3],params[4],sum_amb[11])^(1/params[1]));
    ambiguous[12] = 200*(phi_inv(params[3],params[4],sum_amb[12])^(1/params[1]));
    ambiguous[13] = 200*(phi_inv(params[3],params[4],sum_amb[13])^(1/params[1]));
    ambiguous[14] = 200*(phi_inv(params[3],params[4],sum_amb[14])^(1/params[1]));
    ambiguous[15] = 200*(phi_inv(params[3],params[4],sum_amb[15])^(1/params[1]));
    ambiguous[16] = 200*(phi_inv(params[3],params[4],sum_amb[16])^(1/params[1]));
    ambiguous[17] = 200*(phi_inv(params[3],params[4],sum_amb[17])^(1/params[1]));
    ambiguous[18] = 200*(phi_inv(params[3],params[4],sum_amb[18])^(1/params[1]));
    ambiguous[19] = 200*(phi_inv(params[3],params[4],sum_amb[19])^(1/params[1]));
    ambiguous[20] = 200*(phi_inv(params[3],params[4],sum_amb[20])^(1/params[1]));
    ambiguous[21] = 200*(phi_inv(params[3],params[4],sum_amb[21])^(1/params[1]));
    ambiguous[22] = 200*(phi_inv(params[3],params[4],sum_amb[22])^(1/params[1]));
    ambiguous[23] = 200*(phi_inv(params[3],params[4],sum_amb[23])^(1/params[1]));
    ambiguous[24] = 200*(phi_inv(params[3],params[4],sum_amb[24])^(1/params[1]));
    ambiguous[25] = 200*(phi_inv(params[3],params[4],sum_amb[25])^(1/params[1]));
    ambiguous[26] = 200*(phi_inv(params[3],params[4],sum_amb[26])^(1/params[1]));
    ambiguous[27] = 200*(phi_inv(params[3],params[4],sum_amb[27])^(1/params[1]));
    ambiguous[28] = 200*(phi_inv(params[3],params[4],sum_amb[28])^(1/params[1]));
    ambiguous[29] = 200*(phi_inv(params[3],params[4],sum_amb[29])^(1/params[1]));
    ambiguous[30] = 200*(phi_inv(params[3],params[4],sum_amb[30])^(1/params[1]));
    ambiguous[31] = 200*(phi_inv(params[3],params[4],sum_amb[31])^(1/params[1]));
    ambiguous[32] = 200*(phi_inv(params[3],params[4],sum_amb[32])^(1/params[1]));
    return ambiguous;
  }
}

// The input data is a vector 'y' of length 'N'.
data {
  int<lower=0> N;
  int<lower=0> M;
  matrix[N, M] Y;
}

// The parameters accepted by the model.
parameters {
  real<lower=0.5, upper=1> u;
  real<lower=0.3, upper=1> wr;
  real<lower=5, upper=30> sigma_r;
  real<lower=0, upper=10> theta;
  real<lower=0, upper=3> rho;
  real<lower=0.3, upper=1> wa;
  real<lower=5, upper=30> sigma_a;
  real<lower=1, upper=10> aBR;
  real<lower=1, upper=10> bBR;
  real<lower=1, upper=10> aGP;
  real<lower=1, upper=10> bGP;
  real<lower=0.01, upper=0.25> sigma_e;
}

// The model to be estimated.
model {
  int counter_risky = 0;
  int counter_ambi = 0;
  int counter_exch = 0;
  
  vector[9] cache_risky = [((200^u)*(((0.1)^wr)/(((0.1)^wr+(1-(0.1))^wr)^(1/wr))) + (0^u)*(1-(((0.1)^wr)/(((0.1)^wr+(1-(0.1))^wr)^(1/wr)))))^(1/u),
                           ((200^u)*(((0.3)^wr)/(((0.3)^wr+(1-(0.3))^wr)^(1/wr))) + (0^u)*(1-(((0.3)^wr)/(((0.3)^wr+(1-(0.3))^wr)^(1/wr)))))^(1/u),
                           ((200^u)*(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr))) + (0^u)*(1-(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr)))))^(1/u),
                           ((200^u)*(((0.7)^wr)/(((0.7)^wr+(1-(0.7))^wr)^(1/wr))) + (0^u)*(1-(((0.7)^wr)/(((0.7)^wr+(1-(0.7))^wr)^(1/wr)))))^(1/u),
                           ((200^u)*(((0.9)^wr)/(((0.9)^wr+(1-(0.9))^wr)^(1/wr))) + (0^u)*(1-(((0.9)^wr)/(((0.9)^wr+(1-(0.9))^wr)^(1/wr)))))^(1/u),
                           ((160^u)*(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr))) + (20^u)*(1-(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr)))))^(1/u),
                           ((120^u)*(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr))) + (40^u)*(1-(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr)))))^(1/u),
                           ((100^u)*(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr))) + (0^u)*(1-(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr)))))^(1/u),
                           ((60^u)*(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr))) + (0^u)*(1-(((0.5)^wr)/(((0.5)^wr+(1-(0.5))^wr)^(1/wr)))))^(1/u)]';
  vector[32] cache_ambi = model_ambi([u,wr,theta,rho,wa,aBR,bBR,aGP,bGP]');
  vector[8] cache_exch = [inc_beta_inverse(0.5,aBR,bBR),
                          inc_beta_inverse((beta_cdf(0.6,aBR,bBR)*0.5),aBR,bBR),
                          inc_beta_inverse((beta_cdf(0.6,aBR,bBR)+(1-beta_cdf(0.6,aBR,bBR))*0.5),aBR,bBR),
                          inc_beta_inverse((beta_cdf(0.4,aBR,bBR)+(beta_cdf(0.8,aBR,bBR)-beta_cdf(0.4,aBR,bBR))*0.5),aBR,bBR),
                          inc_beta_inverse(0.5,aGP,bGP),
                          inc_beta_inverse((beta_cdf(0.8,aGP,bGP)*0.5),aGP,bGP),
                          inc_beta_inverse((beta_cdf(0.8,aGP,bGP)+(1-beta_cdf(0.8,aGP,bGP))*0.5),aGP,bGP),
                          inc_beta_inverse((beta_cdf(0.6,aGP,bGP)+(1-beta_cdf(0.6,aGP,bGP))*0.5),aGP,bGP)]';
  
  print(u," ",wr," ",sigma_r," ",theta," ",rho," ",wa," ",sigma_a," ",aBR," ",bBR," ",aGP," ",bGP," ",sigma_e);
  print(cache_risky);
  print(cache_ambi);
  print(cache_exch);
                   
  // priors
  u ~ normal(0.7, 0.1);
  wr ~ normal(0.6, 0.1);
  sigma_r ~ normal(20, 2);
  theta ~ normal(5, 0.5);
  rho ~ normal(2, 0.1);
  wa ~ normal(0.6, 0.1);
  sigma_a ~ normal(20, 2);
  aBR ~ normal(6, 1);
  bBR ~ normal(6, 1);
  aGP ~ normal(6, 1);
  bGP ~ normal(6, 1);
  sigma_e ~ normal(0.1, 0.02);
  
  // likelihood
    for (ir in 1:9) {
      counter_risky += 1;
     // Y[ : , ir ] ~ normal(cache_risky[counter_risky], sigma_r);
      target += normal_lpdf(Y[ : , ir ] | cache_risky[counter_risky], sigma_r);
      print(normal_lpdf(Y[ : , ir ]|cache_risky[counter_risky], sigma_r));
      }
      
    for (ia in 10:41) {
      counter_ambi += 1;
      // Y[ : , ia ] ~ normal(cache_ambi[counter_ambi], sigma_a);
      target += normal_lpdf(Y[ : , ia ] | cache_ambi[counter_ambi], sigma_a);
      print(normal_lpdf(Y[ : , ia ]|cache_ambi[counter_ambi], sigma_a));
      }
    
    for (ie in 42:49) {
      counter_exch += 1;
      // Y[ : , ie ] ~ normal(cache_exch[counter_exch], sigma_e);
      target += normal_lpdf(Y[ : , ie ] | cache_exch[counter_exch], sigma_e);
      print(normal_lpdf(Y[ : , ie ]|cache_exch[counter_exch], sigma_e));
      }
}


Data (18 subjects, 49 responses each, columns 1 to 9 of type r (risky), 10 to 41 of type a (ambiguity), 42 to 49 of type e (exchangeability): Data.csv - Google Drive

Print output:
Chain 1: 0.7 0.6 23 5 2 0.5 20 5 3 8 3 0.08 \ parameters
[18.3686,38.6534,57.057,79.909,120.781,68.4516,70.1088,28.5285,17.1171] \ model values
[57.5764,45.7787,32.387,17.104,103.566,88.3043,62.5247,45.2922,37.4191,29.605,20.7837,10.8787,86.1647,75.9938,50.3539,38.8037,70.0211,55.1646,38.6624,20.2391,113.675,95.0377,69.1429,48.5912,30.4169,23.7486,16.4382,8.48862,80.4828,71.759,46.0751,36.4311] \ model values
[0.635884,0.48966,0.728204,0.624452,0.741425,0.683056,0.857345,0.769137] \ model values
\ log-likelihoods below (not sure why the last 8 are positive)
-75.7156
-76.4277
-83.6016
-91.9167
-101.897
-81.9225
-76.6586
-75.0969
-74.7254
-90.6924
-80.6223
-76.5595
-72.494
-87.839
-78.1986
-76.8485
-73.2168
-82.2842
-76.3244
-73.3673
-73.485
-75.3047
-72.6258
-74.7296
-71.7283
-120.785
-94.0934
-82.0316
-75.0529
-103.39
-84.5975
-81.5045
-73.6528
-83.5774
-74.2217
-74.7448
-71.2936
-79.4868
-74.0253
-73.713
-71.9951
17.0846
16.8288
22.3359
24.3734
9.93335
18.067
27.1445
23.8726

1 Like