[Contents] [TitleIndex] [WordIndex

Other libraries of Distribution Functions

Hypergeometric distribution

Stixbox

Octave

http://fossies.org/dox/octave-3.6.2/hygecdf_8m_source.html

Essentially, the formula used in the code is

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

http://fossies.org/dox/octave-3.6.2/hygeinv_8m_source.html

Essentially the formula used is

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

http://fossies.org/dox/octave-3.6.2/hygepdf_8m_source.html

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

http://fossies.org/dox/octave-3.6.2/hygernd_8m_source.html

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

http://people.sc.fsu.edu/~jburkardt/m_src/asa152/hypergeometric_cdf_values.m

http://people.sc.fsu.edu/~jburkardt/m_src/asa152/chyper.m

* The above 2 algorithms are matlab version of :

http://lib.stat.cmu.edu/apstat/152

R

http://svn.r-project.org/R/trunk/src/nmath/dhyper.c

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

http://svn.r-project.org/R/trunk/src/nmath/phyper.c

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

http://svn.r-project.org/R/trunk/src/nmath/qhyper.c

Algorithm Used

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

http://svn.r-project.org/R/trunk/src/nmath/rhyper.c

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

Boost

http://www.boost.org/doc/libs/1_50_0_beta1/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/hypergeometric_dist.html

Geometric Distribution

Stixbox

Boost

http://www.boost.org/doc/libs/1_50_0_beta1/libs/math/doc/sf_and_dist/html/math_toolkit/dist/dist_ref/dists/geometric_dist.html

R

The help page is at :

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/Geometric.html

http://svn.r-project.org/R/trunk/src/nmath/dgeom.c

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;

http://svn.r-project.org/R/trunk/src/nmath/pgeom.c

Essentially, the formula used is

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

http://svn.r-project.org/R/trunk/src/nmath/qgeom.c

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:

https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=14967

 *    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

http://fossies.org/dox/octave-3.6.2/geocdf_8m_source.html

Essentially, the code uses the formula :

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

http://fossies.org/dox/octave-3.6.2/gaminv_8m_source.html

Essentially, the code uses the formula :

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

http://fossies.org/dox/octave-3.6.2/geopdf_8m_source.html

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

http://fossies.org/dox/octave-3.6.2/geornd_8m_source.html

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

Poisson Distribution

R

http://svn.r-project.org/R/trunk/src/nmath/rpois.c

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.

F distribution

Boost

Burkardt

Stixbox

casci

Other sources

http://www.gnu.org/software/octave/doc/interpreter/Distributions.html


2022-09-08 09:26