Custom beta quantile function stuck in iteration 1

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

You can try this version helpful_stan_functions/inc_beta_inverse.stanfunctions at main · spinkney/helpful_stan_functions · GitHub

1 Like

Thank you for pointing out- I didn’t know the algorithm was implemented for Stan! Perhaps they can add it to the next manual.

I could get my estimation to work using inc_beta_inverse. It is somehow still slower than JAGS (which uses R’s qbeta function- also based on the same algorithm) but I attained convergence within 30 mins or so and that’s perfectly fine!

Great! This is not in Stan. That repo is a personal repo that I put things that are helpful to me. Anyone is free to add Stan code they have found helpful. The incomplete beta inverse should be in the next version of Stan the PR is Incomplete Beta Function Inverse by andrjohns · Pull Request #2637 · stan-dev/math · GitHub. This will be faster b/c it includes the derivatives.

To eek out more speed, you can remove the loop over the observations and cache the calculation of the mean values of the normal.

...
 // likelihood
vector[8] cache = [qbeta(0.5,aBR,bBR),
                   qbeta((beta_cdf(0.6,aBR,bBR)*0.5),aBR,bBR),
                   qbeta((beta_cdf(0.6,aBR,bBR)+(1-beta_cdf(0.6,aBR,bBR))*0.5),aBR,bBR),
                   qbeta((beta_cdf(0.4,aBR,bBR)+(beta_cdf(0.8,aBR,bBR)-beta_cdf(0.4,aBR,bBR))*0.5),aBR,bBR),
                   qbeta(0.5,aGP,bGP),
                   qbeta((beta_cdf(0.8,aGP,bGP)*0.5),aGP,bGP),
                   qbeta((beta_cdf(0.8,aGP,bGP)+(1-beta_cdf(0.8,aGP,bGP))*0.5),aGP,bGP),
                   qbeta((beta_cdf(0.6,aGP,bGP)+(1-beta_cdf(0.6,aGP,bGP))*0.5),aGP,bGP)]';
int counter = 0;

for (i in 42:49) {
   counter += 1;
   Y[ : , i ] ~ normal(cache[counter], sigma_e);
 }
}

The PR is still open (as I see)? I wonder if I could incorporate the code with derivatives into my program before it is added to Stan?

The run-time went down from 25 mins to 5 mins! Is it because of declaring the variable “cache” in advance (so that it is allocated differently in memory)? Thanks for the tip!

Partially but it’s really due to computing qbeta. Before you were computing the same output N times and we know that qbeta is slow. Precomputing it allows the program to just access the value when needed for the n^th row

Oh right! I was translating my code from JAGS, made a blunder!