2020-01-03

ベイズ統計の理論と方法：ノート5.2節その1（平均場近似）

本読み

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらコメント等でご指摘いただけますと幸いです。

ベイズ統計の理論と方法

作者:渡辺澄夫
出版社/メーカー: コロナ社
発売日: 2012/03/01
メディア: 単行本

書籍サポートページ

ベイズ統計の理論と方法：ノート1：主要な量の定義。
ベイズ統計の理論と方法：ノート2：指数型分布の場合の具体的な事後分布とか。
ベイズ統計の理論と方法：ノート2.5：コイン投げの例とか。
ベイズ統計の理論と方法：ノート3： 1章の章末問題と、真の分布と確率モデルの関係。
ベイズ統計の理論と方法：ノート4章その0.1：開集合とか。
ベイズ統計の理論と方法：ノート4：キュムラント母関数。
ベイズ統計の理論と方法：ノート4章その1：多様体の導入だけ。
ベイズ統計の理論と方法： 4章の章末問題【2】（実対数閾値とその多重度）

f:id:cookie-box:20200101101603p:plain:w60

…データから混合指数型分布 $p(x|w) = \sum_{k=1}^K a_k v(x) \exp \bigr( f (b_k) \cdot g(x) \bigl)$ を推定するのに、

「 $k$ の目が出る確率が $a_k v(x) \exp \bigr( f (b_k) \cdot g(x) \bigl)$ であるようなサイコロ」

を仮想的に考えて、

「 $n$ 個のデータ $x^n$ を受け取ったとき、パラメータ $a_{1:K}, b_{1:K}$ は何であったろうか」という問題を
「 $n$ 個のデータ $x^n$ を受け取ったとき、パラメータ $a_{1:K}, b_{1:K}$ と各データのサイコロの目 $y^n$ は何であったろうか」という問題に

読み換えたのですよね？そうすることで事後分布 $P(y^n, a_{1:K}, b_{1:K} | x^n)$ を、独立な分布の積である試行分布 $q(y^n) r(a_{1:K}, b_{1:K})$ で平均場近似することが実現できるというのはわかります。わかりますが、そもそもそんなサイコロを考えなかったらどうなるんです？

f:id:cookie-box:20200101101614p:plain:w60

…まずサイコロを導入した場合はどうなったの？

えっと、そもそも、 $p(w_1, w_2)= \exp \bigl( -\beta H(w_1, w_2) \bigr) / Z(\beta)$ という形で表される目標分布を、独立な分布の積である試行分布 $q(w_1, w_2) = q_1(w_1)q_2(w_2)$ で（カルバック・ライブラー情報量 $D(q || p)$ の意味で）最もよく近似したとき、 $q_1(w_1)q_2(w_2)$ を $p(w_1, w_2)$ の平均場近似とよぶのですね。

それではどのような $q_1(w_1)q_2(w_2)$ が $p(w_1, w_2)$ の平均場近似なのか考えると、これらのカルバック・ライブラー情報量は具体的に、

$D(q || p) = \int q_1(w_1)q_2(w_2) \log \bigl( q_1(w_1)q_2(w_2)/p(w_1, w_2) \bigr) dw$
$\quad = \int q_1(w_1)q_2(w_2) \log q_1(w_1)q_2(w_2) dw - \int q_1(w_1)q_2(w_2) \log p(w_1, w_2) dw$
$\quad = \int q_1(w_1)q_2(w_2) \log q_1(w_1)q_2(w_2) dw + \int q_1(w_1)q_2(w_2) \bigl( \beta H(w_1, w_2) + \log Z(\beta) \bigr) dw$
$\quad = \int q_1(w_1)q_2(w_2) \log q_1(w_1)q_2(w_2) dw + \beta \int q_1(w_1)q_2(w_2) H(w_1, w_2) dw - \beta F(\beta)$
$\quad = \int q_1(w_1) \log q_1(w_1) dw_1 + \int q_2(w_2) \log q_2(w_2) dw_2 + \beta \int q_1(w_1)q_2(w_2) H(w_1, w_2) dw - \beta F(\beta)$
$\quad \geqq 0$

となりますから（ただし、 $F(\beta) = - \beta^{-1} \log Z(\beta)$ は自由エネルギーです）、これを最小化する $q_1(w_1)q_2(w_2)$ が $p(w_1, w_2)$ の平均場近似ということになります。等号成立は $p(w_1, w_2) = q_1(w_1)q_2(w_2)$ のときですが、いま、試行確率分布は成分どうしが独立という制約がありますし、この等号が達成できる保証はありません。ところで、上式の下から2行目の、試行分布に依存しない $- \beta F(\beta)$ を除いた部分を、平均場自由エネルギー $\beta \hat{F}(\beta)$ とよぶとのことです。もし試行分布が目標分布を実現できれば真の自由エネルギーに等しくなる量ですが、一般的には真の自由エネルギー以上になる量ですね。最小化すべき目的関数の別表現といった感じがします。それで、もし $\beta \hat{F}(\beta)$ が最小値になっているならば、この最小化問題のラグランジュ関数

$\int q_1(w_1) \log q_1(w_1) dw_1 + \int q_2(w_2) \log q_2(w_2) dw_2 + \beta \int q_1(w_1)q_2(w_2) H(w_1, w_2) dw_1 dw_2$
$\quad + \lambda_1 ( \int q_1(w_1)dw_1 - 1 ) + \lambda_2 ( \int q(w_2)dw_2 - 1 )$

の $q_1, \, q_2$ に関する1次の変化分が恒等的に $0$ になっているはずです。このことから、

$\quad q_1(w_1) \propto \exp \bigl( - \beta \int q_2(w_2) H(w_1, w_2) dw_2 \bigr) \propto \exp \bigl(\int q_2(w_2) \log p(w_1, w_2) dw_2 \bigr)$
$\quad q_2(w_2) \propto \exp \bigl( - \beta \int q_2(w_1) H(w_1, w_2) dw_1 \bigr) \propto \exp \bigl( \int q_2(w_1) \log p(w_1, w_2) dw_1 \bigr)$

が要請され、これを自己無矛盾条件とよぶと。それで自己無矛盾条件を満たす試行分布をみつけよということですが…一般にラグランジュ関数の停留点は最小点である保証はないのであくまで最小点であることの必要条件になると思うんですが、この場合も極大点であることはないんでしょうか。そもそも必ず最小点があるんでしょうか。テキストの雰囲気的には「自己無矛盾条件を満たす試行分布の中には必ず平均場近似はある」といった感じですが。

あと、平均場自由エネルギー $\beta \hat{F}(\beta)$ も何がエネルギーなのかわかりませんね。いうほど通常の自由エネルギーも何のことかわかりませんが。

自由エネルギー $\beta F(\beta) = -\log Z(\beta)$ って、確率分布が $p(w_1, w_2)= \exp \bigl( -\beta H(w_1, w_2) \bigr) / Z(\beta)$ の形でかかれていることを前提としていて、このとき点 $(w_1, w_2)$ の選択情報量（本当は $p(w_1, w_2)$ は確率じゃなくて確率密度だから選択情報量とはいわないけど）は $\log \bigl( 1/p(w_1, w_2) \bigr)$ だけど、「もし正規化を忘れて選択情報量を $\log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right)$ と考えてしまうと $\beta F(\beta)$ だけ底上げされているよ！ $\beta F(\beta)$ だけ差し引く必要があるよ！」という量なのかな。 $w_1, w_2$ によらずいつも $\log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right)$ から $\beta F(\beta)$ が差し引かれているよね。

$\displaystyle \log \left( \frac{1}{p(w_1, w_2)} \right) = \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) + \log Z(\beta)= \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) - \beta F(\beta)$

選択情報量を「その出来事を知らされたときの驚きの量」みたいに捉えるなら、「必ず起きることを知らされたときの驚きの量をゼロにするために差し引くべきオフセット」が自由エネルギーだといえるかもしれない。じゃあ確率分布 $q_1(w_1)q_2(w_2)$ によるエンコードで同じ点 $(w_1, w_2)$ の選択情報量を測ったものはというと、 $\log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right)$ に比べて $\log \left( 1/\exp \bigl( -\beta H(w_1, w_2) \bigr) \right) + \log q_1(w_1)q_2(w_2)$ だけ差し引かれたものになっているよね。トートロジカルだけど。

$\displaystyle \log \left( \frac{1}{q_1(w_1)q_2(w_2)} \right) = \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) - \log \left( \frac{1}{\exp \bigl( -\beta H(w_1, w_2) \bigr) } \right) - \log q_1(w_1)q_2(w_2)$

この差し引くオフセットは $w_1, w_2$ に依存しているけど、確率分布 $q_1(w_1)q_2(w_2)$ で平均すると平均場自由エネルギー $\beta \hat{F}(\beta)$ になるね。だからやっぱりこれも「驚きの量から差し引くべきオフセットの平均」になっている。このときの驚きの量のエンコーディングは最適ではないかもしれないけど。
$\beta \hat{F}(\beta) = \int q_1(w_1)q_2(w_2) \log \Bigl( q_1(w_1)q_2(w_2)/ \exp \bigl( -\beta H(w_1, w_2) \bigr) \Bigr) dw_1 dw_2$
$\quad \displaystyle =\int q_1(w_1)q_2(w_2) \left( \log \left( \frac{ 1}{ \exp \bigl( -\beta H(w_1, w_2) \bigr)} \right) + \log q_1(w_1)q_2(w_2) \right) dw_1 dw_2$

ただ「平均場自由エネルギー」を「驚きの量から差し引くべき平均的なオフセット」と捉えても、「なぜ平均場自由エネルギーを最小化しなければならないのか」に示唆は得られないな…。いまは「カルバック・ライブラー情報量」つまり「真の驚きの量からの平均的なムダ」を最小化したいから、それと定数 $\beta F(\beta)$ ずれているだけの「（仮の）驚きの量から差し引くべきオフセットの平均」を最小化したいはずだ。エンコードがずれているほど差し引くべきオフセットは平均的にかさむからね。

だからその説明で誰かに伝わりますかね…？まあそれで冒頭の推定をするんですが、つまり、パラメータ $a_{1:K}, \, b_{1:K}$ の分布 $\varphi(a_{1:K}, b_{1:K})$ を更新していくわけです。

パラメータ $a_{1:K}$ の分布はディリクレ分布 ${\rm Dir}(a_{1:K}|\phi)$ ですね。ハイパーパラメータ $\phi \in \mathbb{R}_{+}^K$ を更新していくことになります。ディリクレ分布は $K$ 面あるサイコロの各目が出る割合を生成する確率分布ですね。ある割合における確率密度は「各目が $\phi_k - 1$ 回ずつ出たときにその割合であった確率密度」といえます。
パラメータ $b_{1:K}$ の分布は $\prod_{k=1}^K \bigl(1 / z(\eta_k) \bigr) \exp \bigl( \eta_k \cdot f(b_k) \bigr)$ です。12ページでやった指数型分布の共役な事前分布ですね。ハイパーパラメータ $\eta \in \mathbb{R}^{J \times K}$ を更新していくことになると思います。

12ページの指数型分布のときにはもう事後分布を直接計算しにいきましたが、ここでは謎の確率変数 $Y$ を導入します。これは実質的には $1, \cdots, K$ のいずれかの値を取る確率変数で、元々の確率変数 $X$ と $Y$ の同時分布が混合分布の $k$ 番目の山のみになるようなものです。言い換えると、 $Y$ は $K$ 面サイコロではあるのですが、 $\mathbb{R}^N$ 空間内の点によってサイコロの各面の出る割合が変動します。その割合は、その点における混合分布の $1, \cdots, K$ 番目の山の密度の比になります。 $a_{1:K}, \, b_{1:K}$ から山たちが決まり、山たちから「各点における密度」「各点における各山の密度の比」が決まるので、 $x^n, y^n, a_{1:K}, b_{1:K}$ の同時分布はこうです。 $\displaystyle P(x^n, y^n, a_{1:K}, b_{1:K}) = \varphi(a_{1:K}, b_{1:K}) \prod_{i=1}^n p(x_i, y_i |a_{1:K}, b_{1:K})$ これをさらに観測データで条件付けた $P(y^n, a_{1:K}, b_{1:K} | x^n)$ を $q(y^n) r(a_{1:K}, b_{1:K})$ で近似したいのですが、 $P(y^n, a_{1:K}, b_{1:K} | x^n) \propto P(x^n, y^n, a_{1:K}, b_{1:K})$ なのでこのまま自己無矛盾条件に代入して差し支えありません。そうすると結局、 $q(y^n)$ は、各データの $K$ 面サイコロの目があるパターンのときだけ $1$ をとる確率分布になります（これは誤り）。 $r(a_{1:K}, b_{1:K})$ の方はハイパーパラメータを更新した形になります。それで「各データのサイコロの目のパターン」と「ハイパーパラメータ」を交互に更新していけばよかろうという算段のようなんです。

…この手続きで、「謎の確率変数 $Y$ 」と「交互の更新」がしれっと導入されましたが、そもそもこんなことをしてよいのでしょうか？「謎の確率変数 $Y$ 」は、いわば各データの実家が何番目の山かというようなものになっていますが、本当はどの山が実家とかないでしょう？そこには山の重なりのような分布があるだけで、各データがどの山からきたのかなんて人間の妄想に過ぎないはずなんです。それに、その「各データがどの山からきたのか」と「ハイパーパラメータ」を交互に更新するというのも、直感的には解に向かっていきそうですが、何の根拠があってそんなことをしているのか…。

2020-01-04 追記

「こんなことをしてもよいのか」どうかは、「何をするためにそうしたのか」によるだろう。今は何をしたいんだっけ。

何をしたいんだっけって、データから混合分布を推定したいんですよ。

推定するって何？

推定するって何って…この本では、4ページで定義した事後分布によって確率モデルを平均したものこそが推定分布でしたよね。だからこの定義に沿って事後分布を得ようとするのが筋のはずです。

$\displaystyle \begin{split} p(a_{1:K}, b_{1:K}|X^n) &= \frac{1}{Z_n(1)} \varphi(a_{1:K}, b_{1:K}) \prod_{i=1}^n p(X_i |a_{1:K}, b_{1:K} ) \\ &= \frac{1}{Z_n(1)} \varphi(a_{1:K}, b_{1:K}) \prod_{i=1}^n \sum_{k=1}^K a_k v(X_i) \exp \bigr(f (b_k) \cdot g(X_i) \bigl) \end{split}$ …山が1つのみのときは共役な事前分布がありましたが、山が複数あると山どうしのクロスタームが出てきてしまいますね…少なくとも山の個数だけ共役な事前分布を用意すればよいということにはなりません…。

まあ解析的に求まる例ではないんだよね。だからモンテカルロ法で解くとか、別の方法で解っぽい何かを探すことになる。5.2節の「平均場近似」はその「解っぽい何かを探す別の方法」だから、「4ページで定義した事後分布で確率モデルを得ること」はもはや放棄されている。「こんなことをしてもよいのか」というよりは「こうすることにした」んだね。…といっても、「だからどうしてそうすることにしたのか」の回答にはなっていないな。いま解析的には解けないけど「パラメータ $a_{1:K}, b_{1:K}$ の分布をよりよさそうなものに更新する」ことを達成したい。だから次善の策として平均場近似では別の分布でこう近似することにした。

本当の事後分布とのカルバック・ライブラー距離の極小点を目指すことにする。
極小点が満たすべき条件を更新式として利用できるようにするために、パラメータを2グループ（以上）に分けて、こちらのグループの分布を固定してあちらのグループの分布を更新する → 今度はさっき固定していたグループの分布を更新する → …という繰り返しができるようにする。この目的のため、パラメータの分布を独立な分布の積とする。

…実際に平均場近似が編み出された歴史的経緯はさておき、目的からブレークダウンするなら平均場近似というフレームワークは上のように解釈できるだろう。

なるほど、だから「謎の確率変数」が導入された…説明になっていないんですが。「交互の更新」の方にはそれで説明が与えられましたが。

うん、上のフレームワーク自体は「謎の確率変数の導入」を要請しないからね。ただ今回は、157ページの $\log P$ は「謎の確率変数」を導入しなくては「あるパラメータたちの分布を固定して別のパラメータたちの分布のハイパーパラメータを更新する」という形に持ち込めないはず。「謎の確率変数」を導入しないと結局山どうしのクロスタームが出ちゃうしね。

「謎の確率変数」の導入はそういう背景が…しかし、「じゃあこういう確率変数を導入しよう」というのも突飛な気がしますが…。

そんな突飛かな？「山どうしのクロスターム消えてほしい」→「各データの実家の山は一つまで」って思えば自然に導入できるんじゃない？

いうほど自然に導入できますかね？？それに、クロスタームを消してしまうということはそれだけ正しいやり方を損ねているでしょう？「各データの実家の山は一つまで」というのも結構大胆な仮定にみえるんですが、それが心配というか…。

もう少しつづく

2020-01-01

雑記：サンプリング定理の話（その2）

雑記

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。

参考文献

前回：雑記：サンプリング定理の話

前回は、関数 $x(t)$ が $f_s/2$ 以上の周波数成分を含まないならば、サンプリング周波数 $f_s$ でサンプリングした $x(n/f_s), \, n \in \mathbb{Z}$ のフーリエ変換の $[-f_s/2, f_s/2]$ の範囲を切り取って逆フーリエ変換することで元の $x(t)$ を復元できる、という話をしたね。じゃあぴったり $f_s/2$ の周波数成分についてはどうなるのかを調べると、余弦波は復元できて正弦波は復元できなかった。これは、正弦波はそもそもサンプリングする点にいつも「波」がないのだから捉えられないのも当然といえる。

それで前回の最後に副部長は「余弦波と正弦波をフーリエ変換してみてもわかる」といっていましたね。では余弦波からフーリエ変換してみますか。周波数 $f_0$ の余弦波は $x(t)=\cos(2 \pi f_0 t)$ ですから、フーリエ変換の定義にしたがって、

$\displaystyle \begin{split} X(f) &= \int_{-\infty}^\infty x(t) e^{-2 \pi i t f} dt = \int_{-\infty}^\infty \cos(2 \pi f_0 t) e^{-2 \pi i t f} dt = \int_{-\infty}^\infty \frac{e^{2 \pi i f_0 t} + e^{-2 \pi i f_0 t}}{2} e^{-2 \pi i t f} dt \\ &= \int_{-\infty}^\infty \frac{e^{-2 \pi i (f - f_0) t} + e^{-2 \pi i (f + f_0) t}}{2} dt = \frac{1}{2} \left[ \frac{e^{-2 \pi i (f - f_0) t}}{- 2 \pi i (f - f_0)}\right]_{t=-\infty}^{t=\infty} + \frac{1}{2} \left[ \frac{e^{-2 \pi i (f + f_0) t}}{- 2 \pi i (f + f_0)}\right]_{t=-\infty}^{t=\infty} \end{split}$ …負けました。副部長、余弦波にはフーリエ変換はありません。

投了しないで！余弦波にもフーリエ変換はあるよ！！

だって、上の積分は値をもちませんよ。フーリエ変換の定義は $X(f) = \int_{-\infty}^\infty x(t) e^{-2 \pi i t f} dt$ ですから、この積分ができない $x(t)=\cos(2 \pi f_0 t)$ はフーリエ変換をもたないと結論付けざるをえません。

そのフーリエ変換の定義は $x \in L^1(\mathbb{R}^n)$ である関数に対する定義だよ！ $L^1(\mathbb{R}^n)$ に属さない関数をフーリエ変換したいときは、定義を拡張するんだよ。つまり、その可積分関数に対するフーリエ変換の定義は、 $x(t), y(t)$ が2乗可積分でもあるならば以下のパーセバルの定理を満たす。「フーリエ変換は $L^2$ -内積を保つ」といっても同じだ。

$\displaystyle \int_{-\infty}^\infty x(t) \overline{y(t)} dt = \int_{-\infty}^\infty X(f) \overline{Y(f)} df$ これをフーリエ変換の新しい定義にすることができる。 $x(t)$ が可積分でなくても、可積分な $y(t)$ との内積が定義できる関数なら（というか関数でなくても $y(t)$ に作用してからそれを積分することができるようなものなら）、上式を満たす $X(f)$ を $x(t)$ のフーリエ変換とする（はず）。例えば $x(t) = e^{2 \pi i f_0 t}$ を左辺に代入してみよう。 $\displaystyle \begin{split} \int_{-\infty}^\infty x(t) \overline{y(t)} dt &= \int_{-\infty}^\infty \overline{ y(t) } e^{2 \pi i f_0 t} dt = \int_{-\infty}^\infty \overline{ y(t) e^{-2 \pi i f_0 t} } dt = \overline{ \int_{-\infty}^\infty y(t) e^{-2 \pi i f_0 t} dt } = \overline{Y(f_0)} \end{split}$ よって $\overline{Y(f_0)} = \int_{-\infty}^\infty X(f) \overline{Y(f)} df$ を満たす $X(f)$ が $x(t)$ のフーリエ変換だけど、この $X(f)$ は「積分することで隣の関数のある点での値を抜いてくるもの」、つまりデルタ関数だ。よって $X(f) = \delta(f - f_0)$ に他ならない（まあ以下の303番にもそうかいてあるし…）。

https://en.wikipedia.org/wiki/Fourier_transform#Distributions,_one-dimensional

$x(t) = e^{2 \pi i f_0 t}$ のフーリエ変換が $X(f) = \delta(f - f_0)$ …となると、 $x(t) = \cos(2 \pi f_0 t)$ のフーリエ変換は

$\displaystyle X_{\cos}(f) = \mathcal{F}\bigl[ \cos(2 \pi f_0 t) \bigr] = \frac{1}{2} \mathcal{F}\bigl[e^{2 \pi i f_0 t}\bigr]+ \frac{1}{2}\bigl[e^{-2 \pi i f_0 t}\bigr] = \frac{1}{2}\delta(f - f_0) + \frac{1}{2}\delta(f + f_0)$ ですね…。 $x(t) = \sin(2 \pi f_0 t)$ のフーリエ変換も同様です。 $\displaystyle X_{\sin}(f) = \mathcal{F}\bigl[ \sin(2 \pi f_0 t) \bigr] = \frac{1}{2i} \mathcal{F}\bigl[e^{2 \pi i f_0 t}\bigr]- \frac{1}{2i}\bigl[e^{-2 \pi i f_0 t}\bigr] = \frac{1}{2i}\delta(f - f_0) - \frac{1}{2i}\delta(f + f_0)$ まあ上のリンクの304番と305番にもこうかいてありますし…。しかし、余弦波のフーリエ変換も、正弦波のフーリエ変換も、2箇所に立つデルタ関数でできていて、どちらも似たような感じにみえますが。

余弦波と正弦波の周波数 $f_0$ をサンプリング周波数の半分 $f_s/2$ にして、 $X_{\cos}(f)$ と $X_{\sin}(f)$ を $f_s$ ずつずらして自身の「かげぶんしん」を重ねていくとどうなるかな。

ああ、いまは余弦波／正弦波からのサンプルを離散フーリエ変換することを考えているのでしたね。余弦波はこうです。

$\displaystyle \begin{split} \sum_{k = - \infty}^\infty X_{\cos}(f - k f_s) &= \frac{1}{2} \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s / 2) + \frac{1}{2} \sum_{k = - \infty}^\infty \delta(f - k f_s + f_s / 2) \\ &= \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s/ 2) \end{split}$ これ、2つの項がそれぞれ $1/2$ の高さ（？）のデルタ関数を櫛のように立たせていますが、立たせている場所が同じなので、結局、 $f=f_s/2$ を起点に $f_s$ おきに高さ $1$ のデルタ関数が立っていることになりますね。正弦波の方は、
$\displaystyle \begin{split} \sum_{k = - \infty}^\infty X_{\sin}(f - k f_s) &= \frac{1}{2i} \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s / 2) - \frac{1}{2i} \sum_{k = - \infty}^\infty \delta(f - k f_s + f_s / 2) \\ &= 0 \end{split}$ 2つの項が高さ $1/2i$ のデルタ関数の櫛と高さ $-1/2i$ のデルタ関数の櫛をちょうど同じ場所に建設しているので、相殺して消えますね…なんという負の奇跡…。

うん、だから余弦波の離散フーリエ変換は、区間 $[-f_s/2, f_s/2]$ を切り取って、区間の端は $1/2$ 倍すれば余弦波の離散ではないフーリエ変換の結果に一致する。正弦波の離散フーリエ変換は、区間 $[-f_s/2, f_s/2]$ を切り取る以前に $0$ だ。これじゃどうしようもない。

じゃあ $x(t) = \exp(\pi i f_s t)$ という関数だとどうなるかな？これは高くなったり低くなったりを繰り返す波ではなくて、時間 $2/f_s$ の間に複素数平面の単位円上を反時計回りに一周することを繰り返す、「回転する振り子」のようなものだね。ちなみにこれのフーリエ変換は $X(f) = \delta(f - f_s/2)$ だとさっきやった。

それなら離散フーリエ変換はこうですね。

$\displaystyle \sum_{k = - \infty}^\infty X_{\exp}(f - k f_s) = \sum_{k = - \infty}^\infty \delta(f - k f_s - f_s/ 2)$ ん？これ余弦波の離散フーリエ変換と一緒ですね…。複素数平面上を回転するものを時間 $1/f_s$ ごとに観測した結果と、実軸上を往復するものを時間 $1/f_s$ ごとに観測した結果とで、抽出した周波数成分が同じになってしまう…？いやでも、時刻 $n/f_s$ における座標がどちらも $\bigl( \cos(n \pi), \, 0 \bigr)$ で一致しているのでそれも道理な気はしますが。しかし、今回は正弦波のときのような相殺が起きたわけではありませんよね？？それなのに情報が欠落してしまうのでしょうか…？

サンプリング定理の手続きには、暗黙裡に「ちょうど $f_s/2$ の周波数成分については、『反時計回りに回転する成分』と『時計回りに回転する成分』が等しい（区間 $[-f_s/2, f_s/2]$ の端は $1/2$ 倍して切り取る）」という仮定が置かれている。だから、ちょうど $2/f_s$ の周期で複素数平面上を反時計回りに回転する振動子は、振幅の半分を「時計回りに回転する成分」にもぎとられてしまう。その結果、実軸上を往復する振動子に同一視されてしまう…って感じだと思うけどな。区間が重なる端っこの取り扱いはなんか決め打たないといけないし、対称だと考えるのが自然だろうしね。

結局、以下の関数をサンプリング周波数 $f_s$ でサンプリングして元の信号を復元するとこうなるね。正弦波は「消える」。

$x(t) = \cos(\pi f_s t)$ → $x(t) = \cos(\pi f_s t)$
$x(t) = \sin(\pi f_s t)$ → $x(t) = 0$
$x(t) = \exp(\pi i f_s t) = \cos(\pi f_s t) + i \sin(\pi f_s t)$ → $x(t) = \cos(\pi f_s t)$

サンプリング周波数のちょうど半分 $f_s/2$ の周波数成分については、「時刻 $0$ で複素数平面上のどこにいるか」の方向に射影した1次元の振動の成分しか生き残らない、とまとめることができるんじゃないかな。

しかし、サンプリング周波数のちょうど半分の周波数の正弦波が消えてしまうとはややこしいですね。いったいどうするのがいいのか…。

サンプリング周波数大きく取ればいいんじゃないかな。

はい。

おわり

2019-12-31

雑記：サンプリング定理の話

雑記

キャラクターは架空のものです。おかしい点がありましたらご指摘いただけますと幸いです。

参考文献

Nyquist–Shannon sampling theorem - Wikipedia

f:id:cookie-box:20190101155733p:plain:w60

Nyquist–Shannon sampling theorem - Wikipedia に以下のようにあります。

関数 $x(t)$ が $B$ ヘルツを超える周波数を含まないならば^※、 $x(t)$ は $1/(2B)$ ごとにサンプリングした値から完全に復元することができる。

※ ただし、ちょうど $B$ ヘルツの正弦波を含んではならない。「 $B$ ヘルツ以上の周波数を含まないならば」とされることも多い。

…サンプリングしたのに元の関数が「完全に復元できる」とは不思議ではありませんか？

あ、ちなみにこの定理、ウィキペディアの記事名は「ナイキスト‐シャノンのサンプリング定理」となっていますが、後にコテルニコフさんという方が先にこの定理を発表していたのが発見されたので、「ナイキスト‐シャノン‐コテルニコフのサンプリング定理」ともいうようですね。まるでカルーシュ‐クーン‐タッカー条件のようです。

f:id:cookie-box:20190101160814p:plain:w60

何にせよ $x(t)$ を完全に復元できるサンプリング周波数 $2B$ をナイキストレートというんだね。他方、実際のサンプリング周波数 $f_s$ の半分 $f_s/2$ をナイキスト周波数というのか…ややこしいな。サンプリング周波数 $f_s$ でサンプリングしたときは $f < f_s/2$ を満たす周波数成分は再構築されることが保証されると。だからなぜ保証されるかが知りたいんだけど…先にエイリアシングという節で次の式が出てくる。

$\displaystyle X_s(f) := \sum_{k=-\infty}^\infty X(f - k f_s) = \sum_{n=-\infty}^\infty \frac{1}{f_s} \cdot x \left( \frac{n}{f_s} \right) e^{- 2 \pi i n f / f_s}$ これは、「もし関数 $x(t)$ を $f_s$ ごとにサンプリングしたサンプル $x(n/f_s)$ が無限に手元にあったら、それに $e^{- 2 \pi i n f / f_s}$ をかけて足したものは、周波数の世界で間隔 $f_s$ ごとに元の関数のフーリエ変換 $X(f)$ を無限に重ねていったものと等しい」という式かな。これ自体はポアソンの和公式にサンプルたちを当てはめただけだ。中辺で表される周期関数から出発してフーリエ級数展開した係数を求めたらこうなったともいえる。

Poisson の和公式 - 記憶帳

えっと、だから何なんですか？私たちはいま、サンプル $x(n/f_s)$ たちから元の関数 $x(t)$ を復元できるのかに興味があると思うんですが。

だからだよ。最右辺をみたらこれは元の関数の「離散フーリエ変換」だ。上式の主張は、「ある関数を離散フーリエ変換すると、素直に周波数成分が出てきてくれるのではなく、周波数成分が無限に『かげぶんしん』したものが出てきてしまう」といっているに他ならない。ウィキペディアの下図のようにね。

AliasedSpectrum.png

私たちは「サンプルから元の関数が復元できるのか」に興味があるけど、これは現在の文脈でもっと正確にいうと「サンプルたちから元の関数の周波数成分の情報を得て、そこから元の関数が復元できるのか」だね。そして上図が示唆することは、「サンプリング周波数 $f_s$ の半分を超える周波数成分は復元できない」だ。なぜなら、周波数成分のグラフ（青色）の $f_s/2$ を超える部分は自身の「かげぶんしん」（黄緑色）と重なって一般的には識別できなくなってしまうからね。

ちなみに、画像や音声をデジタル化するときに高周波成分が低周波なノイズになってあらわれる現象を「エイリアシング」というみたいだね。これを防ぐために予め元の情報をローパスフィルタにかけるらしいよ。

折り返し雑音 - Wikipedia

「サンプリング周波数 $f_s$ の半分を超える周波数成分は復元できない」ですか…想像するに、1秒光って1秒消えるような周期2秒の点滅を繰り返すランプがあったとして、これを1秒おきに観測するなら「1秒おきに光っていそうだ」とわかりそうですが、2秒おきに観測するならもう点滅しているかもわからなさそうですね…。他方、対象の関数 $x(t)$ の周波数成分がサンプリング周波数 $f_s$ の半分を超えない範囲に収まっていれば、 $X(f)$ が識別できるのでそれをフーリエ逆変換すればよいということなんですね。

そうなんだけど、元の関数を直接サンプルたちで表す式も続きにあるね。上で定義した「かげぶんしん付き」フーリエ変換 $X_s(f)$ からかげぶんしんを排除するため、矩形関数で $f < |f_s/2|$ の部分だけ切り取る。

$\displaystyle \begin{split} X(f) &= {\rm rect}\left( \frac{f}{f_s} \right) X_s(f) \\ & = {\rm rect}\left( \frac{f}{f_s} \right) \sum_{n=-\infty}^\infty \frac{1}{f_s} \cdot x \left( \frac{n}{f_s} \right) e^{- 2 \pi i n f / f_s} \\ & = \sum_{n=-\infty}^\infty x \left( \frac{n}{f_s} \right) \cdot \frac{1}{f_s} {\rm rect}\left( \frac{f}{f_s} \right) e^{- 2 \pi i n f / f_s} \end{split}$ 矩形関数はsinc関数のフーリエ変換だから、両辺をフーリエ逆変換することで以下を得る。 $\displaystyle x(t) = \sum_{n=-\infty}^{\infty} x \left( \frac{n}{f_s} \right) \cdot {\rm sinc}\left( f_s t - n\right)$ sinc関数や矩形関数の定義は（そもそもフーリエ変換の定義も）、冒頭の記事からもリンクがあるけど、以下に準拠するよ。

https://en.wikipedia.org/wiki/Fourier_transform#Square-integrable_functions,_one-dimensional

「矩形関数はsinc関数のフーリエ変換だから」で片付けられても…本当にそのsinc関数をフーリエ変換すると矩形関数になりますか？

$\displaystyle \begin{split} \mathcal{F}[{\rm sinc}(f_s t - n)](f) &= \int_{-\infty}^\infty {\rm sinc}\left( f_s t - n\right) e^{-2 \pi i t f}dt \\ &= \int_{-\infty}^\infty {\rm sinc}\left( f_s t' \right) e^{-2 \pi i t' f} e^{-2 \pi i n f / f_s} dt' \qquad (t' = t - n/f_s) \\ &= e^{-2 \pi i n f / f_s} \int_{-\infty}^\infty {\rm sinc} \left( f_s t' \right) e^{-2 \pi i t' f} dt' \\ &= e^{-2 \pi i n f / f_s} \frac{1}{f_s} \int_{-\infty}^\infty {\rm sinc} \left( t'' \right) e^{-2 \pi i t'' f / f_s } dt'' \qquad (t'' = f_s t') \\ &= e^{-2 \pi i n f / f_s} \frac{1}{f_s} \int_{-\infty}^\infty \frac{\sin (\pi t'')}{\pi t''} e^{-2 \pi i t'' f / f_s} dt'' \end{split}$ この積分が ${\rm rect}(f/f_s)$ になるはずですが…難しそうですね。逆に矩形関数をフーリエ逆変換しますか。 $\displaystyle \begin{split} \mathcal{F}^{-1}[{\rm rect}(f)](t) &= \int_{-\infty}^\infty {\rm rect}(f) \cdot e^{2 \pi i t f} df \\&= \int_{-1/2}^{1/2} e^{2 \pi i t f} df \\&= \left[ \frac{e^{2 \pi i t f}}{2 \pi i t} \right]_{f=-1/2}^{f=1/2}\\&=\frac{e^{\pi i t } - e^{-\pi i t}}{2 \pi i t} \\ &=\frac{\sin (\pi t)}{\pi t} = {\rm sinc}(t) \end{split}$ …確かに矩形関数のフーリエ逆変換がsinc関数なので、sinc関数のフーリエ変換は矩形関数ですね。

ここまでの話を整理すると、

サンプリング周波数 $f_s$ で取ってきた関数の値のサンプルたちから、周波数成分の情報を得ることを経由して元の関数を復元したい。＝サンプルたちをフーリエ変換してから逆フーリエ変換することで元の関数を復元したい。
サンプルたちを素直に離散フーリエ変換すると $[-f_s/2, f_s/2]$ の範囲が延々と繰り返されたものになってしまうことがわかる。ので、元の関数の周波数成分がこの範囲に収まるようにサンプリング周波数 $f_s$ を大きくとり、かつ、サンプルたちの離散フーリエ変換に対して矩形関数で $[-f_s/2, f_s/2]$ の範囲を切り取ってから逆フーリエ変換する必要がある。
「サンプルたちの離散フーリエ変換に矩形関数をかけて逆フーリエ変換したもの」は、「サンプルたちにその位置でてっぺんになるsinc関数をかけて足したもの」に他ならない。

ということですか…ん、ちょうど $f_s/2$ の周波数成分はどうなるんです？ウィキペディアの Introduction のところに、「ちょうど $B$ ヘルツの正弦波を含んではならない」というようにありましたが、正弦波限定でだめなんですか？では余弦波は大丈夫だと？

…じゃあ部長に演習問題だ。以下の関数をサンプリング周波数 $f_s$ でサンプリングしてから復元してほしい。

$x(t) = \cos(\pi f_s t)$
$x(t) = \sin(\pi f_s t)$

ええ…1. は $x( n /f_s) = \cos(n \pi)$ なので、これは $n$ が偶数のとき $1$ 、 $n$ が奇数のとき $-1$ ですね。

$\displaystyle \begin{split} x(t) &=\sum_{n=-\infty}^\infty (-1)^n {\rm sinc}\left( f_s t - n\right)=\sum_{n=-\infty}^\infty (-1)^n \frac{\sin \bigl( \pi( f_s t - n) \bigr)}{\pi( f_s t - n)} \\&=\sum_{n=-\infty}^\infty (-1)^n \frac{\sin ( \pi f_s t) \cos(n \pi ) - \cos( \pi f_s t) \sin(n \pi ) }{\pi( f_s t - n)} = \sum_{n=-\infty}^\infty \frac{\sin ( \pi f_s t)}{\pi( f_s t - n)} \\&= \sin ( \pi f_s t) \sum_{n=-\infty}^\infty \frac{1}{\pi f_s t - n \pi} \\&= \sin ( \pi f_s t) \frac{\cos(\pi f_s t) }{\sin(\pi f_s t)} = \cos(\pi f_s t) \end{split}$ ここで最後の行への式変形は以下の記事の極展開によりました。

ミッタク＝レフラーの定理 - Wikipedia

$f_s t$ が整数に一致するような $t$ ではこの変形はできませんが、逆に $f_s t$ が整数に一致しているならば最初の行からただちに $x(t)=1$ （偶数に一致しているとき）または $x(t)=-1$ （奇数に一致しているとき）なので結局如何なる $t$ でも $x(t)=\cos(\pi f_s t)$ といえるのではないでしょうか。つまり、「サンプリング周波数の半分ぴったりの余弦波は完全に復元できる」といえそうです（自信はありませんが…）。

それでは 2. は、 $x( n /f_s) = \sin(n \pi)$ なので…え、これサンプルが全て $0$ ですね。だったら復元しても $x(t)=0$ にしかなりません…しかし、なぜ余弦波と正弦波でこのような違いが… 慢心、環境の違い…。

原点を起点にサンプリングしていくなら、正弦波は波のてっぺんの位置がちょうど悪くて（サンプリングするインターバルのちょうど中間にてっぺんがきてしまって、サンプリングする点ではいつも波が全くない）、余弦波は逆にちょうどいいってことだと思う。あと、余弦波と正弦波をフーリエ変換してみてもなんでこの違いが出るのかわかると思うよ。

つづくかも

f:id:cookie-box:20191231173321p:plain:w390

2019-12-22

雑記：レーティングのベイズ的解釈の話

雑記

キャラクターは架空のものです。お気付きの点がありましたらご指摘いただけますと幸いです。

一昨日、佐藤和俊七段が渡辺明三冠に勝利して、以下のサイトでの両者のレーティングが $14$ 変化しましたね。以下のサイトでのイロレーティングの定数値 $K$ は $16$ ですから、一回のレーティング変化の理論上限が $16$ であり、 $14$ というのは相当大きな変化ということになります。

将棋連盟棋士別成績一覧（レーティング）

今回の変化を表にまとめると以下です。

	12月20日対局前	12月20日対局後
渡辺明三冠レーティング $R_W$	$R_W^{(0)} = 2016$	$R_W^{(1)} = R_W^{(0)} - 16 W_{WS}^{(0)} = 2002$
佐藤和俊七段レーティング $R_S$	$R_S^{(0)} = 1696$	$R_S^{(1)} = R_S^{(0)} + 16 W_{WS}^{(0)} = 1710$
佐藤和俊七段が渡辺明三冠に勝利する確率 $W_{SW}$	$W_{SW}^{(0)} = 0.1368$	$W_{SW}^{(1)} = 0.1570$
渡辺明三冠が佐藤和俊七段に勝利する確率 $W_{WS}$	$W_{WS}^{(0)} = 0.8632$	$W_{WS}^{(1)} = 0.8430$

勝利確率はレーティングから算出される勝利確率 $\displaystyle W_{SW}^{(t)} = \frac{1}{10^{(R_W^{(t)} - R_s^{(t)})/400} + 1}$ ですね。 $400$ はただのスケーリング係数です。レーティングの更新規則も表の中に示しました。 $K=16$ は上記のサイトで採用されている定数値、いわば更新の度合いのハイパーパラメータですね。大きいほど直近の対局結果を反映するレーティングになります。つまりレーティングとは、負けた方から勝った方へ数値を移動させる、どれだけ数値を移動させるかは「勝った方が勝てなかったであろう確率」に比例させる、というパイの奪い合いです。そのため、パイを奪って引退する人が多いとレーティングはデフレし、奪われて引退する人が多いとレーティングはインフレします。将棋のプロ棋士の場合、引退間際には棋力が衰えてパイを奪われた状態で引退するケースが多いと思われるのでインフレ傾向にあると思われますが。無論、対局が強い対局者同士／弱い対局者同士に偏っていたり、棋士によって対局の頻度が違ったり、更新の度合いが状況の変化に追いつかないとレーティングは実態から乖離しますし、そもそも棋力は1次元に射影できるものでもありませんが…。

…さておき、上の表をみると、佐藤和俊七段が渡辺明三冠に勝利する確率（またはその逆）が、12月20日の対局の結果を受けて 13.68% から 15.70% に更新されていますが、まるでベイズ推測のようではないですか？であれば、更新の度合い $K$ とは、ベイズ逆温度 $\beta$ に相当するのではないでしょうか？？

…まぎれもなく12月20日の対局の結果を受けて更新されてはいるよね。ただそれがベイズ的かどうかは、ベイズ推定の枠組みにはまるかどうかという気がするけど…「佐藤和俊七段が渡辺明三冠に勝つか負けるか」の確率モデルはベルヌーイ分布だよね。なら共役事前分布はベータ分布 ${\rm Beta}(a, b)$ になる。以下の日記でやったコイン投げの例だね。

ベイズ統計の理論と方法：ノート2.5 - クッキーの日記

$\varphi(w) = {\rm Beta}(a, b)$ とすれば、まだデータを1つも観測していないときは $p^{\ast}( {\rm head} ) = a/ (a + b)$ 、 $p^{\ast}( {\rm tail} ) = b/ (a + b)$ だね。「表（勝ち）」のデータを1つ観測したときは以下のように更新される。

$Z_1(\beta) = \int_{0}^{1} w^{a - 1 + \beta} (1-w)^{b - 1} dw = \Gamma(a + \beta) \Gamma(b) / \Gamma(a + b + \beta)$
$p(w|X^1) = \bigl\{ \Gamma(a + b + \beta) / \Gamma(a + \beta) \Gamma(b) \bigr\} w^{a - 1 + \beta} (1-w)^{b - 1}$
$p^{\ast}( {\rm head} ) = \bigl\{ \Gamma(a + b + \beta) / \Gamma(a + \beta) \Gamma(b) \bigr\} \cdot \int_{0}^{1} w^{a - 1 + \beta + 1} (1-w)^{b - 1} dw = (a + \beta) / (a + b + \beta)$
$p^{\ast}( {\rm tail} ) = \bigl\{ \Gamma(a + b + \beta) / \Gamma(a + \beta) \Gamma(b) \bigr\} \cdot \int_{0}^{1} w^{a - 1 + \beta} (1-w)^{b - 1 + 1} dw = b/ (a + b + \beta)$

$\varphi(w) = {\rm Beta}(a, b)$ を事前分布にしてベイズ推測するのは、予め「表（勝ち）」を $a$ 回、「裏（負け）」を $b$ 回観測しておいたと考えて、観測したデータを $\beta$ 倍に「水増し」して最尤推定した結果と同じになるね。…それで、上の結果をイロレーティングと対比するとこうかな（ここで、本質的でない係数や対数の底を適当に変換したよ）。

	ベルヌーイ分布のベイズ推測	イロレーティング
対局前の勝利確率	$\displaystyle \frac{a}{a + b}$	$\displaystyle \frac{1}{\exp \bigl( R_W - R_S \bigr) + 1}$
対局後の勝利確率	$\displaystyle \frac{a + \beta}{a + b + \beta}$	$\displaystyle \frac{1}{\exp \left( R_W - R_S - \frac{2K}{\exp(R_S - R_W) + 1} \right) + 1}$

ちょっとみづらいから $r_{W} := \exp(R_W)$ 、 $r_{S} := \exp(R_S)$ とおくか。ウィキペディアのイロレーティングの記事をみるとわかるけど、レーティングのエクスポネンシャルは意味としては平均的なプレイヤーに対するオッズ（勝つ確率／負ける確率）だね。どちらかというと、オッズの対数がレーティングなんだけど。

	ベルヌーイ分布のベイズ推測	イロレーティング
対局前の勝利確率	$\displaystyle \frac{a}{a + b}$	$\displaystyle \frac{r_S}{ r_S + r_W}$
対局後の勝利確率	$\displaystyle \frac{a + \beta}{a + b + \beta}$	$\displaystyle \frac{r_S}{r_S + \exp \left( - \frac{2 K r_W}{r_S + r_W} \right) r_W}$ $\displaystyle = \frac{r_S \exp \left( \frac{2 K r_W}{r_S + r_W} \right) }{r_S \exp \left( \frac{2 K r_W}{r_S + r_W} \right) + r_W}$ $\displaystyle = \frac{r_S + r_S \left[ \exp \left( \frac{2 K r_W}{r_S + r_W} \right) -1 \right] }{r_S + r_W + r_S \left[ \exp \left( \frac{2 K r_W}{r_S + r_W} \right) -1 \right]}$

無理やり形を合わせにいってみた…変なことをしていない自信はない…。上の式変形を認めるなら、ベイズ推測とイロレーティングの対応は以下だな。

	ベルヌーイ分布のベイズ推測における各パラメータの役割	イロレーティングをベイズ更新と解釈すると対応するもの
事前分布 ${\rm Beta}(a, b)$ の $a$	予め相手に $a$ 回勝っていたと考える	自分の平均的プレイヤーへのオッズ $r_S$
事前分布 ${\rm Beta}(a, b)$ の $b$	予め相手に $b$ 回負けていたと考える	相手の平均的プレイヤーへのオッズ $r_W$
ベイズ逆温度 $\beta$	対局結果を $\beta$ 倍に水増しする	$\displaystyle r_S \left[ \exp \left( \frac{2 K r_W}{r_S + r_W} \right) -1 \right]$ $K$ は予め適当に決めた定数 $r_W/(r_S + r_W)$ は事前の負け確率

だから「更新の度合い $K$ はベイズ逆温度 $\beta$ に相当するのか」という問いへの答えは、ぴったり対応はしていない。でも、以下の極限での挙動は一致している。

$\beta \to 0$ も $K \to 0$ も、対局結果を全く取り入れないことに相当する。
$\beta \to +\infty$ も $K \to +\infty$ も、対局結果にしたがい勝った方の勝率を $1$ にすることに相当する。

イロレーティングにおいてベイズ逆温度 $\beta$ に相当するもの $\displaystyle r_S \left[ \exp \left( \frac{2 K r_W}{r_S + r_W} \right) -1 \right]$ は、更新の度合い $K$ の他に自分のオッズ $r_S$ と事前の負け確率 $r_W/(r_S + r_W)$ が入っているね。

自分のオッズ $r_S$ が大きいほどそれに比例して対局結果を水増しする。
更新の度合い $K$ が大きいほど指数関数的に対局結果を水増しする。
（勝ったとき）事前の負け確率 $r_W/(r_S + r_W)$ が大きいほど指数関数的に対局結果を水増しする。

ん？ $K$ や事前の負け確率が大きいほどレーティングが大きく変動するのはわかります。そのような定義式ですし。でも、自分のオッズ $r_S$ が大きいほどレーティングは大きく変化するということはありませんよ？レーティング変化は水準によらず、 $K$ に事前の負け確率をかけた数値です。

レーティングはね。でもレーティングと勝利確率は違うよ。例えば冒頭の対局結果は両対局者のレーティングが 400 ずつ低かったとしても対局前後の勝利確率は同じになるけど（勝利確率はレーティング差にしか依存しないからね）、レーティングが低い（オッズが小さい）のに同じだけの勝利確率の変化をさせたかったら、対局結果をオッズをかけて取り込まないといけないって感じかな。ただそもそも上の定義だと「予め相手に $a$ 回勝っていたと考える」の $a$ が自分のオッズだから自分のオッズが大きいほど強い事前信念をもっている感じになっちゃうんだよね…全体を自分のオッズで割ってからベイズ推測とイロレーティングを対応させた方がよかったかな…。

	ベルヌーイ分布のベイズ推測	イロレーティング
対局前の勝利確率	$\displaystyle \frac{a}{a + b}$	$\displaystyle \frac{1}{ 1 + \frac{r_W}{r_S}}$
対局後の勝利確率	$\displaystyle \frac{a + \beta}{a + b + \beta}$	$\displaystyle \frac{1}{1 + \exp \left( - \frac{2 K r_W}{r_S + r_W} \right) \frac{r_W}{r_S}}$ $\displaystyle = \frac{\exp \left( \frac{2 K r_W}{r_S + r_W} \right) }{\exp \left( \frac{2 K r_W}{r_S + r_W} \right) + \frac{r_W}{r_S}}$ $\displaystyle = \frac{1 + \left[ \exp \left( \frac{2 K r_W}{r_S + r_W} \right) -1 \right] }{1 + \frac{r_W}{r_S} + \left[ \exp \left( \frac{2 K r_W}{r_S + r_W} \right) -1 \right]}$

つづかない

2019-12-19

論文読みメモ： Latent Ordinary Differential Equations for Irregularly-Sampled Time Series（その1）

論文読み

以下の論文を読みます。

Yulia Rubanova, Tian Qi Chen, David K. Duvenaud. Latent Ordinary Differential Equations for Irregularly-Sampled Time Series. In Advances in Neural Information Processing Systems 32, 2019.
https://papers.nips.cc/paper/8773-latent-ordinary-differential-equations-for-irregularly-sampled-time-series

※ キャラクターは架空のものです。解釈の誤りは筆者に帰属します。おかしい点がありましたらご指摘ください。

前回：雑記： NeurIPS 2019 Proceedings の「時系列」を含むタイトル／次回：まだ

前回のあらすじです。

「不規則にサンプリングされた時系列のための隠れ常微分方程式」？不規則にサンプリングされた時系列は通常のRNNで扱うのは難しいと。それはそうですね。そこでRNNを、常微分方程式による連続的な隠れ状態をもつものに魔改造されたのですね。なるほど、だから「隠れ常微分方程式」ですか…。そもそもこれ以前に Latent ODE という手法が提案されているのですね。この ODE-RNNs も Latent ODE もポアソン過程から不規則にサンプリングされた系列を上手くモデリングしたようですね。この ODE を用いた手法は既存の RNN ベースの手法よりアウトパフォームしたということです。

ところでRNNには「等間隔でない」系列の特徴を学ぶのは難しいんですか？ RNNに入力するのって必ずしも時間方向に増えていく系列ではないと思うので等間隔というのがよくわからないんですが、例えば文章を形態素ごとに入力していくのとかは「等間隔」なんですか？

RNNは「最新の入力」と「直前までの特徴」から「最新の特徴」を出すモデルだよね。「直前までの特徴」を「最新の特徴」に更新できるだけの情報が各 $x_i$ に含まれていればいいと思うから、不等間隔というよりはそれで情報が欠落するのが問題なんじゃないかな。もちろん、情報の欠落がなかったとしても、あまりにばらばらな間隔の情報がどんどんモデルに投げ込まれてきたら、パターンが爆発して、モデル側も学習しづらいと思うけど。

うーん、それで、不等間隔なデータを取り扱う常套手段として、「等間隔になるようにプレ処理する」がありますが、これは情報を損ねると。特に情報が到着したタイミングなどを。それはそうですね。だから連続時間に対応できるモデルが望ましいといっています。それで、近年では観測時点の間で内部状態が指数関数的に減衰するような RNN（RNN-Decay）が発表されているそうです。指数関数的に減衰させるというのはそれはそれでいいんでしょうか？まあ、あまりに長い間観測値が入ってこなかったらRNNさんもいまどんな状況なのか全然わからなくなりそうなので、そういう意味ではもっている情報が減衰していくことになるんでしょうか…。また別のモデルとして、著者の方々のグループが過去に Neural ODEs というのを発表しているようですね。これは Figure 1 をみるに隠れ状態の各次元の時間変化を柔軟に記述しているようですが、どのようなものなんでしょうか。これは RNN ではなさそうですね。今回はその Neural ODEs で構築した、ODE による連続的な状態変化を RNN に組み込むといった感じなのですね。それが ODE-RNN だと…。それで2節に入ると、不等間隔なデータをRNNで取り扱う最もシンプルな方法は、前回の入力からの時間差をモデルに入れてしまうことだとありますね。確かに、これならRNNさんはいまの入力をどのように受け止めるべきか知ることができるかもしれません。

$h_i = {\rm RNNCell}(h_{i-1}, \Delta_t, x_i)\tag{1}$ しかしこれだと観測間の隠れ状態をどう定義するかに疑問が残る…？ RNNさんが出力するのは、「これまで受け取った系列の特徴はこうです」ですよね？観測間の隠れ状態を定義する必要があるんですか？？

だから、判定を常に出しっぱなしにしてほしいんだろうね。例えば、不定期的にジムにトレーニングに来るお客さんがいたとして、その人がもうジムに来なくなってしまう確率を毎日出しっぱなしにしたいとか？きっと日々変化していくよね。最後に来てから日が空くほど徐々に離脱確率は高まっていきそうだし、来てくれたらまた下がりそうだし。離脱確率が高い人が判定できれば確率が高い人たちだけに無料クーポンを配布するとかできそうだよね…まあそんな風に行動履歴をつかえるのかとか、クーポンが配布されなかったお客さんに不公平感か出ないかとかはありそうだけど。…これだと単純そうだから医療の例とかの方がいいのかな。色々な症状が不定期的に出て来院する患者さんが向こう1年に重大な病気を発症する可能性、とか？

急に中途半端に現実的になりましたね…ともかく、いまの問題設定は通常のRNNの想定と少しずれている気がします。通常のRNNってこうでしょう？

「最新の入力」と「直前までの特徴」から「最新の特徴」を出す。そういう仕組みで現在までに順番に入力されてきた入力たちの特徴を出す。

しかし、いまの文脈だとこうなのですね？

不定期的に入力される入力を受け取って常に特徴を出し続ける。

これだけ状況が変わると、まずリカレントにやるべきなんでしょうか…まあそれで、状態を連続的に出し続ける1つの方法は、ずっと同じ状態を出し続けることですね。しかしそれではよくないだろうということで次の RNN-Decay ですか。 $h_i = {\rm RNNCell} \bigl( h_{i-1} \cdot \exp(- \tau \Delta_t), x_i\bigr)\tag{2}$ 左辺が $h_i$ ですが連続時間の特徴を考えるなら添え字を $t$ にした方がいい気もしますが、とにかく $x_i$ が入力されたときに出力される状態（特徴）が $h_i$ ですね。しかし RNN-Decay では予測性能がよくならなかったと…何のタスクでなんでしょうか。そもそも RNN-Decay 自体何がやりたいのかよくわからないし…。

何がしたいのかも含めて読んでいけばわかるんじゃないかな。次に Neural ODEs の話が出てくるけど、隠れ状態 $h(t)$ を次の ODE の解とするんだね。

$\displaystyle \frac{dh(t)}{dt} = f_\theta \bigl( h(t), t \bigr), \; h(t_0) = h_0 \tag{3}$ この $f_\theta$ を学びさえすれば好きな時刻の $h(t)$ を手に入れられると。これはどのタイミングで学習してどのタイミングで利用するモデルなのかな…。まあそれで今回提案する ODE-RNN は3ページ目の上にある Algorithm 1 をみればいいのかな。つまり、各時刻 $t_i$ に、状態をまず時間発展させて、時間発展させた状態を通常のRNNのようにつかうんだね。

「最新の入力」と「直前の特徴をその時刻まで時間発展させたもの」から「最新の特徴」を出す。そういう仕組みで現在までに順番に入力されてきた入力たちの特徴を出す。

入力がない時刻についても時間発展させた状態が取り出せそうかな。RNN-Decay は減衰しかできなかったのを、より柔軟な時間発展が学習できるようにしたって感じだね。

（その2があれば）つづく