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