statsmodels: error in kde on a list of repeated values

Question

I am using Python's statsmodels module to plot a violin/bean plot of some data. I get the error "LinAlgError: singular matrix" in the KDE calculation whenever a single violin plot is drawn from a list of repeated instances of the same number. Regular boxplots work just fine.

Is this a bug or is there a reason behind this?

/home/me/myscript.py in plot_data_in_violinplot()
   3435         plot_opts={'violin_fc':(0.8, 0.8, 0.8), 'cutoff':True, 
                         'bean_color':'#FF6F00', 'bean_mean_color':'#009D91',
                         'bean_median_color':'b', 'bean_median_marker':'+'}
-> 3436         sm.graphics.beanplot(data_to_boxplot, ax=ax, jitter=jitter, plot_opts=plot_opts)


/usr/local/lib/python2.7/dist-packages/statsmodels/graphics/boxplots.py in beanplot(data, ax, labels, positions, side, jitter, plot_opts)
    333     for pos_data, pos in zip(data, positions):
    334         # Draw violins.
--> 335         xvals, violin = _single_violin(ax, pos, pos_data, width, side, plot_opts)
    336 
    337         if jitter:

/usr/local/lib/python2.7/dist-packages/statsmodels/graphics/boxplots.py in _single_violin(ax, pos, pos_data, width, side, plot_opts)
    170     pos_data = np.asarray(pos_data)
    171     # Kernel density estimate for data at this position.
--> 172     kde = gaussian_kde(pos_data)
    173 
    174     # Create violin for pos, scaled to the available space.

/usr/local/lib/python2.7/dist-packages/scipy/stats/kde.pyc in __init__(self, dataset, bw_method)
    186 
    187         self.d, self.n = self.dataset.shape
--> 188         self.set_bandwidth(bw_method=bw_method)
    189 
    190     def evaluate(self, points):

/usr/local/lib/python2.7/dist-packages/scipy/stats/kde.pyc in set_bandwidth(self, bw_method)
    496             raise ValueError(msg)
    497 
--> 498         self._compute_covariance()
    499 
    500     def _compute_covariance(self):

/usr/local/lib/python2.7/dist-packages/scipy/stats/kde.pyc in _compute_covariance(self)
    507             self._data_covariance = atleast_2d(np.cov(self.dataset, rowvar=1,
    508                                                bias=False))
--> 509             self._data_inv_cov = linalg.inv(self._data_covariance)
    510 
    511         self.covariance = self._data_covariance * self.factor**2

/usr/local/lib/python2.7/dist-packages/scipy/linalg/basic.pyc in inv(a, overwrite_a, check_finite)
    381         inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1)
    382     if info > 0:
--> 383         raise LinAlgError("singular matrix")
    384     if info < 0:
    385         raise ValueError('illegal value in %d-th argument of internal '
LinAlgError: singular matrix

Answer 1

I've never used the Python statsmodels package, and I'm not familiar with it, but based upon the error messages, I think I have a pretty good guess as to what is probably going on, and it's not a bug--the problem is with your input. According to wikipedia, a key step in Kernal Density Estimation is bandwidth estimation. As described in the link, for Gaussian basis functions (which, based upon your error message referencing gaussian_kde(), appears to be what you are in fact using), one common choice for estimating the bandwidth, in the special case of one dimensional KDE, requires the sample standard deviation $\hat{\sigma}$ as an input.

Your mention that the error arises, specifically, in cases where your input is a list of repeated instances of the same number. Imagine calculating the sample standard deviation for a list of identical numbers, what do you get? Well, essentially what you are simulating in that case is a dirac delta function, so effectively your sample standard deviation is $\hat{\sigma} = 0$. Based upon the fact the next error down in the stack trace, below the exception caught by gaussian_kde(), occurs within a method called set_bandwidth(), I would say that what appears to be happening is that you are feeding the code a distribution whose standard deviation is zero, and the code is attempting to use this value to calculate an initial guess for the KDE bandwidth parameter, and it's choking because zero isn't really a valid value.

"O.K.", you reply, "but your explanation doesn't mention anything about linear algebra or singular matrices--why does the error manifest itself within a linear algebra routine, specifically?" Good question. I'm not absolutely certain, but here's what I suspect is happening. The concept of standard deviation, or its square, the variance, is really an inherently one dimensional concept. The more general concept, valid for multi-variate distributions, is the covariance matrix. The code that you are using is likely designed to be as general as possible, in order to be able to handle a case where the user feeds it a multivariate distribution. In fact, you'll notice as you work your way further down the stack trace, that the next method down, below set_bandwidth(), is one called compute_covariance(). If you know much about covariance matrices, it turns out that a popular way of analyzing and thinking about them is to reduce them to what are known as principal components. The effect of principal component analysis, after it is performed, is to diagonalize the initial covariance matrix, creating an effectively equivalent new matrix which has been transformed in such a way that it consists exclusively of a set of one-dimensional variances lined up only along the diagonal. These variances may be identified as the eigenvalues of the original, non-diagonal matrix, and it also turns out that, in linear algebra, one of the properties of matrices is that those which are related by having identical eigenvalues also inevitably have identical determinants.

So, what I suspect is happening in your case is that, by giving the code repeated instances of the same values as input, you are creating a covariance matrix which has at least one eigenvalue equal to zero, and this condition means that the determinant is zero also, since in the special case of a diagonalized matrix, the determinant will be simply the product of all the values along the diagonal. So, what do we call it when a matrix has a determinant of zero? According to the definition of an invertible matrix, "A square matrix that is not invertible is called singular or degenerate. A square matrix is singular if and only if its determinant is 0." And that's why you are getting the error at the bottom of the stack trace--at some point, the code needs to invert the covariance matrix for whatever reason (you can see that the final method at the bottom of the trace is called _data_inv_cov()) but it can't do it because the matrix is singular and therefore non-invertible.

Bottom line, what's effectively happening is that, by giving the code repeated instances of the same number to use as input, you are basically generating the linear algebra equivalent of a divide-by-zero error.

Answer 2

Float32 only has 1e-6 precision in numpy, therefore, if you are manipulating small numbers, similar instances could become identical (or very close) therefore producing singular or badly scaled matrices. This issue is particularly tricky to as there are no algebric reason for the desired inversion not to be possible. Two easy way to solve this issue could be to scale you dataset by a power of 10 (i.e. 10e6) or to use float64 or double point precision which have a 1e-15 precision in numpy.

dataset.astype('float64') is an easy way to go about the change.