Подгонка гауссовского KDE в numpy / scipy в Python

Question

Я подгоняю гауссовский оценщик плотности ядра к переменной, которая является разностью двух векторов, называемых "diff", следующим образом: gaussian_kde_covfact (diff, smoothing_param) - где gaussian_kde_covfact определяется как:

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()

Это работает, но есть граничный случай, когда diff является вектором всех 0. В этом случае я получаю ошибку:

 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

Какой способ обойти это? В этом случае я бы хотел, чтобы он возвращал плотность, которая по существу полностью достигла пика при разнице 0, при этом нигде больше нет массы.

спасибо.

Answer 1

Пиковая плотность, масса которой в одной точке не является гауссовой, поэтому, строго говоря, то, что вы хотите сделать, не определено (и такое распределение не имеет конечной ковариации).

Теперь, в вашем случае, для вектора, который равен нулю, вы можете использовать его в особом случае, минуя всю инфраструктуру. Простой способ обнаружить случай - это вычислить максимальное значение разности и сравнить его с eps (numpy.finfo (x.dtype) .eps для вектора x). Вы также можете просто обнаружить его, перехватив LinalgError, но вы должны быть осторожны, чтобы различать случаи, когда ковариация плохо определена и 0 записей.