原始 C 中的不完整 beta 函数答案

【问题标题】：Incomplete beta function in raw C原始 C 中的不完整 beta 函数
【发布时间】：2012-06-07 07:58:03
【问题描述】：

我的一个朋友需要 MatLAB 的 betainc 函数的类似物，用于可编程逻辑器件 (PLD) 中的一些统计计算（我不是硬件人，还不知道他项目的任何细节）。

因此，不能选择使用预编译库。考虑到三个参数中的每一个都是可变的，她需要在原始 C 中实现。

网络上有什么好的吗？

非常感谢您！

【问题讨论】：

查看 C++/TR1 的 math.h 中的实现并在必要时复制它。移植到 C 应该相当容易。
或者你可以从matlab源码中复制betainc.m的算法。不是那么容易，但更有启发性。
这是一个很难实现的函数，具体取决于参数。如果您对需要支持的 beta 发行版的参数范围有所了解，则可以删减很多代码。

标签： c function math statistics implementation

【解决方案1】：

我知道我回答迟了，但您当前接受的答案（使用“数字食谱”中的代码）的许可证很糟糕。此外，它不会帮助尚未拥有这本书的其他人。

这是在 Zlib 许可下发布的不完整 beta 功能的原始 C99 代码：

#include <math.h>

#define STOP 1.0e-8
#define TINY 1.0e-30

double incbeta(double a, double b, double x) {
    if (x < 0.0 || x > 1.0) return 1.0/0.0;

    /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/
    if (x > (a+1.0)/(a+b+2.0)) {
        return (1.0-incbeta(b,a,1.0-x)); /*Use the fact that beta is symmetrical.*/
    }

    /*Find the first part before the continued fraction.*/
    const double lbeta_ab = lgamma(a)+lgamma(b)-lgamma(a+b);
    const double front = exp(log(x)*a+log(1.0-x)*b-lbeta_ab) / a;

    /*Use Lentz's algorithm to evaluate the continued fraction.*/
    double f = 1.0, c = 1.0, d = 0.0;

    int i, m;
    for (i = 0; i <= 200; ++i) {
        m = i/2;

        double numerator;
        if (i == 0) {
            numerator = 1.0; /*First numerator is 1.0.*/
        } else if (i % 2 == 0) {
            numerator = (m*(b-m)*x)/((a+2.0*m-1.0)*(a+2.0*m)); /*Even term.*/
        } else {
            numerator = -((a+m)*(a+b+m)*x)/((a+2.0*m)*(a+2.0*m+1)); /*Odd term.*/
        }

        /*Do an iteration of Lentz's algorithm.*/
        d = 1.0 + numerator * d;
        if (fabs(d) < TINY) d = TINY;
        d = 1.0 / d;

        c = 1.0 + numerator / c;
        if (fabs(c) < TINY) c = TINY;

        const double cd = c*d;
        f *= cd;

        /*Check for stop.*/
        if (fabs(1.0-cd) < STOP) {
            return front * (f-1.0);
        }
    }

    return 1.0/0.0; /*Needed more loops, did not converge.*/
}

取自Github repo。还有一篇关于how it works here的非常详尽的文章。

希望对您有所帮助。

【讨论】：

同意可怕的许可证。好书，糟糕的律师。

【解决方案2】：

或者您可以阅读“C 中的数字食谱”并找到完整的源代码。您将不得不担心许可问题，但它会清楚地解释该功能及其实现的含义。

【讨论】：