はじめに

最近，ずっとTheoretical Neuroscienceを読んで計算論的神経科学の勉強を深めています．
インプットばかりでは自分の力にはならないと思うので当然アウトプットが必要なんですが，その方法にも

実際にプログラムを書いて実装してみる
勉強したことをどこかで書き綴る

の2択ぐらいがあるんではないかと．
前者はCRCNSという，実際に研究で使われたデータセットをアーカイブしているウェブサイトがあるのでそこからデータを落としてきていろいろ試すことができます．
後者はせっかくブログを書いているのでここをノートにしようかと思っています．

この分野に精通している方からしたら当然みたいなことを書いているかもしれませんが，まぁ備忘録ということで書かせていただければなぁと思っています．

ただし，参照している書籍が英語で書かれており，かつ私は日本語でこのブログの記事を書いているので，日本語で何て言うのかわからないような単語は英語のまま書いたりするつもりです．

発火率の種類

まず前提として以下では，時間 $T$ にわたってニューロンの発火を観測することを考えます．
さらにその試行は， $N$ 回繰り返されるものとしましょう．

Firing rate

「発火率」を直訳すると"firing rate"になってしまうので，
発火率の一種としてこれを記すのは違和感たっぷりですが，書籍の中では"firing rate"と呼ばれている数値です．

時刻 $t$ の関数として $\mathsf{r}(t)$ と書くのですが， $\mathsf{r}(t)$ は言い換えると，時刻 $t$ における発火の確率密度 $p [t]$ です．
如何せん時間というのは連続なものですから，厳密な $t$ に対して発火率を求めることは当然不可能です( $\mathsf{r}(t) \simeq 0$ )．
なので通常は， $t \sim t + \Delta t$ の間での発火率を求めます．
試行1回に対してこの値を求めようとすると，「発火したか発火してないか」の2値になってしまうので，複数回の試行の平均として求めます．

実験対象のニューロンの，時刻に依存する発火率の変遷などを見ることができます．

Spike-count rate

直訳すると「数え上げ発火率」とかでしょうか？
これは1回の試行についての値です．

時間 $T$ の試行全体での発火率を求めたいので，引数に時刻を取らない定数 $r$ として求められます．
ある試行中で $n$ 回発火があったとすると， $r$ は

${ \displaystyle r = \frac{n}{T} }$

ということになります．

実験対象のニューロンが与えた刺激に対してどのくらい反応したかを，試行ごとに見ることができます．

Trial average rate

これは，「firing rateの時刻についての平均」とも言えますし，「spike-count rateの試行平均」とも言える値です．
前述のふたつの数値，firing rateとspike-count rateを，前者は試行時間全体で，後者は試行回数でそれぞれ割ると，どちらもTrial average rateという数値になります．

「このニューロンにこの刺激を与えたらたいていこのくらい反応する」というざっくりした様子を見ることができます．

発火率の測り方

一応，ここからが本題です．
なぜ「発火率の測り方」なんていうものを考える必要があるかというと，発火率を時間の関数として記述すると，論ずることができる内容もかなり広がるからです． spike-count rateやtrial average rateは時間に依存しないので簡単に求めることができるのに対し， firing rateについては $t$ を厳密に定めて $\mathsf{r}(t)$ を求めることができません．
なので，時間の関数としての発火率を測るためにはいろいろ方法を工夫する必要があるよなぁ，というお話です．

発火率を図るためには，当然，発火のデータが必要です．
試しにプロットしてみたので，こちらを使ってみましょう．

f:id:arthurs_kurt:20190202135245p:plain — スパイクデータ

横軸が時間（秒）に対応していて， 3秒間の観測を想定しています．
縦軸は気にしないでください．このデータは，別に実際のニューロンを選んで観測したものではなく，平均の発火率を70Hzとして確率論的にプロットしたものです．
これについても，もしかしたらそのうち記事を書くかもしれません．

時間をビンにわける方法

最も単純な方法は，試行時間（この場合は3秒）を，いくつかのビン(=Time bins)に分ける方法です．

ここでは，0.1秒のビン×30に分けるとします．
1つのビンの中でスパイクの回数をカウントし，その値をビンのサイズで割ることで，発火率を時間関数として記述します．

f:id:arthurs_kurt:20190202143522p:plain — Time binsによる計測

横軸は先ほどと同じで，時間（秒）．縦軸は，発火率（Hz）です．
スパイクのグラフと見比べてみると，うまく発火率が測れていることがわかると思います．

ただし，これはビン内では一定の発火率となっており，得られる結果はビンのサイズに大きく左右されます．
不連続な関数として記述されるので，真に「時間によって変化する発火率」として表すことのできる指標ではないのかもしれません．

ちなみに，このビンサイズを試行時間と同じにすると，上述のSpike-count Rateになります．

当然ながらSpike-count Rateと比べれば時間解像度はずいぶん高いのですが，後述の方法ならばもっともっと時間解像度の高い $\mathsf{r}(t)$ を求めることができます．

Rectangular window

試行時間をビンにわける方法は，ビンのサイズに依存すると同時に，ビンの位置にも依存しています．

こうしたビンの位置による任意性を排除するための方法としては，一定の大きさ（たとえば前述のビンと同じ0.1秒）のウインドウを用意して，それを時間軸上でスライドさせていく方法があります．
$t$ を中心に置いたウインドウ内で観測されているスパイクの数を，ウインドウのサイズで割ったものを $\mathsf{r}(t)$ とします． 下の図は，ウインドウサイズを $\Delta t = 0.1$ とした場合のRectangular windowによる計測の結果です．

f:id:arthurs_kurt:20190202161009p:plain — Rectangular windowによる計測

計測結果を見ると，ビンの位置による任意性が除かれているのがわかると思います．
今回は，利用しているスパイクデータが，サンプリング周波数＝1000Hzとしているので，かなりギザギザな見た目になっていますが，もう少し高いサンプリング周波数に対してこの方法で発火率を図れば，もう少し見た目も滑らかなグラフになります．

Gaussian window

上述のRectangular windowによる計測でも十分なのですが，もう少し踏み込んでみましょう．
当然のことですが， $t$ における発火率を図るのに， $t$ と $t+\Delta t$ の発火を同等に扱うのは不自然な気もします．
ですので， $t$ での発火には大きめの重みをかけ，周辺に離れるに従って徐々にその重みを小さくしていくことを考えてみると，ガウス分布によるウインドウが尤もらしく見えます．
このようなウインドウをGaussian windowと呼びます．

Gaussian windowは，
${\displaystyle w(\tau) = \frac{1}{\sqrt{2\pi} \sigma_w} \exp( -\frac{\tau ^2}{2 \sigma ^2_w} ) }$
という式で定義されます．

このウインドウ関数を用いて
${\mathsf r}(t) = \int^{\infty}_{-\infty} w(\tau) \rho (t-\tau) d\tau$
として発火率を求めます．
$\rho(z)$ は， $z$ の瞬間にスパイクが発生していれば $1$ を，それ以外では $0$ を返します．
中身はディラックのデルタ関数とかを使って定義されています．

この様なインテグラルをLinear filterと呼びます．
ちなみに，Rectangular windowも， $r(\tau)$ に異なる関数を選択しているだけで，やっている処理はこれと同じで，線形フィルタリングです．

f:id:arthurs_kurt:20190202183345p:plain — Gaussian windowによる計測

Gaussian windowによって ${\mathsf r}(t)$ を計測すると，このようになります．
使っているウインドウ関数が滑らかなのも相まって，グラフ自体もかなり滑らかになったのがわかると思います．

Alpha window

ガウスウインドウでは，遠くに離れるほど重みを減らそうというアイデアがありました．

さらにもう少し突き詰めたのが，Alpha windowです．
どういう風に突き詰めるかというと，離れるほど重みを減らすのは同じなのですが， ある時点での発火率には，その時点より前の発火しか加味しない というアイデアです．

ウインドウ関数には
$w(\tau) = \mathsf{ReLU}(\alpha ^2 \tau \exp (- \alpha \tau))$
を用います．
ニューラルネットワークの活性化関数としてもおなじみ，ReLU（＝ Rectified Linear Unit）が出てきましたね．
ReLUは $\mathsf{ReLU}(z) = \max (0, z)$ と書くこともできて，負の値はすべて0にしてくれる関数です．

実行した結果はこのようになります．

f:id:arthurs_kurt:20190203103157p:plain — Alpha windowによる計測

概形はGaussian windowによるものと似ていますね．
特徴としては，それより前の点でスパイクが観測されていない $t = 0$ 付近の点で $\mathsf r (t) = 0$ になっていることです．
また，滑らかさはGaussianには劣るものの，そちらと比べると値の変動が激しく，より細かな特徴をとらえているように見えます．