Other libraries of Distribution Functions
Contents
Hypergeometric distribution
Stixbox
http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/dhypg.sci
http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/phypg.sci
http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/rhypg.sci
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
PDF : http://www.boost.org/doc/libs/1_48_0/boost/math/distributions/detail/hypergeometric_pdf.hpp
CDF : http://www.boost.org/doc/libs/1_48_0/boost/math/distributions/detail/hypergeometric_cdf.hpp
Inverse CDF : http://www.boost.org/doc/libs/1_48_0/boost/math/distributions/detail/hypergeometric_quantile.hpp
Geometric Distribution
Stixbox
Boost
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
PDF : http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/df.sci
CDF : http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/pf.sci
Inverse CDF : http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/qf.sci
RNG : http://forge.scilab.org/index.php/p/stixbox/source/tree/master/macros/rf.sci
casci
PDF : http://forge.scilab.org/index.php/p/casci/source/tree/master/macros/pdffisher.sci
CDF : http://forge.scilab.org/index.php/p/casci/source/tree/master/macros/cdffisher.sci
Inverse CDF : http://forge.scilab.org/index.php/p/casci/source/tree/master/macros/idffisher.sci
Other sources
http://www.gnu.org/software/octave/doc/interpreter/Distributions.html