Kullback-Leibler divergence

//*****************************************************************************
// KLDivergence
// Created 20110702
// 
// Shawn Driscoll
// Salk Institute for Biological Studies
//
// Calculates the Kullback-Leibler divergence between the probability 
// distributions P and Q.  P and Q should have the same number of points
// and the distributions should have been created using a common min and max
// range to both.
//*****************************************************************************
function KLDivergence(P,Q)
    //--------------------------------------------------------------------------
    // parameters
    //--------------------------------------------------------------------------
    wave P,Q;
 
    //--------------------------------------------------------------------------
    // variables
    //--------------------------------------------------------------------------
    variable npnts = numpnts(P),i;
 
    //--------------------------------------------------------------------------
    // main function
    //--------------------------------------------------------------------------
 
    duplicate/free P,Pn;
    duplicate/free Q,Qn;
    
    ////
    // normalize distributions by their sums
    
    wavestats/q Pn;
    Pn /= v_sum;
    
    wavestats/q Qn;
    Qn /= v_sum;
 
    ////
    // create free wave for the results of the point by point
    // calculations
    make/free/n=(numpnts(P)) wresult = 0;
    
    ////
    // divide normalized P by normalized Q
    matrixop/free wdiv = Pn/Qn;
    
    ////
    // perform point by point calculations. I put this in a normal 
    // loop because for some reason the shorthand technique wasn't 
    // working.
    for(i=0; i<npnts; i+=1)
        wresult[i] = Pn[i]*(log(wdiv[i])/log(2));
        if(numtype(wresult[i]) > 0)
            wresult[i] = 0;
        endif
    endfor
    
    ////
    // find sum of result wave
    wavestats/q wresult;
 
    // set global var
    variable/g v_KLDivergence = v_sum;
    
    // return value
    return v_sum;
    
end