I have the following code from the numerical recipes in C which calculates the incomplete beta function using continuous fraction and Lentz method.
float betacf(float m1, float m2, float theta){
void nrerror(char error_text[]);
int k, k2, MAXIT;
float aa, c, d, del, t, qab, qam, qap;
qab = m1 + m2;
qap = m1 + 1.0;
qam = m1 - 1.0;
c = 1.0;
d = 1.0 - (qab * theta)/qap;
if (fabs(d) < FPMIN) d = FPMIN;
d = 1.0/d;
t = d;
for (k = 1; k <= MAXIT; k++) {
k2 = 2 * k;
aa = k * (m2 - k) * theta/((qam + k2) * (m1 + k2));
d = 1.0 + aa * d;
if (fabs(d) < FPMIN) d = FPMIN;
c = 1.0 + aa/c;
if (fabs(c) < FPMIN) c = FPMIN;
d = 1.0/d;
t *= d * c;
aa = -(m1 + k) * (qab + k) * theta/((m1 + k2) * (qap + k2));
d = 1.0 + aa * d;
if (fabs(d) < FPMIN) d = FPMIN;
c = 1.0 + aa/c;
if (fabs(c) < FPMIN) c=FPMIN;
d = 1.0/d;
del = d * c;
t *= del;
if (fabs(del - 1.0) < EPS) break;
}
if (k > MAXIT) nrerror("m1 or m2 too big, or MAXIT too small in betacf");
return t;
}
/* Returns the incomplete beta function Ix(a, b) */
float betai(float m1, float m2, float theta){
void nrerror(char error_text[]);
float bt;
if (theta < 0.0 || theta > 1.0){
nrerror("Bad x in routine betai");
}
if (theta == 0.0 || theta == 1.0){
bt = 0.0;
}
else {
bt = exp(gammaln(m1+m2)-gammaln(m1)-gammaln(m2)+m1*log(theta)+m2*log(1.0-theta));
}
if (theta < (m1 + 1.0)/(m1 + m2 + 2.0))
{
return (bt * betacf(m1, m2, theta)/m1);
}
else {
return (1.0 - bt * betacf(m2, m1, 1.0 - theta)/m2);
}
}
Then I write a main code where I throw in theta as input and get a value for incomplete beta function.
Now I need to obtain a distribution for theta = [0,1]. Is there a way to write it in way where I don't change anything in this code. I mean just add a for loop in my main function for theta and get the output of the incomplete beta function. I tried doing this, but it throws an error "Incompatible types, expected 'double' but argument is of type 'double *' . I understand the error is because I try to get the output as an array but in my function it is defined to be a single value. Is there a work around this where I don't have to declare theta as an array in my function.
Failing main function
int main() {
float *theta, *result;
.....
.....
printf("Enter number of points required to describe the PDF profile:", N);
scanf("%d", &N);
theta = (float *)malloc(N*sizeof(float));
for (j = 1; j < N; j++)
theta[j] = (float)(j)/ ((float)(N) - 1.0);
result[j] = betai(m1, m2, theta);
printf("%f %f", theta[j], result[j]);
}
}
Thank you