Drum roll please....
I figured it out!!
Gave up translating strange engineering code (read LISP and Fortran) and decided to do it with pencil and paper. Too make a long story short, I poked around on the net and found some course notes for a Math class on linear equations at a graduate Math program in New Zealand. I then wrote out, using the aid of those notes, a solution to the problem using a sample data set. I have spent the last 3 days turning those 20+ pages of iterative math into less than 50 lines of PHP code.
And here it is, free for all to use however you want. I hope you can make some money from this (as I am), but give a little credit (to ynoth) where credit is due if you use it. Noone can copyright math, no matter what the crazies behind the DMCA have to say, so this is free for all to use as they choose.
I'm working on a class that incorporates this algorithm into something more substantial. But that's a while (days not weeks) away from an initial release. I want to include a Jacobi Method option in addition to the Gauss-Seidel Method. In addition, functions for matrix validation, matrix modification, error reporting, etc. are in various stages of completion.
I have chosen to not use the Matrix class in the Math_Matrix PEAR package as a foundation. It has some serious bugs/problems regarding Matrix multiplication among other things. While I think you could hack it together, I see no reason or obligation to base this working algorithm on another, buggy, Math package.
In addition, I'm putting together a series of notes on practical applications of this code. Hopefully, some or all of those notes will accompany the class itself. For future work, I'm extending the algorithm to solve for certain variations of the Bradley-Terry model, such as those involving multiple (non-pairwise) comparisons, weighting factors, and non-exclusive comparison results (draws/ties/equivocations/etc.).
And so, here is the just the basic algorithm piece in case anyone wants to see it or use it. (I wish I could get my little sigmas in the comment to show up in the forum, but you'll have settle for SUM instead.)
Code: Select all
/**
* @Author ynoth
* Remove author credit and suffer pain of karmic justice.
*
* Iterative solution of the Bradley-Terry Model using the Gauss-Seidel Method.
* Bradley-Terry Model:
* Ri = Wi / (SUMj Nij / (Ri + Rj))
* for i = 1, 2... n
* where i != j
* and where n = the number of columns of the pairwise matrix
*
* Gauss-Seidel Method applied to the Bradley-Terry Model:
* (k+1) (k) (k)
* Ri = Wi / (SUMj Nij / ( Ri + Rj ))
* for i = 1, 2... n
* where i != j
* and where n = the number of columns of the pairwise matrix
* and where k = 1, 2... kmax
*
* Convergence is obtained if:
* (k) (k)
* epsilon > max | Ri - Rj |
* for all i,j = 1, 2... n
* where i != j
* and where n = the number of columns of the pairwise matrix
*
* Matrix notes:
* Must be a matrix containing the results of pairwise comparisons and
* therefore must be square (#rows = #cols).
* Not all matrices will converge. See Bradley-Terry assumptions for
* specifics.
*
* Sample matrix form (this matrix will converge).
* array(0 => array( 0, 3, 1),
* 1 => array(0, 0, 2),
* 2 => array(1, 1, 0));
*
* Anticipated results for the above sample matrix:
* (2.280469, 0.456094, 0.760156)
*/
function solveModel($matrix,
$k_max = 1000,
$epsilon = .0000001,
$Ri_estimate = array())
{
$matrix_cols = count($matrix);
// Populate $Ri_estimate if none was passed.
if (count($Ri_estimate) == 0) {
$Ri_estimate = array_fill(0, $matrix_cols, 1);
}
// Variables needed for the algorithm.
$row_sum = array();
foreach ($matrix as $key=>$val) {
array_push($row_sum, array_sum($val));
}
$Ri = $Ri_estimate;
$convergence = array_fill(0, $matrix_cols, 0);
for ($k = 1; $k <= $k_max; $k++) {
// Maintain two arrays to allow for convergence checks.
$Ri_last = $Ri;
for ($i = 0; $i < $matrix_cols; $i++) {
$sum = 0;
for ($j = 0; $j < $matrix_cols; $j++) {
// i != j
if ($i == $j) {
continue;
}
// SUMj Nij / (Ri + Rj)
$sum += ($matrixї$i]ї$j] + $matrixї$j]ї$i]) /
($Riї$i] + $Riї$j]);
}
// Wi / (SUMj Nij / (Ri + Rj))
$Riї$i] = $row_sumї$i] / $sum;
// Monitor convergence.
if (abs($Ri_lastї$i] - $Riї$i]) < $epsilon) {
$convergenceї$i] = 1;
}
}
// Test for convergence after each complete iteration.
if (array_sum($convergence) == $matrix_cols) {
break;
}
// If there is no convergence after max_iterations, return false.
if ($k == $k_max) {
return FALSE;
}
}
return $Ri;
}
peace
eat more toast.