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

saratirvenracs · February 25, 2022, 6:38pm

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

saratirvenracs · July 17, 2022, 8:24pm

I solved the above issue by discretizing the beta distribution differently. Above, I discretized a continuous beta(x,y) using probability densities at 0, 0.01, 0.02, … 0.99, 1. I now discretize it using the cumulative densities at each of these points (for example, I assign to 0.01 the difference of cumulative densities at 0.015 and 0.005). Code below:

 matrix beta_xy(real x, real y){
    matrix[100,3] mat;
    mat = rep_matrix(0,100,3);
    for (i in 1:100){
    mat[i,1] = (i-1)*0.01;
    mat[i,2] = beta_cdf(mat[i,1]+0.005, x, y) - beta_cdf(0.005, x, y);
    mat[i,3] = beta_cdf(mat[i,1]+0.005, y, x) - beta_cdf(0.005, y, x);
    }
    return mat;
  }

Bob_Carpenter · July 18, 2022, 6:17pm

Thanks for providing the solution. You don’t need the rep_matrix line here as you fill in the entire matrix in the loop. Also, it’s not clear why that isn’t a 100 x 2 matrix instead of 100 x 3. Why store the first entires with the step sizes?

saratirvenracs · July 18, 2022, 6:28pm

Indeed, the rep_matrix is redundant (I was translating from R code). Column 1 can be done away with too. In fact, columns 2 and 3 also contain the same sequence in opposite directions. So, column 2 is all that matters. I didn’t have any computation issues when maximizing LogL. I ought to be more efficient with variables when using MCMC methods- thank you.

Topic		Replies	Views
Gradient evaluated at the initial value is not finite although the lp is finite Modeling cognitive-science	12	3778	April 13, 2020
Gradient evaluated at the initial value is not finite Modeling rstan , fitting-issues	1	544	January 1, 2022
Please delete Modeling	4	505	May 14, 2018
'Gradient evaluated at the initial value is not finite' with very simple multinomial model! Modeling fitting-issues	2	387	August 3, 2023
Rejecting initial value - Gradient evaluated at the initial value is not finite Modeling fitting-issues , specification	1	33	August 22, 2024

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

Related topics