# Other libraries of Distribution Functions

## Hypergeometric distribution

### Octave

Essentially, the formula used in the code is

`   cdf(i) = discrete_cdf (x(i), v, hygepdf (v, t(i), m(i), n(i)));`

Essentially the formula used is

`  inv(i) = discrete_inv (x(i), v, hygepdf (v, t(i), m(i), n(i)));`

Essentially the formula used is

```        if (isscalar (t) && isscalar (m) && isscalar (n))
pdf(k) = (bincoeff (m, x(k)) .* bincoeff (t-m, n-x(k))
/ bincoeff (t, n));
else
pdf(k) = (bincoeff (m(k), x(k)) .* bincoeff (t(k)-m(k), n(k)-x(k))
./ bincoeff (t(k), n(k)));```

Essentially the formula used is

```      v = 0:n;
p = hygepdf (v, t, m, n);
rnd = v(lookup (cumsum (p(1:end-1)) / sum (p), rand (sz)) + 1);
rnd = reshape (rnd, sz);```

### Burkardt's codes

* The above 2 algorithms are matlab version of :

### R

Essentially, the formula used is

```    p = ((double)n)/((double)(r+b));
q = ((double)(r+b-n))/((double)(r+b));

p1 = dbinom_raw(x,  r, p,q,give_log);
p2 = dbinom_raw(n-x,b, p,q,give_log);
p3 = dbinom_raw(n,r+b, p,q,give_log);

return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );```

Algorithm used:

```/*
* Calculate
*
*          phyper (x, NR, NB, n, TRUE, FALSE)
*   [log]  ----------------------------------
*             dhyper (x, NR, NB, n, FALSE)
*
* without actually calling phyper.  This assumes that
*
*     x * (NR + NB) <= n * NR
*
*/```

Algorithm Coded

```    d  = dhyper (x, NR, NB, n, log_p);
pd = pdhyper(x, NR, NB, n, log_p);

return log_p ? R_DT_Log(d + pd) : R_D_Lval(d * pd);```

Algorithm Used

```    /* Goal:  Find  xr (= #{red balls in sample}) such that
*   phyper(xr,  NR,NB, n) >= p > phyper(xr - 1,  NR,NB, n)
*/```

``` *  REFERENCE
*
*    V. Kachitvichyanukul and B. Schmeiser (1985).
*    ``Computer generation of hypergeometric random variates,''
*    Journal of Statistical Computation and Simulation 22, 127-145.```

## Geometric Distribution

### R

The help page is at :

Essentially, the source uses the formula

```/* prob = (1-p)^x, stable for small p */
prob = dbinom_raw(0.,x, p,1-p, give_log);
return p*prob;```

Essentially, the formula used is

```x = log1p(-p) * (x + 1);
return lower_tail ? -expm1(x) : exp(x);```

Essentially, the formula which is used is :

```/* add a fuzz to ensure left continuity */
ceil(R_DT_Clog(p) / log1p(- prob) - 1 - 1e-7)```

This formula is perhaps the cause of the accuracy bug report:

``` *    We generate lambda as exponential with scale parameter
*    p / (1 - p).  Return a Poisson deviate with mean lambda.
rpois(exp_rand() * ((1 - p) / p));```

### Octave

Essentially, the code uses the formula :

`cdf(k) = 1 - ((1 - p(k)) .^ (x(k) + 1));`

Essentially, the code uses the formula :

`inv(k) = max (ceil (log (1 - x(k)) ./ log (1 - p(k))) - 1, 0);`

`pdf(k) = p(k) .* ((1 - p(k)) .^ x(k));`

`rnd = floor (- rande (sz) ./ log (1 - p));`

## Poisson Distribution

### R

The source uses the paper :

``` *    Ahrens, J.H. and Dieter, U. (1982).
*    Computer generation of Poisson deviates
*    from modified normal distributions.
*    ACM Trans. Math. Software 8, 163-179.```

2022-09-08 09:26