HOME PAGE

Spike Density Estimation

- Histogram Method

- Kernel Method

Share it with your friend

Tweet

補遺 統計的推測

統計的推測の理論と実践はカール・ピアソン(Karl Pearson)とロナルド・エイルマー・フィッシャー(R.A. Fisher)によって築かれた. ピアソンはデータの背後に潜む統計モデルという概念を明確にし, 統計モデルの母数をデータから決定して科学的な検証に使うことを推進した. フィッシャーはデータから得られる統計量と統計モデルの母数とを明確に区別し，統計量（推定量）の望ましい性質を論じた．そして，繰り返し行われる実験によりデータを集めて, 統計モデルの母数を推定する最尤法を推進し統計的推測の理論を確立した.     フィッシャーの打ち立てた統計的推測の理論の強力な仮定は, 母数が定数であることである. 計測期間中に母数が変動したり, 試行ごとに母数が変動することは想定しない. 神経科学における実験パラダイムも, これに従い同一実験下で 繰り返し実験を行うのが主流である. 現代はこのような制約から離れ，現実世界で得られるデータから統計的な判断を下す必要に迫られているために, ベイズ統計学が隆盛である.

平均二乗誤差

$\hat{\theta}$をパラメータ$\theta$の推定量とする. 例をあげてみよう. $X_{i}$が平均$\theta$の指数分布に従う乱数とすると, $\theta$の推定量 $\hat{\theta}%$の例としてまず考えられるのが, 平均 $\hat {\theta}_{n}=\frac{1}{n}\sum X_{i}$, がある. 推定量は他にも考えられる. 例えば初めのサンプルだけを使う $\hat{\theta}_{1}=X_{i}$も$\theta$の推定量である. では, 望ましい推定量とはどのような性質を持つべきだろうか. 平均二乗誤差はその一つの基準を与えると考えられる．

$$=E[(\hat{\theta}-\theta)^{2}]%$$ (4.1)

二乗損失関数、損失関数を真の分布で期待値をとった関数をリスク関数という. 狭義の意味では, リスク関数の最小化が統計学や学習理論のテーマである. 推定量の期待値を $E\hat{\theta}=\bar{\theta}$とすると,

$$=E[(\hat{\theta}-\bar{\theta}+\bar{\theta}-\theta )^{2}]$$ MSE

$$=E[(\hat{\theta}-\bar{\theta})^{2}+2(\hat{\theta}-\bar{\theta})(\bar {\theta}-\theta)+(\bar{\theta}-\theta)^{2}]$$

$$=E[(\hat{\theta}-\bar{\theta})^{2}]+(\bar{\theta}-\theta)^{2} %$$ (4.2)

第一項 $E[(\hat{\theta}-\bar{\theta})^{2}]$は推定分散. 第二項はバイアスの二乗である. 従ってバイアスがなく, 分散が小さい推定量が望ましい推定量である. 隔たりのある推定量が小さい分散を持ち, そのMSEが可能な不偏推定量のMSEよりも小さくなることはありうる. これが望ましい推定量であるのは間違いないが, 歴史的には不偏性, 有効性は別々に吟味されてきた. まず普遍性のある推定量を見つけて, 次に不偏推定量の中で分散最小のものを選ぶという手順を踏む. 古典的な推定量についてその性質を見てみよう.

問題 13   Gauss分布に従う確率変数に対する統計量 $S_{n-1}% ,S_{n},S_{n+1}$のMSEを求めよ．

ただし $S_{n-1}^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(X_{i}-\bar {X})^{2}$, $S_{n}^{2}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}$, $S_{n+1}^{2}=\frac{1}{n+1}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}$ .

どれが一番小さいか．

解法 14   MSE $(S_{n+1})<$MSE $(S_{n})\leq$MSE$(S_{n-1})$

$$(S_{n-1}) =E(S_{n-1}-\sigma^{2})^{2}=\frac{2}{n-1}\sigma^{4}$$

$$(S_{n}) =E(S_{n}-\sigma^{2})^{2}=\frac{2n-1}{n^{2}}\sigma^{4}$$

$$(S_{n+1}) =E(S_{n+1}-\sigma^{2})^{2}=\frac{2}{n+1}\sigma^{2}%$$

最尤推定

尤度(likelihood)を最大にするパラメータを最尤推定量という．

$$l(\theta)=\log p\left( X\vert\theta\right) =\sum\log p\left( X_{i}% \vert\theta\right)$$

例    ポアソン分布の場合

$$l(\theta) =\sum_{i=1}^{n}\log p\left( X_{i}\vert\theta,\Delta\right)$$

$$=\sum_{i=1}^{n}\log\frac{\left( \theta\Delta\right) ^{X_{i}}}{X_{i}% !}e^{-\theta\Delta}$$

$$=\sum_{i=1}^{n}[X_{i}\log\left( \theta\Delta\right) -\theta\Delta-\log X_{i}!]$$

最尤推定は

$$\frac{\partial l(\theta)}{\partial\theta} =\sum_{i=1}^{n}[\frac{X_{i}% }{\theta}-\Delta]$$

$$=\frac{\sum_{i=1}^{n}X_{i}}{\theta}-n\Delta$$

$\frac{\partial l(\theta)}{\partial\theta}=0$より

$$\hat{\theta}_{n}=\frac{\sum_{i=1}^{n}X_{i}}{n\Delta}=\frac{\bar{X}}{\Delta}%$$

ただし平均値 $\bar{X}=\frac{1}{n}\sum X_{i}$

最尤推定値の期待値は

$$E\hat{\theta}_{n}=\frac{\sum_{i=1}^{n}EX_{i}}{n\Delta}=\theta$$

分散は

$$E(\hat{\theta}_{n}-E\hat{\theta}_{n})^{2} =E\left[ \frac{\sum_{i=1}% ^{n}X_{i}}{n\Delta}-\theta\right] ^{2}$$

$$=\frac{1}{(n\Delta)^{2}}\sum_{i=1}^{n}E(X_{i}-\Delta\theta)^{2}$$

$$=\frac{\theta}{n\Delta}%$$

最後の等式にはポアソン分布の性質 $E(X_{i}-\Delta\theta)^{2}=EX_{i}=\Delta\theta$を用いた.

十分統計量

発火率$\lambda$の定常ポアソン過程に従うスパイク時系列を得たとしよう. スパイク時系列をごく短い$\Delta$秒の$n$個のビンに区切り, $t$番目のビンにスパイクが存在すれば$X_{t}=0$, しなければ$X_{t}=1$として, $X=\left( X_{1},X_{2},\ldots,X_{n}\right)$をベルヌーイ過程と見なす. このとき統計量 $T=\sum_{t=1}^{n}% X_{i}$はベルヌーイ過程の$\lambda$に関する十分統計量である.

$$P\left( X\vert\lambda,\Delta\right) =% {\displaystyle\prod\limits_{t=1}^{n}} P\left( X_{i}\vert\lambda,\Delta\right)$$

$$=% {\displaystyle\prod\limits_{t=1}^{n}} (\lambda\Delta)^{X_{i}}\left( 1-\lambda\Delta\right) ^{1-X_{i}}$$

$$=(\lambda\Delta)^{X_{1}+X_{2}+\cdots+X_{n}}(1-\lambda\Delta)^{n-(X_{1}% +X_{2}+\cdots+X_{n})}%$$

$\Delta$は決まっているとすれば, $\lambda$の値を推定するのに, ここのデータ $\left( X_{1},X_{2},\ldots,X_{n}\right)$は必要でなく, 統計量 $T=\sum_{t=1}^{n}X_{i}$を用いれば十分であることがわかる.

今度は定常ポアソン過程をやや大きめの$n$個の観測区間$\Delta$で区切ったとしよう. このとき各区間内のスパイクの個数 $X=\left( X_{1},X_{2},\ldots,X_{n}\right)$は平均 $\lambda\Delta$のポアソン分布に従う. このとき $T=\sum_{t=1}^{n}X_{i}%$は$\lambda$に関する十分統計量である.

$$P\left( X\vert\lambda,\Delta\right) =% {\displaystyle\prod\limits_{i=1}^{n}} P\left( X_{i}\vert\lambda,\Delta\right)$$

$$=% {\displaystyle\prod\limits_{t=1}^{n}} \frac{\left( \lambda\Delta\right) ^{X_{i}}}{X_{i}!}e^{-\lambda\Delta}$$

$$=\frac{1}{X_{1}!X_{1}!\cdots X_{n}!}\cdot\left( \lambda\Delta\right) ^{X_{1}+X_{2}+\cdots+X_{n}}e^{-n\lambda\Delta}%$$

最尤推定を考えてみよう. 因子 $\frac{1}{X_{1}!X_{1}!\cdots X_{n}% !}$は$\lambda$に関係がない. $\lambda$に関係する第二因子は $\left( \lambda\Delta\right) ^{T}e^{-n\lambda\Delta }$と$T$の関数として書けるから, 個々のデータを知る必要はないことがわかる.

一般に密度分布のモデル $f\left( x\vert\theta\right)$が

$$f\left( x\vert\theta\right) =h\left( x\right) g\left( T\left( x\right) \vert\theta\right)$$ (4.3)

の形に因子分解できる場合, かつその場合に限り, 統計量$T$はモデルパラメータ$\theta$の十分統計量である.    これをフィッシャー・ネイマンの因子分解定理という. 密度分布のモデル $f\left( x\vert\theta\right)$を$\theta$に依存する因子と依存しない因子に分解できて, 依存する因子が統計量 $T\left( x\right)$の関数で表すことができるとすれば, $\theta$の推定には $T\left( x\right)$を用いれば十分である. ベルヌーイ過程の例では $h\left( x\right) =1$, $g\left( T\left( x\right) \vert\theta\right) =\theta^{T}(1-\theta)^{n-T}$と見て, $T=\sum_{t=1}^{n}X_{i}%$が十分統計量であることがわかる. またポアソン分布は $h\left( x\right) =\frac{1}{X_{1}!X_{1}!\cdots X_{n}!}$, $g\left( T\left( x\right) \vert\theta\right) =\left( \lambda\Delta\right) ^{T}e^{-n\lambda\Delta}%$となり, $T=\sum_{t=1}^{n}X_{i}$が十分統計量であることがわかる.

定義 15   データ $X=\left( X_{1},X_{2},\ldots,X_{n}\right)$が $T\left( X\right)$が与えられた下で$\theta$と独立な時, 関数 $T\left( X\right)$をモデル $f\left( x\vert\theta\right)$の十分統計量という.

$$P\left( X=x\vert T\left( X\right) =t,\theta\right) =P\left( X=x\vert T\left( X\right) =t\right)$$

フィッシャー情報量

スコアと呼ばれる統計量を導入する

$$V\left( X;\theta\right) =\frac{\partial}{\partial\theta}\log f\le... ...rt\theta\right) }\frac{\partial f\left( X\vert\theta\right) }{\partial\theta}%$$ (4.4)

スコアの期待値は

$$E[V] =\int V\left( x;\theta\right) f\left( X\vert\theta\right) dx$$

$$=\int\frac{\partial f\left( x\vert\theta\right) }{\partial\theta}dx$$

$$=\frac{\partial}{\partial\theta}\int f\left( x\vert\theta\right) dx$$

$$=0$$

スコアの分散は

$$E[V^{2}]=E\left[ \frac{\partial}{\partial\theta}\log f\left( X\vert\theta \right) \right] ^{2}%$$

スコアの分散をフィッシャー情報量とよび, $J\left( \theta\right) =E[V^{2}]$で表す.

コーシー・シュワルツの不等式(Cauchy-Schwarz inequality)

$$\int\left\vert f\left( x\right) \right\vert \,^{2}dx\cdot\int\lef... ...}dx\geq\left\vert \int f\left( x\right) g\left( x\right) dx\right\vert \,^{2}%$$ (4.5)

を用いれば次の不等式が得られる.

$$E[V^{2}]\cdot E[(\hat{\theta}-\bar{\theta})^{2}]\geq\left\vert E\left[ V(\hat{\theta}-\bar{\theta})\right] \right\vert ^{2}%$$

ここで左辺の第一因子はフィッシャー情報量 $E[V^{2}]=J\left( \theta\right)$. 第二因子は推定量の分散, $var\hat{\theta}=E[(\hat{\theta}-\bar{\theta})^{2}% ]$である. 右辺についてはは二乗の中身が

$$EV\hat{\theta} =% {\displaystyle\int} V(\hat{\theta}-\bar{\theta})f\left( x\right) dx$$

$$=\int\frac{1}{f\left( x\vert\theta\right) }\frac{\partial f\left(... ...ght) }{\partial\theta}(\hat{\theta}-\bar{\theta})f\left( x\vert\theta\right) dx$$

$$=\int(\hat{\theta}-\bar{\theta})\frac{\partial f\left( x\vert\theta\right) }{\partial\theta}dx$$

$$=\int\hat{\theta}\frac{\partial f\left( x\vert\theta\right) }{\pa... ...-\bar{\theta}\frac{\partial}{\partial\theta}\int f\left( x\vert\theta\right) dx$$

$$=\frac{\partial}{\partial\theta}\int\hat{\theta}f\left( x\vert\theta\right) dx$$

$$=\frac{\partial}{\partial\theta}E\hat{\theta}%$$

となる. ゆえに統計量の分散の下限はフィッシャー情報量の逆数で与えられる.

$$var\hat{\theta}\geq\frac{\left( \frac{\partial E\hat{\theta}}{\partial\theta }\right) ^{2}}{J\left( \theta\right) }%$$ (4.6)

この不等式をクラマー・ラオ不等式という. フィッシャー情報量は統計量をある特定のモデル $f\left( x\vert\theta\right)$のパラメータ$\theta$の推定値に用いる場合の分散の下限を与える.    不偏推定量 $E\hat{\theta}=\theta$の場合

$$var\hat{\theta}\geq\frac{1}{J\left( \theta\right) }%$$ (4.7)

バイアスがある場合は $E\hat{\theta}=\theta+b\left( \theta\right)$

$$var\hat{\theta}\geq\frac{1+b^{\prime}\left( \theta\right) }{J\left( \theta\right) }%$$ (4.8)

フィッシャー情報量は次のようにも書ける.

$$J\left( \theta\right) =E\left[ \left( \frac{\partial}{\partial\th... ...rac{\partial^{2}% }{\partial\theta^{2}}\log f\left( x\vert\theta\right) \right]$$ (4.9)

問題 16   フィッシャー情報量の等式を証明しなさい.

解法 17   右辺の中身について,

$$-\frac{\partial^{2}}{\partial\theta^{2}}\log f\left( x;\theta\rig... ...}}}{f^{2}}=-\frac{f^{^{\prime\prime}}}% {f}+\frac{(f^{^{\prime}})^{2}}{f^{2}}%$$

の関係があるから, 期待値をとって

$$E\left[ -\frac{f^{^{\prime\prime}}}{f}+\left( \frac{f^{^{\prime}}... ...e}}}{f}\right] ^{2}=E\left[ \frac{\partial}{\partial\theta}\log f\right] ^{2}%$$

ポアソン分布のフィッシャー情報量. 定常ポアソンスパイク時系列を$n$個の観測区間$\Delta$で区切ったとしよう. このとき各区間内のスパイクの個数 $X=\left( X_{1},X_{2}% ,\ldots,X_{n}\right)$は平均 $\lambda\Delta$のポアソン分布

$$\Pr\left( X=\left( X_{1},X_{2},\ldots,X_{n}\right) \vert\lambda,\... ...1}^{n}} \frac{\left( \lambda\Delta\right) ^{X_{i}}}{X_{i}!}e^{-\lambda\Delta}%$$

に従う. フィッシャー情報量は

$$J\left( \lambda\right) =-E\left[ \frac{\partial^{2}}{\partial \lambda^{2}}\log f\left( X;\lambda,\Delta\right) \right]$$

$$=-E\left[ \frac{\partial^{2}}{\partial\lambda^{2}}\sum_{i=1}^{n}\left( X_{i}\log\lambda\Delta-\lambda\Delta-\log X_{i}!\right) \right]$$

$$=-E\left[ \sum_{i=1}^{n}\frac{\partial}{\partial\lambda}\left( X_{i}% \frac{1}{\lambda}-\Delta\right) \right]$$

$$=\frac{\sum_{i=1}^{n}EX_{i}}{\lambda^{2}}$$

$$=\frac{n\lambda\Delta}{\lambda^{2}}=\frac{n\Delta}{\lambda}%$$

従って分散の下限は $\lambda/n\Delta$.

問題 18   ポアソン分布 $p\left( X\vert\lambda ,\Delta\right) =% {\displaystyle\prod\limits_{t=1}^{n}} \frac{\left( \lambda\Delta\right) ^{X_{i}}}{X_{i}!}e^{-\lambda\Delta}%$の$\theta$に関するフィッシャー情報量を公式 $J\left( \theta\right) =E\left[ \left( \frac{\partial }{\partial\theta}\log f\left( x\vert\theta\right) \right) ^{2}\right]$を用いて導出せよ.

問題 19  

指数分布のフィッシャー情報量を求めよ$\quad$

解法 20   平均$1/\lambda$の指数分布 $f\left( x;\lambda\right) =\lambda e^{-\lambda x}$の$\lambda$に関するフィッシャー情報量は

$$J\left( \lambda\right) =-E\left[ \frac{\partial^{2}}{\partial\lam... ...tial\lambda}\left( \frac{1}{\lambda}-x\right) \right] =\frac {1}{\lambda^{2}}%$$

ゆえに分散の下限は $\lambda^{2}$．

ベイズ推定（はやわかり）

Bayesの定理

$$p\left( \theta\vert x,w\right) =\frac{p(x\vert\theta)\pi(\theta\vert w)}{p(x\vert w)}%$$

$p(x\vert\theta)$を尤度(likelihood), $\pi(\theta)$を事前分布(prior distribution)， $p\left( \theta\vert x\right)$ を事後分布(posterior distribution)という．

分母は正規化項で分子の総和をとって

$$p(x\vert w)=\int p\left( x,\theta\vert w\right) d\theta$$

で与えられる．周辺尤度関数・エビデンスという．

事後分布による推定値

事後損失関数(posterior loss function, posterior expected loss) を定義し，事後分布による期待値: 事後リスク関数(posterior risk function)    

$$L^{\text{Bayes}}(\delta,\theta;X,w)\equiv E^{\theta\vert X,w}[L(\... ...eta)]=\int L\left( \delta,\theta\right) p\left( \theta\vert X,w\right) d\theta$$

を最小化するパラメータを推定量とする

$$\hat{\theta}_{\text{Bayes}}(X,w)=\arg\underset{\theta}{\min}L^{\text{Bayes}% }(\delta,\theta;X,w)$$

損失関数として支持関数(indicator function)$L_{0}$

$$L_{0}(\delta,\theta)=1_{\delta,\theta}%$$

を選ぶとMAP推定値(maximum a posteriori estimate)が得られる．

$$\hat{\theta}_{\text{MAP}}(X,w)=\arg\underset{\theta}{\max}p\left( \theta\vert X,w\right)$$

損失関数として絶対誤差$L_{1}$

$$L_{1}(\delta,\theta)=\vert\delta-\theta\vert$$

を選ぶと事後中央値推定(posterior median)が得られる．

損失関数として二乗誤差$L_{2}$

$$L_{2}(\delta,\theta)=(\delta-\theta)^{2}%$$

を用いると事後平均(posterior mean)

$$\hat{\theta}_{\text{PM}}(X,w)=\int\theta p\left( \theta\vert X,w\right) d\theta$$

が得られる．事後平均を特にベイズ推定値(Bayes estimate)と呼ぶことがある．

問題 21   事後損失 $L^{Post}(\delta,\theta)$最小化の原理からMAP推定値, 事後平均値を求めよ.

経験ベイズ推定量

ベイズ推定量は超パラメータに依存する．超パラメータを最尤推定で求めベイズ推定量に用いることを経験ベイズ法と呼ぶ(empirical Bayes estimator)．

周辺尤度

$$p(x\vert w)=\int p\left( x,\theta\vert w\right) d\theta=\int p(x\vert\theta)\pi (\theta\vert w)d\theta$$

周辺尤度の最尤推定値(周辺尤度最大化)

$$\hat{w}=\arg\underset{w}{\max}p(x\vert w)$$

経験ベイズ推定量は最適化した超パラメータ$\hat{w}$を用いて

$$\hat{\theta}_{\text{EB}}(X)=\hat{\theta}_{\text{Bayes}}(X,\hat{w})$$

で与えられる．