Implement Bayesian inference using PHP: Part 2, IBM developerWorks, April 13, 2004.
In this second article on Bayesian inference, Paul Meagher examines how you can use Bayes methods to solve parameter estimation problems. Relevant concepts are explained in the context of Web survey analysis using PHP and JPGraph.
Last updated: April 14, 2004 - 2:30 AM AST
Demos:
<?php
/**
* @author Jaco van Kooten
* @author Paul Meagher
* @license PHP v 3.0
* @version 0.1
*
* Port of JSci methods found in SpecialMath.java. This was originally
* a Java class - was easier to convert it to a function library for
* inclusion in PHP scripts. Note that newer versions of this port move
* many of the constants into the appropriate function instead of leaving
* them global.
*
* @see http://jsci.sourceforge.net/
*/
define("EPS", 2.22e-16);
define("MAX_VALUE", 1.2e308);
define("LOG_GAMMA_X_MAX_VALUE", 2.55e305);
define("SQRT2PI", 2.5066282746310005024157652848110452530069867406099);
define("XMININ", 2.23e-308);
define("MAX_ITERATIONS", 150);
define("PRECISION", 8.88E-016);
/*
* The SpecialMath lib implements the following functions.
*
* logBeta($p, $q)
* logGamma($x)
* Beta($p, $q)
* incompleteBeta($x, $p, $q)
* betaFraction($x, $p, $q)
*/
// Function cache for logBeta
$logBetaCache_res = 0.0;
$logBetaCache_p = 0.0;
$logBetaCache_q = 0.0;
/**
* The natural logarithm of the beta function.
* @param p require p>0
* @param q require q>0
* @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
* @author Jaco van Kooten
*/
function logBeta($p, $q) {
global $logBetaCache_res, $logBetaCache_p, $logBetaCache_q;
if ($p != $logBetaCache_p || $q != $logBetaCache_q) {
$logBetaCache_p = $p;
$logBetaCache_q = $q;
if ($p <= 0.0 || $q <= 0.0 || ($p + $q) > LOG_GAMMA_X_MAX_VALUE)
$logBetaCache_res = 0.0;
else
$logBetaCache_res = logGamma($p) + logGamma($q) - logGamma($p + $q);
}
return $logBetaCache_res;
}
/**
* logGamma function
*
* @version 1.1
* @author Jaco van Kooten
*
* Original author was Jaco van Kooten. Ported to PHP by Paul Meagher.
*
* The natural logarithm of the gamma function. <br />
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz <br />
* Applied Mathematics Division <br />
* Argonne National Laboratory <br />
* Argonne, IL 60439 <br />
* <p>
* References:
* <ol>
* <li>W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural
* Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.</li>
* <li>K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.</li>
* <li>Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.</li>
* </ol>
* </p>
* <p>
* From the original documentation:
* </p>
* <p>
* This routine calculates the LOG(GAMMA) function for a positive real argument X.
* Computation is based on an algorithm outlined in references 1 and 2.
* The program uses rational functions that theoretically approximate LOG(GAMMA)
* to at least 18 significant decimal digits. The approximation for X > 12 is from
* reference 3, while approximations for X < 12.0 are similar to those in reference
* 1, but are unpublished. The accuracy achieved depends on the arithmetic system,
* the compiler, the intrinsic functions, and proper selection of the
* machine-dependent constants.
* </p>
* <p>
* Error returns: <br />
* The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
* The computation is believed to be free of underflow and overflow.
* </p>
* @return MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
*/
// Log Gamma related constants
$lg_d1 = -0.5772156649015328605195174;
$lg_d2 = 0.4227843350984671393993777;
$lg_d4 = 1.791759469228055000094023;
$lg_p1 = array( 4.945235359296727046734888,
201.8112620856775083915565, 2290.838373831346393026739,
11319.67205903380828685045, 28557.24635671635335736389,
38484.96228443793359990269, 26377.48787624195437963534,
7225.813979700288197698961 );
$lg_p2 = array( 4.974607845568932035012064,
542.4138599891070494101986, 15506.93864978364947665077,
184793.2904445632425417223, 1088204.76946882876749847,
3338152.967987029735917223, 5106661.678927352456275255,
3074109.054850539556250927 );
$lg_p4 = array( 14745.02166059939948905062,
2426813.369486704502836312, 121475557.4045093227939592,
2663432449.630976949898078, 29403789566.34553899906876,
170266573776.5398868392998, 492612579337.743088758812,
560625185622.3951465078242 );
$lg_q1 = array( 67.48212550303777196073036,
1113.332393857199323513008, 7738.757056935398733233834,
27639.87074403340708898585, 54993.10206226157329794414,
61611.22180066002127833352, 36351.27591501940507276287,
8785.536302431013170870835 );
$lg_q2 = array( 183.0328399370592604055942,
7765.049321445005871323047, 133190.3827966074194402448,
1136705.821321969608938755, 5267964.117437946917577538,
13467014.54311101692290052, 17827365.30353274213975932,
9533095.591844353613395747 );
$lg_q4 = array( 2690.530175870899333379843,
639388.5654300092398984238, 41355999.30241388052042842,
1120872109.61614794137657, 14886137286.78813811542398,
101680358627.2438228077304, 341747634550.7377132798597,
446315818741.9713286462081 );
$lg_c = array( -0.001910444077728,8.4171387781295e-4,
-5.952379913043012e-4, 7.93650793500350248e-4,
-0.002777777777777681622553, 0.08333333333333333331554247,
0.0057083835261 );
// Rough estimate of the fourth root of logGamma_xBig
$lg_frtbig = 2.25e76;
$pnt68 = 0.6796875;
// Function cache for logGamma
$logGammaCache_res = 0.0;
$logGammaCache_x = 0.0;
function logGamma($x) {
global $lg_d1, $lg_d2, $lg_d4;
global $lg_p1, $lg_p2, $lg_p4;
global $lg_q1, $lg_q2, $lg_q4;
global $lg_c, $lg_frtbig, $pnt68;
global $logGammaCache_res;
global $logGammaCache_x;
if ($x == $logGammaCache_x) {
return $logGammaCache_res;
}
$y = $x;
if ($y > 0.0 && $y <= LOG_GAMMA_X_MAX_VALUE) {
if ($y <= EPS) {
$res = -log(y);
} else if ($y <= 1.5) {
// ---------------------
// EPS .LT. X .LE. 1.5
// ---------------------
if ($y < $pnt68) {
$corr = -log($y);
$xm1 = $y;
} else {
$corr = 0.0;
$xm1 = $y - 1.0;
}
if ($y <= 0.5 || $y >= $pnt68) {
$xden = 1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; $i++) {
$xnum = $xnum * $xm1 + $lg_p1[$i];
$xden = $xden * $xm1 + $lg_q1[$i];
}
$res = $corr + $xm1 * ($lg_d1 + $xm1 * ($xnum / $xden));
} else {
$xm2 = $y - 1.0;
$xden = 1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; $i++) {
$xnum = $xnum * $xm2 + $lg_p2[$i];
$xden = $xden * $xm2 + $lg_q2[$i];
}
$res = $corr + $xm2 * ($lg_d2 + $xm2 * ($xnum / $xden));
}
} else if ($y <= 4.0) {
// ---------------------
// 1.5 .LT. X .LE. 4.0
// ---------------------
$xm2 = $y - 2.0;
$xden = 1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; $i++) {
$xnum = $xnum * $xm2 + $lg_p2[$i];
$xden = $xden * $xm2 + $lg_q2[$i];
}
$res = $xm2 * ($lg_d2 + $xm2 * ($xnum / $xden));
} else if ($y <= 12.0) {
// ----------------------
// 4.0 .LT. X .LE. 12.0
// ----------------------
$xm4 = $y - 4.0;
$xden = -1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; $i++) {
$xnum = $xnum * $xm4 + $lg_p4[$i];
$xden = $xden * $xm4 + $lg_q4[$i];
}
$res = $lg_d4 + $xm4 * ($xnum / $xden);
} else {
// ---------------------------------
// Evaluate for argument .GE. 12.0
// ---------------------------------
$res = 0.0;
if ($y <= $lg_frtbig) {
$res = $lg_c[6];
$ysq = $y * $y;
for ($i = 0; $i < 6; $i++)
$res = $res / $ysq + $lg_c[$i];
}
$res /= $y;
$corr = log($y);
$res = $res + log(SQRT2PI) - 0.5 * $corr;
$res += $y * ($corr - 1.0);
}
} else {
// --------------------------
// Return for bad arguments
// --------------------------
$res = MAX_VALUE;
}
// ------------------------------
// Final adjustments and return
// ------------------------------
$logGammaCache_x = $x;
$logGammaCache_res = $res;
return $res;
}
/**
* Beta function.
*
* @author Jaco van Kooten
*
* @param p require p>0
* @param q require q>0
* @return 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
*/
function beta($p, $q) {
if ($p <= 0.0 || $q <= 0.0 || ($p + $q) > LOG_GAMMA_X_MAX_VALUE) {
return 0.0;
} else {
return exp(logBeta($p, $q));
}
}
/**
* Incomplete beta function
*
* @author Jaco van Kooten
* @author Paul Meagher
*
* The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
* @param x require 0<=x<=1
* @param p require p>0
* @param q require q>0
* @return 0 if x<0, p<=0, q<=0 or p+q>2.55E305 and 1 if x>1 to avoid errors and over/underflow
*/
function incompleteBeta($x, $p, $q) {
if ($x <= 0.0)
return 0.0;
else if ($x >= 1.0)
return 1.0;
else if ($p <= 0.0 || $q <= 0.0 || ($p + $q) > LOG_GAMMA_X_MAX_VALUE)
return 0.0;
else {
$beta_gam = exp( -logBeta($p, $q) + $p * log($x) + $q * log(1.0 - $x) );
if ($x < ($p + 1.0) / ($p + $q + 2.0))
return $beta_gam * betaFraction($x, $p, $q) / $p;
else
return 1.0 - ($beta_gam * betaFraction(1.0 - $x, $q, $p) / $q);
}
}
/**
* Evaluates of continued fraction part of incomplete beta function.
* Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
* @author Jaco van Kooten
*/
function betaFraction($x, $p, $q) {
$c = 1.0;
$sum_pq = $p + $q;
$p_plus = $p + 1.0;
$p_minus = $p - 1.0;
$h = 1.0 - $sum_pq * $x / $p_plus;
if (abs($h) < XMININ) {
$h = XMININ;
}
$h = 1.0 / $h;
$frac = $h;
$m = 1;
$delta = 0.0;
while ($m <= MAX_ITERATIONS && abs($delta-1.0) > PRECISION ) {
$m2 = 2 * $m;
// even index for d
$d = $m * ($q - $m) * $x / ( ($p_minus + $m2) * ($p + $m2));
$h = 1.0 + $d * $h;
if (abs($h) < XMININ) {
$h = XMININ;
}
$h = 1.0 / $h;
$c = 1.0 + $d / $c;
if (abs($c) < XMININ) {
$c = XMININ;
}
$frac *= $h * $c;
// odd index for d
$d = -($p + $m) * ($sum_pq + $m) * $x / (($p + $m2) * ($p_plus + $m2));
$h = 1.0 + $d * $h;
if (abs($h) < XMININ) {
$h = XMININ;
}
$h = 1.0 / $h;
$c = 1.0 + $d / $c;
if (abs($c) < XMININ) {
$c = XMININ;
}
$delta = $h * $c;
$frac *= $delta;
$m++;
}
return $frac;
}
?>