Installing Gaussian KDE in numpy / scipy in Python
I map the density estimate of a Gaussian kernel to a variable that is the difference of two vectors called "diff" as follows: gaussian_kde_covfact (diff, smoothing_param) - where gaussian_kde_covfact is defined as:
class gaussian_kde_covfact(stats.gaussian_kde):
def __init__(self, dataset, covfact = 'scotts'):
self.covfact = covfact
scipy.stats.gaussian_kde.__init__(self, dataset)
def _compute_covariance_(self):
'''not used'''
self.inv_cov = np.linalg.inv(self.covariance)
self._norm_factor = sqrt(np.linalg.det(2*np.pi*self.covariance)) * self.n
def covariance_factor(self):
if self.covfact in ['sc', 'scotts']:
return self.scotts_factor()
if self.covfact in ['si', 'silverman']:
return self.silverman_factor()
elif self.covfact:
return float(self.covfact)
else:
raise ValueError, \
'covariance factor has to be scotts, silverman or a number'
def reset_covfact(self, covfact):
self.covfact = covfact
self.covariance_factor()
self._compute_covariance()
This works, but there is an edge case where diff is a vector of all 0s. In this case, I get the error:
File "/srv/pkg/python/python-packages/python26/scipy/scipy-0.7.1/lib/python2.6/site-packages/scipy/stats/kde.py", line 334, in _compute_covariance
self.inv_cov = linalg.inv(self.covariance)
File "/srv/pkg/python/python-packages/python26/scipy/scipy-0.7.1/lib/python2.6/site-packages/scipy/linalg/basic.py", line 382, in inv
if info>0: raise LinAlgError, "singular matrix"
numpy.linalg.linalg.LinAlgError: singular matrix
What's the way to get around this? In this case, I would like it to return a density that has essentially peaked with a difference of 0, with no mass all over the place.
thanks.
a source to share
The density that peaks at a single point is not Gaussian, so strictly speaking, what you want to do is undefined (and such a distribution has no finite covariance).
Now, in your case, for a vector that is zero, you can traverse the entire infrastructure in a special case. An easy way to detect a case is to calculate the maximum diff value and compare that to eps (numpy.finfo (x.dtype) .eps for vector x). You can also just catch it by catching a LinalgError, but you need to be careful to distinguish between cases where covariance is poorly defined and 0 entries.
a source to share