LRbinom.c

Go to the documentation of this file.

/*
 * Copyright    2019    R.K. Owen, Ph.D.
 * License      see lgpl.md (Gnu Lesser General Public License)
 */
#ifdef __cplusplus
extern "C" {
#endif

#include <math.h>
#include "libran.h"

/* int */
int LRi_binomial_RAN(LR_obj *o) {
        float   zero = 0.0, one = 1.0, u, sum = zero, p = one;
        int     i;
        if (isnan(o->q))
                o->q = -logf(one - o->p);

        for (i = 0; i < o->n; i++) {
                do {
                        u = o->uf(o);
                } while (u == zero);
                sum += (-logf(u)/(o->n - i));
                if (sum > o->q)
                        break;
        }

        return i;
}

float LRi_binomial_PDF(LR_obj *o, int k) {
        float zero = 0.0, one = 1.0, p = o->p, lp, lq, lt;

        if (isnan(o->q))
                o->q = -logf(one - o->p);

        if (k < 0 || k > o->n)  return zero;

        lp = logf(o->p);
        lq = - o->q;

        lt = k*lp + (o->n - k)*lq;
        lt += (lgammaf(o->n + 1) - lgammaf(k + 1) - lgammaf(o->n - k + 1));

        return expf(lt);
}

float LRi_binomial_CDF(LR_obj *o, int k) {
        float zero = 0.0, one = 1.0, sum = zero, p = o->p, lp, lq, lt;

        if (isnan(o->q))
                o->q = -logf(one - o->p);

        if (k < 0 || k > o->n)  return zero;

        lp = logf(o->p);
        lq = - o->q;

        for (int i = 0; i <= k; i++) {
                lt = i*lp + (o->n - i)*lq;
                lt += (lgammaf(o->n + 1) - lgammaf(i+1) - lgammaf(o->n - i+1));
                sum += expf(lt);
        }

        return sum;
}

#ifdef __cplusplus
}
#endif