The matrix k contains counts in the form of a contingency table. More about
that in a moment.
The rowSums and colSums functions compute the row and column sums of the
matrix k. The result is a vector.
When R applies log to a matrix, it does element-wise evaluation of the
function to each element of the matrix, likewise with * (matrix
multiplication is %*%). The sum function adds up all elements.
The attached image has it in mathematical notation.
I have attached a Java implementation as well (which you are free to use,
but which needs to have its dependency on Jama removed).
On Mon, Jun 2, 2008 at 6:08 PM, Sean Owen <[EMAIL PROTECTED]> wrote:
> Sorry I am still not clear on the algorithm. What is rowSums,
> columnSums, and what is k? k seem to be used as both a scalar, and
> vector (argument to a function sum()?)
>
> Maybe it is just that I don't know R. Can you give me a hint or
> forward any Java implementations? I can figure it out from there.
> Don't want to be a bother but I do want to make sure I understand what
> you're saying.
>
> On Mon, Jun 2, 2008 at 8:56 PM, Ted Dunning <[EMAIL PROTECTED]> wrote:
> > Hmmm... may not have sent that to you.
> >
> > function(k) {
> > return (2 * (h(k) - h(rowSums(k)) - h(colSums(k))));
> > }
> >
> > function(k) {
> > p = k / sum(k) + (k == 0);
> > return (sum(k * log(p)));
> > }
> >
> > This is the pithiest implementation that I know of.
> >
> > On Mon, Jun 2, 2008 at 5:50 PM, Sean Owen <[EMAIL PROTECTED]> wrote:
> >
> >> PS what's the R code? sorry if I am missing something basic.
> >>
> >> On Sun, Jun 1, 2008 at 8:32 PM, Ted Dunning <[EMAIL PROTECTED]>
> wrote:
> >> > The R code is instructive here. Everywhere you are about to take the
> log
> >> of
> >> > a zero probability, you should add 1. This is equivalent to saying
> that
> >> >
> >> > lim_{x -> 0} x log x = 0
> >>
> >> OK true but...
> >>
> >> > One easy way to do this is to wrap your log function:
> >> >
> >> > safeLog = {x -> return (x==0)?0:Math.log(x)}
> >>
> >> ... log(x) approaches -infinity as x approaches 0 so why is 0 a good
> >> substitute for log(0)?
> >>
> >> But I would buy that, for purposes of this function, "log(p)" can be
> >> replaced with "log(p+1)" so at least the output value is >= 0 and less
> >> than 1?
> >>
> >
> >
> >
> > --
> > ted
> >
>
--
ted
/**
* Copyright 2004, Ted Dunning
* Distribution allowed under GPL v2. Other terms available on request.
* @author Ted Dunning
*/
package com.tdunning.stat;
import Jama.Matrix;
/**
* Defines several matrix functions related to Shannon entropy.
*/
public class Entropy {
/**
* Returns the generalized log-likelihood ratio test for row and
* column independence for a contingency table k. Typically, k
* represents observed counts for multinomial processes. Each row
* represents a separate process and each column represents a
* symbol. This test gives an indication of how likely these
* observations were to have come from an identical multinomial
* distribution. The output in the case that the observations
* come from the same distribution is asymptotically chi^2
* distributed with (n-1) (m-1) degrees of freedom where n is the
* number of rows and m the number of columns.
* @param k
* @return The result of the llr test.
*/
public static double llrMultinomial(Matrix k) {
return 2 * (entropy(k.rowSum()) + entropy(k.columnSum()) - entropy(k));
}
/**
* Returns a signed version of the generalized log-likelihood ratio test
* for a two-dimensional table of counts. Each row contains counts for
* a different context that some condition might hold, each column contains
* the number of times the condition held or not. The result of this
* test is positive if the first row of the table had more times where
* the condition held than would be expected compared to the second row.
*
* Fractional counts are allowed because the counts may have been estimated
* by some odd process rather than actually counted.
* @param k1a Number of times condition or event "a" occured in context 1.
* @param k1notA Number of times something other than "a" occured in context 1.
* @param k2a Number of times condition or event "a" occured in context 2.
* @param k2notA Number of times something other than "a" occured in context 2.
* @return A signed score that is positive when "a" occurs more often in condition 1 than in condition 2. If
* the probability of getting "a" is the same in either condition and whatever difference we see is due purely
* to chance, then this score will be asymptotically normally distributed with mean 0 and variance 1.
*/
public static double signedLlrFrequencyTest(double k1a, double k1notA, double k2a, double k2notA) {
double llr = Math.sqrt(llrMultinomial(new Matrix(new double[]{k1a, k2a, k1notA, k2notA}, 2)));
if (k1a / (k1a + k1notA) > k2a / (k2a + k2notA)) {
return llr;
} else {
return -llr;
}
}
/**
* Returns the generalized log-likelihood ratio test for a
* contingency table k with row times in the last column.
* Typically, k represents observed counts for mixed Poisson
* processes. Each row represents a separate process and each
* column represents a symbol except for the last which contains
* the total observation time for the row. This test gives an
* indication of how likely these observations were to have come
* from an identical mixed Poisson distribution. This test can be
* used on normal Poisson processes by giving it a matrix with
* exactly two columns. The output in the case that the
* observations are taken from the same mixed Poisson distribution
* is chi^2 distributed with (n-1) (m-1) degrees of freedom where
* n is the number of rows and m is the number of columns in the
* matrix (note that one of these columns is not a countAll).
* @param k
* @return The results of the llr test.
*/
public static double llrPoisson(Matrix k) {
int rows = k.getRowDimension();
int cols = k.getColumnDimension();
Matrix counts = k.getMatrix(0, rows-1, 0, cols-2);
Matrix times = k.getMatrix(0, rows-1, cols-1, cols-1);
return 2 * (entropy(times) + entropy(k.columnSum()) - entropy(counts));
}
/**
* Computes the entropy of a matrix. Strictly speaking, this
* returns sum_ij k_ij log p_ij where p_ij = k_ij / sum_ij k_ij.
* If sum_ij k_ij = 1, then this does no harm and it is convenient
* other computations.
*
* @param k
* @return The entropy in base e. Divide by log(2) to get bits.
*/
public static double entropy(Matrix k) {
double total = k.columnSum().rowSum().get(0, 0);
if (k.columnMin().rowMin().get(0, 0) < 0) {
throw new IllegalArgumentException("Cannot have negative count");
}
if (total <= 0) {
throw new IllegalArgumentException("Must have non-zero count");
}
// normalize p and add one where zero. Adding one will have no effect since we
// multiply by the original value later. Doing this just avoids log(0) errors.
Matrix p = k.times(1/total).plusEquals(k.lessEqual(0.0));
return -k.arrayTimes(p.log()).rowSum().columnSum().get(0, 0);
}
}
