Implementing Bayesian Inference Using PHP: Part 2 - Graphing the likelihood distribution

(Page 5 of 10 )

You've computed the MLE of p by applying the max()function to the likelihood array (a more accurate root-finding method such as Newton-Raphson should be used when more accuracy is required). Selecting the max of the likelihood array only gives you a point estimate of what value of p is most likely given the results. Other values of p are also possible but less likely to produce the observed results.

To get a sense of how likely various values of p are, we should examine the relationship between different p values and their corresponding likelihood values (denoted as l(p) in the following graph).

JPGraph is the leading PHP-based package for creating professional graphs that can be displayed on the Web or in other media. The following code is used to create a graph of the likelihood distribution for p. One interesting feature of this code is that it demonstrates usage of the cubic spline function for interpolating values and creating smooth curves.

Listing 3. Creating the likelihood distribution graph

  <?php
/**
* Script to graph likelihood distribution of the binomial
* parameter p (for example, probability of success per trial).
*
* Much of the plotting code below has been adapted from
* example code written by JPGraph author Johan Persson.
*
* To get this working locally, you must install JPGraph
* and set the PHP_MATH constant to the folder where the
* JPGraph source code resides.
*
* @see http://www.aditus.nu/jpgraph/index.php
*/

require_once "../config.php";

include_once JPGRAPH . "/src/jpgraph.php";
include_once JPGRAPH . "/src/jpgraph_line.php";
include_once JPGRAPH . "/src/jpgraph_scatter.php";
include_once JPGRAPH . "/src/jpgraph_regstat.php";

include_once "../functions/binomial.php";

$n = 5; // num events
$k = 1; // num success events

$i = 0; // counter
for($p = 0.00; $p <= 1.00; $p += 0.05 ) {
  $likelihoods[$i] = binomial($n, $k, $p);
  $parameters[$i]  = $p;
  $i++;
} 

$mle = max($likelihoods);
$p   = $parameters[array_search($mle, $likelihoods)];

$spline = new Spline($parameters, $likelihoods);
list($newx,$newy) = $spline->Get(50);

$graph = new Graph(450, 350);
$graph->SetMargin(60, 20, 40, 30);

$graph->title->Set("Maximum Likelihood Estimate");
$graph->subtitle->Set(" MLE = P( $k/$n | $p ) = $mle ");
$graph->subtitle->SetColor('darkred');

$graph->SetMarginColor('lightblue');

$graph->SetScale('linlin');

$graph->xaxis->SetLabelFormat('%1.2f');
$graph->yaxis->SetLabelFormat('%1.2f');

$graph->xaxis->SetTitle("p","center");
$graph->yaxis->SetTitleMargin(40);
$graph->yaxis->SetTitle("l(p)", "center");

$graph->xaxis->title->SetFont(FF_FONT1,FS_BOLD);
$graph->yaxis->title->SetFont(FF_FONT1,FS_BOLD);

$splot = new ScatterPlot($likelihoods, $parameters);
$splot->mark->SetFillColor('red@0.3');
$splot->mark->SetColor('red@0.5');

$lplot = new LinePlot($newy,$newx);
$lplot->SetColor('navy');

$graph->Add($lplot);
$graph->Add($splot);
$graph->Stroke();

?>

 

This code produces the following graph:

Figure 1. The likelihood distribution graph

The y-axis represents the likelihood of p (denoted l(p)) and was computed using the binomial formula. The subtitle of the graph tells you that the likelihood achieves a maximum (technically where the derivative is 0) of .4096 when p = 0.20, which is equal to the observed proportion of sample successes.

Why do you need to do all this work to estimate p when we could have used common sense to arrive at the same result? The fact that the MLE procedure agrees with common sense helps to convince you that you can also use this technique when estimating parameters that are not so easy to determine through common sense. In those cases you can proceed by:

1. Finding a way to express the likelihood of the results as a function of the parameters.
   
2. Computing the likelihood of the result with respect to the parameter.
   
3. Selecting the parameter value that yields the maximum likelihood value.

The use of maximum likelihood principle is pervasive in statistical reasoning along with other principles such as the Bayesian principle of maximizing the posterior. Other notable principles include the least-squared error criterion, maximum entrophy, minimum description-length, variational inference, and various energy minimization principles.

High-level statistical reasoning is more about welding these principles effectively to estimate parameters, test hypothesis, and such, than it is about algebraic cleverness in deriving new formulas. A bit of algebraic cleverness, however, can come in handy.

Next: Algebraic cleverness >>
 

More PHP Articles
More By developerWorks

Comments On: Implementing Bayesian Inference Using PHP: Part 2

· Who uses Bayesian inference these days? It is computationally harder,and it also...    · what do they use then? i remember that this, was the first things i learned in my...   >>> Post your comment now!