Regularized incomplete beta function implementation

Modeling
1

Hi,
I’m trying to fit a model where we assume that the data have a relationship based on the regularized incomplete beta function.

y(x)=kI_{cx}(a,b), where k, c, a and b are parameters used for the fit.

Is there an existing implementation of this function in Stan? I know there is an inc_beta function but this doesn’t appear to work for vector input. At the moment I’m running my own metropolis samplers to do this, but I’m trying to generate model weights and WAIC/LOO seem to be pretty unstable with my homemade samplers.

2

I don’t think we have that specific function, but you can compute your likelihood manually using the 1d integrator, docs here: https://mc-stan.org/docs/2_21/stan-users-guide/integrate-1d.html and here: https://mc-stan.org/docs/2_21/functions-reference/functions-1d-integrator.html

3

I have a solution that I think should work which is to pass each element in my data vector to inc_beta and divide by the beta function which should be inc_beta(a,b,1). Unfortunately when I do so, I get initialization failed, even when I initialize at the ‘real’ parameter values for a test dataset, any idea why this would be the case? My model is below.

data {
    int<lower=0> N;
    real y[N];
    real x[N];
    real sigma;
}
parameters {
    real<lower=0,upper=1> U;
    real<lower=0> V;
}
transformed parameters {
    real<lower=0> a= U*V;
    real<lower=0> b= V-a;
    real y_pred[N];
    for (n in 1:N) {
      y_pred[n]= inc_beta(a,b,x[n])/inc_beta(a,b,1);
    } 
}
model {
    U ~ uniform(0,1);
    V ~ gamma(1,20);
    y ~ normal(y_pred,sigma);
}
4

Lol oops I did a poor job reading your question and the docs.

I guess make sure all the constraints are satisfied.

Is x greater than zero and less than one? If there are values exactly equal to 0 and 1 maybe force them to be greater than 0 or less than one and see if the error goes away.

5

0<x<1 is definitely satisfied for the test data I am generating. I am still getting an initialization failure. Not sure as to the cause of this.

6

Can you post the model/data?

It could be difficulty evaluating one of the inc_beta terms. Can you try commenting them out and seeing if it initializes?

7

The data I’m currently testing this on are simply xs generated from a uniform distribution on the range 0 to 1 and ys generated from exactly the model here with known a and b parameters.

It seems that dividing by the complete beta function causes the problem as the model runs if I take this out (although the model is then wrong). Not sure why this would be the case, but perhaps it’s because the denominator here has inc_beta(a,b,1), which is on the edge of the acceptable range for this function. Interestingly, if I set to a value close to 1 but not quite at 1 I get a different error:
“k (internal counter) exceeded 100000 iterations, hypergeometric function gradient did not converge.”

8

Oh okay I guess write it in terms of gamma functions then? Looks like that’s possible.