Hi, I am using the code below to elicit the posterior distribution of 5 parameters (4 beta distribution parameters and one mean-zero Gaussian error term)- (aBR, bBR), (aGP, bGP), sigma_e (Error) that explain the choices of N=18 individuals in an experiment. There are 8 data points for each individual contained in columns 42 to 49. The beta distributions model beliefs of the individuals about an unknown urn.

To define the likelihood function, I needed the quantile function of beta distribution. Currently, Stan only has a built-in CDF beta_cdf that gives the cumulative density given x, alpha, beta (where 0 <= x <= 1). I tried to combine this built-in function with a binary search to define qbeta- the inverse of the CDF (i.e. it returns 0 <= x <= 1 for a given cumulative density, alpha, beta). I have tested my code repeatedly in R and it works; it is not the most efficient solution but it is only 3 times slower than Râ€™s qbeta function (stats package). However, when I run the code, the NUTS algorithm doesnâ€™t go beyond 1 iteration. It seems to be stuck in an endless loop.

Iâ€™m quite sure Stan should work because I already tried this estimation in JAGS and I got back the parameters I expected (JAGS has the advantage of built-in pbeta and qbeta functions). So, I believe the error is linked to the user-defined function that I am using. I would be grateful if someone can point out my mistake. Iâ€™ve tried debugging with print in several different positions and couldnâ€™t figure out why the code doesnâ€™t work. The 18 rows x 8 columns of data (which form columns 42 to 49 of my original dataset) are available here: ExchData.csv - Google Drive

```
functions{
real qbeta(real p, real alpha, real beta){
real tol = 0.5; // initiate tolerance
real mid = beta_cdf(0.5, alpha, beta); // initiate median
if(p == 1){return(1);} // return 1 for quantile (1)
if(p == 0){return(0);} // return 0 for quantile (0)
// print(p, " ", alpha, " ", beta," ",mid); // for debugging
if(mid == p){return 0.5;} // trivial case
else{
real x = 1; // upper limit of binary search
real y = 0; // lower limit of binary search
while (tol > 0.000001){ // run till required accuracy
if(mid > p){
x = (x+y)/2; // modify upper limit if median is too large
mid = beta_cdf((x+y)/2, alpha, beta); // use in-built beta_cdf to compute new point for binary search
tol = fabs(mid - p); // update deviation for new point
}
if(mid < p){
y = (x+y)/2; // modify lower limit if median is too large
mid = beta_cdf((x+y)/2, alpha, beta); // use in-built beta_cdf to compute new point for binary search
tol = fabs(mid - p); // update deviation for new point
}
}
return((x+y)/2); // return required quantile
}
}
}
// 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=1, upper=50> aBR;
real<lower=1, upper=50> bBR;
real<lower=1, upper=50> aGP;
real<lower=1, upper=50> bGP;
real<lower=0, upper=0.5> sigma_e;
}
// The model to be estimated.
model {
// priors
aBR ~ normal(10, 1);
bBR ~ normal(10, 1);
aGP ~ normal(10, 1);
bGP ~ normal(10, 1);
sigma_e ~ normal(0.1, 0.03);
// likelihood
for (i in 1:N){
Y[ i , 42 ] ~ normal(qbeta(0.5,aBR,bBR), sigma_e);
Y[ i , 43 ] ~ normal(qbeta((beta_cdf(0.6,aBR,bBR)*0.5),aBR,bBR), sigma_e);
Y[ i , 44 ] ~ normal(qbeta((beta_cdf(0.6,aBR,bBR)+(1-beta_cdf(0.6,aBR,bBR))*0.5),aBR,bBR), sigma_e);
Y[ i , 45 ] ~ normal(qbeta((beta_cdf(0.4,aBR,bBR)+(beta_cdf(0.8,aBR,bBR)-beta_cdf(0.4,aBR,bBR))*0.5),aBR,bBR), sigma_e);
Y[ i , 46 ] ~ normal(qbeta(0.5,aGP,bGP), sigma_e);
Y[ i , 47 ] ~ normal(qbeta((beta_cdf(0.8,aGP,bGP)*0.5),aGP,bGP), sigma_e);
Y[ i , 48 ] ~ normal(qbeta((beta_cdf(0.8,aGP,bGP)+(1-beta_cdf(0.8,aGP,bGP))*0.5),aGP,bGP), sigma_e);
Y[ i , 49 ] ~ normal(qbeta((beta_cdf(0.6,aGP,bGP)+(1-beta_cdf(0.6,aGP,bGP))*0.5),aGP,bGP), sigma_e);
}
}
```

**Mod edit:** For syntax highlighting. In the future use ``` to start and end the code block. To get Stan highlighting add stan as in

```
```stan
```