0

我从 C 中的数值配方中获得了以下代码，它使用连续分数和 Lentz 方法计算不完整的 beta 函数。

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);
    }
}

然后我写了一个主代码，我把 theta 作为输入，并得到一个不完整的 beta 函数的值。

现在我需要获得 theta = [0,1] 的分布。有没有办法以我不更改此代码中的任何内容的方式编写它。我的意思是在我的主函数中为 theta 添加一个 for 循环并获得不完整 beta 函数的输出。我试过这样做，但它抛出一个错误“不兼容的类型，预期的'double'但参数的类型是'double *'。我理解这个错误是因为我试图将输出作为一个数组但是在我的函数中它被定义了是一个单一的值。有没有一种解决方法，我不必在我的函数中将 theta 声明为数组。

主函数失败

 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]);
    }
}

谢谢

4

0 回答 0