At the moment, our Kraskov calculators assume that they have k values within 
epsilon in the full joint space. Where the distribution is composed of delta 
functions (i.e. somewhat discretised) this won't be the case - there could be 
many more than k. As such, it may be more correct to use the actual k for each 
observation in contributing digamma(k) - 1/k, then averaging, rather than 
assume the parameter k.
This may add significantly to the runtime however so we could either:
- just issue a warning were this to occur, or
- track the actual k values in an array, then compute digamma(k) - 1/k once for 
each distinct value and add these in to the result.