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