カーネル密度推定におけるバンド幅はいくつかの選び方があり、silvermanのrule of thumbが有名。

これは単峰な分布向けで、データの散布度をサンプルサイズや何やらで補正した値になる。計算に用いられるデータの散布度はソフトウェアによって異なり、sdだったりmin(sd,IQR/1.34)だったりする。

外れ値がある場合は後者の方がマシな結果を返すのだが、最近使ったscipy.stats.gaussian_kdeは前者だった。そのため外れ値を持つデータにナイーブに適用すると大変なだらかなpdfが得られ、嬉しくない。

scipy.stats.gaussian_kdeはweightsが指定できるのが便利で、こいつのバンド幅をIQRを使う定義に変更したい。ただし1次元の時だけに対応する。2次元以上への対応は面倒*1だし、今は用がない。

というわけで作った。

と言っても公式の実装をコピペして書き換えただけ*2。計算のキーになるweighted quantileはstatsmodels.stats.weightstatsを使っている。道具が整備されてる言語は本当に便利。

from scipy import stats
from scipy import linalg
from numpy import (atleast_2d, pi, cov)
import numpy as np
from statsmodels.stats import weightstats

class gaussian_kde_1d_v2(stats.gaussian_kde):
    def _compute_covariance(self):
        if self.d !=1:
            raise ValueError("1-D array only.")
        self.factor = self.covariance_factor()
        if not hasattr(self, '_data_cho_cov'):
            _d = weightstats.DescrStatsW(self.dataset[0],weights=self.weights)
            _q = _d.quantile([0.75,0.25], return_pandas=False)
            _cv = cov(self.dataset, rowvar=1, bias=False, aweights=self.weights)
            self._data_covariance = atleast_2d(np.min([(_q[0]-_q[1])/1.34, _cv]))
            self._data_cho_cov = linalg.cholesky(self._data_covariance, lower=True)

        self.covariance = self._data_covariance * self.factor**2
        self.cho_cov = (self._data_cho_cov * self.factor).astype(np.float64)
        self.log_det = 2*np.log(np.diag(self.cho_cov * np.sqrt(2*pi))).sum()
    
    @property
    def inv_cov(self):
        self.factor = self.covariance_factor()
        _d = weightstats.DescrStatsW(self.dataset[0],weights=self.weights)
        _q = _d.quantile([0.75,0.25], return_pandas=False)
        _cv = cov(self.dataset, rowvar=1, bias=False, aweights=self.weights)
        self._data_covariance = atleast_2d(np.min([(_q[0]-_q[1])/1.34, _cv]))
        return linalg.inv(self._data_covariance) / self.factor**2

動かす

外れ値を混入させた配列を作る。生成分布はN(0,1)にN(60,9)が少しだけ混ざってる感じ。

n = 70
x =np.random.normal(size=n)
y = np.random.normal(loc=60,size=3)
data = np.concatenate([x,y],0)
w = np.random.randint(1,3,size=n+3)

改良していないkdeのdensityはこれ↓。バンド幅が大きい雰囲気が出ている。

改良した方のkdeのdensityはこれ↓。結構マシになった。

*1:MCD使えばできそう？

*2:covariance_factorを作る部分を差し替えただけ。

べんきのにっき

いろいろと垂れ流します。

バンドかえようぜ

動かす