title: Continuous mutual information in Python
body: <p>[Frontmatter] <em>(skip this if you just want the question)</em>: </p>
  
  <p>I'm currently looking at using <a href="http://www.scholarpedia.org/article/Mutual_information">Shannon-Weaver Mutual Information</a> and <a href="https://en.wikipedia.org/wiki/Mutual_information#Normalized_variants">normalized redundancy</a> to measure the degree of information masking between bags of discrete and continuous feature values, organized by feature. Using this method, it is my goal to construct an algorithm that looks very similar to <a href="https://en.wikipedia.org/wiki/ID3_algorithm">ID3</a>, but instead of using <a href="https://en.wikipedia.org/wiki/Entropy_%28information_theory%29">Shannon entropy</a>, the algorithm will seek (as a loop constraint) to maximize or minimize shared information between a single feature and a collection of features based on the complete input feature space, adding new features to the latter collection if (and only if) they increase or decrease mutual information, respectively. This, in effect, moves ID3's decision algorithm into pairspace, stapling an ensemble approach to it with all of the expected time and space complexities of both methods.</p>
  
  <p>[/Frontmatter]</p>
  
  <p><br/></p>
  
  <p>On to the question: I'm trying to get a continuous <a href="https://en.wikipedia.org/wiki/Integral">integrator</a> working in Python using <a href="http://www.scipy.org/">SciPy</a>. Because I'm working with comparisons of discrete and continuous variables, my current strategy for each comparison for feature-feature pairs is as follows:</p>
  
  <ul>
  <li><p>Discrete feature versus discrete feature: use the discrete form of mutual information. This results in a double summation of the probabilities, which my code handles  without issue.</p></li>
  <li><p>All other cases (discrete versus continuous, the inverse, and continuous versus continuous): use the continuous form, using a <a href="https://en.wikipedia.org/wiki/Kernel_density_estimation">Gaussian estimator</a> to smooth out the <a href="https://en.wikipedia.org/wiki/Probability_density_function">probability density functions</a>.</p></li>
  </ul>
  
  <p>It is possible for me to perform some kind of discretization for the latter cases, but because my input data sets are not inherently linear, this is potentially needlessly complex.</p>
  
  <p><br/></p>
  
  <p>Here's the salient code:</p>
  
  <pre><code>import math
  import numpy
  import scipy
  from scipy.stats import gaussian_kde
  from scipy.integrate import dblquad
  
  # Constants
  MIN_DOUBLE = 4.9406564584124654e-324 
                      # The minimum size of a Float64; used here to prevent the
                      #  logarithmic function from hitting its undefined region
                      #  at its asymptote of 0.
  INF = float('inf')  # The floating-point representation for "infinity"
  
  # x and y are previously defined as collections of 
  # floating point values with the same length
  
  # Kernel estimation
  gkde_x = gaussian_kde(x)
  gkde_y = gaussian_kde(y)
  
  if len(binned_x) != len(binned_y) and len(binned_x) != len(x):
      x.append(x[0])
      y.append(y[0])
  
  gkde_xy = gaussian_kde([x,y])
  mutual_info = lambda a,b: gkde_xy([a,b]) * \
             math.log((gkde_xy([a,b]) / (gkde_x(a) * gkde_y(b))) + MIN_DOUBLE)
  
  # Compute MI(X,Y)
  (minfo_xy, err_xy) = \
      dblquad(mutual_info, -INF, INF, lambda a: 0, lambda a: INF)
  
  print 'minfo_xy = ', minfo_xy
  </code></pre>
  
  <p>Note that overcounting exactly one point is done deliberately to prevent a singularity in SciPy's <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html">gaussian_kde</a> class. As the size of x and y mutually approach infinity, this effect becomes negligible. </p>
  
  <p><br/></p>
  
  <p>My current snag is in trying to get <a href="http://en.wikipedia.org/wiki/Multiple_integral">multiple integration</a> working against a <a href="https://en.wikipedia.org/wiki/Kernel_density_estimation">Gaussian kernel density estimate</a> in SciPy. I've been trying to use SciPy's <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html"><em>dblquad</em></a> to perform the integration, but in the latter case, I receive an astounding spew of the following messages. </p>
  
  <p>When I set <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.seterr.html"><code>numpy.seterr ( all='ignore' )</code></a>:</p>
  
  <blockquote>
    <p>Warning: The ocurrence of roundoff error is detected, which prevents
         the requested tolerance from being achieved.  The error may be
         underestimated.</p>
  </blockquote>
  
  <p>And when I set it to <code>'call'</code> using an error handler:</p>
  
  <blockquote>
    <p>Floating point error (underflow), with flag 4</p>
    
    <p>Floating point error (invalid value), with flag 8</p>
  </blockquote>
  
  <p>Pretty easy to figure out what's going on, right? Well, almost: <a href="https://en.wikipedia.org/wiki/IEEE_754-2008">IEEE 754-2008</a> and SciPy only tell me what's going on here, not <em>why</em> or <em>how to work around it</em>.</p>
  
  <p><br/></p>
  
  <p>The upshot: <code>minfo_xy</code> generally resolves to <code>nan</code>; its sampling is insufficient to prevent information from becoming lost or invalid when performing Float64 math.</p>
  
  <p>Is there a general workaround for this problem when using SciPy? </p>
  
  <p>Even better: <em>if there is a robust, canned implementation of continuous mutual information for Python with an interface that takes two collections of floating point values or a merged collection of pairs, it would resolve this complete problem. Please link it if you know of one that exists.</em> </p>
  
  <p>Thank you in advance.</p>
  
  <hr>
  
  <p><strong>Edit:</strong> this resolves the <code>nan</code> propagation issue in the example above:</p>
  
  <pre><code>mutual_info = lambda a,b: gkde_xy([a,b]) * \
      math.log((gkde_xy([a,b]) / ((gkde_x(a) * gkde_y(b)) + MIN_DOUBLE)) \
          + MIN_DOUBLE)
  </code></pre>
  
  <p>However, the question of roundoff correction remains, as does the request for a more robust implementation. Any help in either domain would be greatly appreciated.</p>