統計的学習の導入ぶぶん

\(y\) は決定的な値ではなく, 偶然によってふらつく値であると考える.
全く同じ条件で観測しても, あるときは 5 と観測され, あるときは 7 と観測されたりする.
一つひとつの \(y\) の実現値はそんな値.
すなわち, (\(x\) を固定して)全く同じ条件で何回も観測すれば, \(y\) のばらつき具合(分布)を見ることが出来る筈である.
実際は, 全く同じ条件で何度も観測できてる状態はまれなので, \(x\) を固定しない \(y\) のデータ全体から, \(x\) を固定した状態での \(y\) の分布を想定することになる(?).

例えば, 大体いつも出てくる典型的な想定は \(y\) は正規分布をする, というもの.
観測値 \(y\) のばらつきがいろんな細かい条件の積み重ねに由来するものであれば, \(y\) は正規分布すると考えて問題無いことが多かったりする(中心極限定理. だったっけ？).

例: \(y\) は正規分布してる, と考えた場合, \(y\) の確率密度 \(p(y)\) は以下のようにかける.
$$ p(y) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp(-\frac{(y-m)^2}{2\sigma^2}) $$

正規分布の場合, 分布の具体的な形を決めるためのパラメータは 2 つあって, 一つは期待値 \(m\), もう一つは標準偏差 \(\sigma\) である.

想定した分布のパラメータ(と \(y\) の代表値との関係)を把握したあと次に進む.

想定される分布は, \(y\) がどのような値をとりえるかによってある程度制限される.
\(y\) の真の分布はわからないので(うまくデータと合うようにできれば実際はなんでもよい？(注)),
\(y\) の値の特性によって用いる分布に必要条件的に制限がかけられている, というイメージでかんがえるとよい.

(注)
Stan と R でベイズ統計モデリング p80 より

    ほとんどの場合, 現実のデータがある確率分布に従っていると考えるというのは単なる仮定にすぎない.

以下によく見る分布をしめす. 各々の分布についての詳細は省略する.

\(y\) が連続値をとる場合:

\(y\) の範囲	分布	備考
\((-\infty,\infty)\)	正規分布
\((0,\infty)\)	ガンマ分布	待ち時間とか
\((0,\infty)\)	対数正規分布	体重とか
\((0,1)\)	ベータ分布	比率
\([0,\infty)\)	指数分布	イベントが起こるまでの時間

\(y\) が \(K\) 次元のベクトルの場合:

\(y_{k}\) の範囲	\(\vec{y}\) の条件	分布	備考
\((-\infty,\infty)\)		多変量正規分布
\((0,1)\)	\(\sum y_k = 1\)	ディリクレ分布	ベータ分布の多次元版

\(y\) が離散値をとる場合:

y の範囲	分布	備考
\([0,1]\)	ベルヌーイ分布	コイン投げ
\([0,N]\)	二項分布	複数のコイン投げで表の出た数
\([0,\infty)\)	ポアソン分布	イベント回数/単位時間
\([1,K]\)	カテゴリカル分布	サイコロ投げで出た目

\(y\) が \(K\) 次元のベクトルの場合:

\(y_{k}\) の範囲	\(\vec{y}\) の条件	分布	備考
\([0,N]\)	\(\sum y_{k}=N\)	多項分布	複数のサイコロ投げで \(k(k=1,..K)\) が出た数の分布

(注):
ベルヌーイ分布(0,1)の多次元版(1..K)がカテゴリカル分布
二項分布(0,1 をN 回)の多次元版(1..K をN 回)が多項分布
ベルヌーイ分布(0,1)を多数回やったのが二項分布(0,1 を N 回)
カテゴリカル分布(1..K)を多数回やったのが多項分布(1..K をN回)
参照:
ややこしい離散分布に関するまとめ - 作って遊ぶ機械学習。

\(y\) の確率分布 \(p(y)\) は, パラメータ \(m\), \(\sigma\) (正規分布の場合) の具体的な値が決まれば決るので, この具体的な値を求めることが目標になる.
\(y\) の分布を正規分布と想定することで, 求めるべきものは \(y\) そのものではなくパラメータ \(m\), \(\sigma\) になった.

ところで, 変数 \(y\) は 変数 \(x\) と関係している筈だからそれらの関係性を知りたい, というのが元の問いである.
上で, 考えるのは \(y\) そのものではなく分布のパラメータ \(m\), \(\sigma\) に移行しているので, \(m\), \(\sigma\) と \(x\) との関係性について考えることになる.

\(m\) は \(y\) の期待値という代表値の一つを表しており, これが \(x\) の値と関係するのでは, と考えるのはまぁ自然なことかなーと思うけどどうだろうか？
例えばこんな感じ.

例: \(y\) の代表値(期待値) \(m\) が \(x\) が大きくなるにつれて大きくなる(あるいは小さくなる)と考えた場合:
$$ m = \alpha x + \beta $$

この場合, パラメータ \(m\) は データ(\(x\)) 毎に異なった値と仮定されて, そのパラメータは \(\alpha\), \(\beta\) となっている.

この段階で, 求めるべき変数は \(m\), \(\sigma\) から \(\alpha\), \(\beta\), \(\sigma\) の 3 つになった(増えちゃった..).
上の例のように, 普通は, \(y\) の代表値を示すようなパラメータと \(x\) との関係を考えることが多いようである.

上の例では正規分布の平均値を表すパラメータ \(m\) と \(x\) が線形の関係であると仮定して,
$$ m = \alpha x + \beta $$
とした.
\(m\) の値の範囲は \((-\infty,\infty)\) なので, 右辺をどんなふうに仮定しもとりあえず矛盾することはないかと思われる.

しかし, 正規分布以外の場合には, 推定するパラメータの値の制限によって, \(x\) との関係性(右辺)はある程度制限されることがある.

以下によく見る分布とパラメータの制限, 代表的な \(x\) との関係を表す関数についてまとめておく.

分布	パラメータの制限	\(x\) との関係を表す関数	備考(glm)
正規分布	\(m=(-\infty,\infty)\), \(\sigma \geq 0\)	平均値 \(m=ax+b\)	identity
ガンマ分布	\(\alpha=(0,\infty)\), \(\beta =(0,\infty)\)	平均値 \(\frac{\alpha}{\beta} =\exp(ax+b)\)	log (注1)
対数正規分布	\(m=(-\infty,\infty)\), \(\sigma=(0,\infty)\)	平均値 \(m=ax+b\) ?	(注2)
ベータ分布	\(\alpha=(0,\infty)\), \(\beta =(0,\infty)\)	平均値 \(\frac{\alpha}{\alpha+\beta}=\frac{1}{1+\exp(-(ax+b))}\) ?	(注3)
指数分布	\(\beta=(0,\infty)\)	平均値 \(\frac{1}{\beta}=\frac{1}{exp(ax+b)}\) ?	(注3)
多変量正規分布	\(m_k=(-\infty,\infty)\)	\(m_k = a_k x + b_k\) (平均値 \(\vec{m}\))
ディリクレ分布	\(\alpha_k=(0,\infty)\)	\(\alpha_k=\exp(a_k x + b_k)\) ? (平均値 \(\frac{\alpha_k}{\sum \alpha_k}\))	(注3)
ベルヌーイ分布	\(\theta=[0,1]\)	確率 \(\theta=\frac{1}{1+\exp(-(ax+b))}\)	logit
二項分布	\(\theta=[0,1]\)	確率 \(\theta=\frac{1}{1+\exp(-(ax+b))}\)	logit
ポアソン分布	\(\lambda=(0,\infty)\)	平均値(分散) \(\lambda=\exp(ax+b)\)	log (注4)
カテゴリカル分布	\(\theta_k=[0,1]\), \(\sum \theta_k =1\)	確率 \(\theta_k = \frac{\exp(a_k x+b_k)}{\sum \exp(a_k x+b_k)}\)	softmax
多項分布	\(\theta_k=[0,1]\), \(\sum \theta_k =1\)	確率 \(\theta_k = \frac{\exp(a_k x+b_k)}{\sum \exp(a_k x+b_k)}\)	softmax

(注1):
ガンマ分布の場合, \(\alpha(shape)\) は一定値で, \(\beta(rate)\) が \(x\) と関係する, と考えることが多い?
glm の default の link 関数は inverse (http://www.statmethods.net/advstats/glm.html)
(注2):
対数正規分布の場合, glm では, family=gaussian で link 関数を log にする(でいいのかな?).
http://qiita.com/gigamori/items/acf6b45b4a5cc2f604ce
(注3): よくわかんないので適当.
以下の URL も参照
重回帰分析＆一般化線形モデル選択まわりの再まとめ 六本木で働くデータサイエンティストのブログ
Stanコードの書き方　中級編:
(注4): \(\log(\lambda)=ax+b\)

パラメータ値の制限と使われる典型的な関数の関係

パタメータ制限	関数
\((-\infty,\infty)\)	\(ax+b\)
\((0,\infty)\)	\(\exp(ax+b)\)
\([0,1]\)	\(\frac{1}{1+\exp(-(ax+b))}\)
\(\sum = 1\)	\(\frac{\exp(a_k x+b_k)}{\sum \exp(a_k x+b_k)}\)

まずは 1 つ目のデータ \((X^{(1)}, Y^{(1)})\) が得られる確率(密度)は, \(p(y)\) が(今回の例では)以下のように書けると仮定してるので,
$$ p(y) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp(-\frac{(y-m)^2}{2\sigma^2}) $$

\(y\) に具体的な値 \(Y^{(1)}\) を代入して
$$ p(y=Y^{(1)}) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp(-\frac{(Y^{(1)}-m)^2}{2\sigma^2}) $$
一方 \(m\) については
$$ m = \alpha x + \beta $$
と考えているので,
\(x\) に具体的な値 \(X^{(1)}\) を代入して,
$$ m = \alpha X^{(1)} + \beta $$
となって, 結局
$$ p(y=Y^{(1)}) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp(-\frac{(Y^{(1)}-\alpha X^{(1)} -\beta)^2}{2\sigma^2}) $$
となった.

データ \(2, 3, \ldots N\) についても同様に計算できる.

\(p(y=Y^{(1)})\) の式には, 未知な変数(パラメータ)として \(\alpha\), \(\beta\), \(\sigma\) がある.
\(Y^{(1)}\), \(X^{(1)}\) は具体的な値で, 例えば 3 とか 5 とかが入っている.
この式は, \(\alpha\), \(\beta\), \(\sigma\) の関数として見ることができる.
つかわれてる数(変数, データ)を明示的にするため以下のように書いておく.

$$ l(\alpha,\beta,\sigma, y=Y^{(1)}, x=X^{(1)}) = p(y=Y^{(1)}) $$
\(y\) についてみるとこの値は確率(密度)である(\(y\) について積分すると 1 になる).
しかし \(\alpha\), \(\beta\), \(\sigma\) から見ると確率ではなく普通の関数である(規格化されてない)とみなせる.
\(\alpha\), \(\beta\), \(\sigma\) の関数としてみたとき, この値は「尤度(ゆうど)」と呼ばれている.

さて, \(p(y=Y^{(1)})\) は データ \(Y^{(1)}\) が得られる確率(密度)だが, データはさらに \(Y^{(2)}\), \ldots \(Y^{(N)}\) まで, 全部で \(N\) 個ある.
ここで, \(N\) 個のデータが得られる同時確率(密度)について考えてみる.
各々のデータ \(Y^{(i)}\) はそれぞれ独立(\(Y^{(1)}\) の値が \(Y^{(2)}\) の値の出方(確率)に影響しない) であると普通は考えられるので, 同時確率(密度) \(q\) は
$$ q = \Pi_{i=1}^{N} p(y=Y^{(i)}) $$
となる.
$$ p(y=Y^{(i)}) = l(\alpha,\beta,\sigma, y=Y^{(i)}, x=X^{(i)}) $$
なので \(q\) は以下のように書ける.
$$ q = \Pi_{i=1}^{N} l(\alpha,\beta,\sigma, y=Y^{(i)}, x=X^{(i)}) $$
使われてる数(パラメータ, データ)を明示して書くと以下のようになる.
$$ q(\alpha,\beta,\sigma, Y^{(1)},X^{(1)},...,Y^{(N)},X^{(N)}) = \Pi_{i=1}^{N} l(\alpha,\beta,\sigma, Y^{(i)}, X^{(i)}) $$
\(Y^{(1)}\), \(X^{(1)}\) … には具体的な数値がはいっている.

\(q\) は掛け算の繰り返し(\(\Pi\))があって取扱いにくいので, 計算の際には普通は上式の対数をとったものが使われる.
対数をとると掛け算(\(\Pi\)) は足し算(\(\sum\)) になる.
\(L = \log q\) とかいて,

$$ L(\alpha,\beta,\sigma, Y^{(1)},X^{(1)},...,Y^{(N)},X^{(N)}) = \sum_{i=1}^{N} \log l(\alpha,\beta,\sigma, Y^{(i)},X^{(i)}) $$
これを「対数尤度」をいったりするらしい.
基本これを元に色々な計算を行っていく.

ベイズ推定では, パラメータの「真の値」を推定するのではなくパラメータの「分布」(事後分布)を求める.

ベイズの公式:
$$ P(A,B) = P(A|B)P(B) = P(B|A)P(A) $$

$$ P(A|B) = \frac{P(B|A)P(A)}{P(B)} $$

\(A\) を パラメータ(\(\alpha\),\(\beta\),\(\sigma\)), \(B\) をデータ(\(X^{(i)}\),\(Y^{(i)}\) \(i=1..N\)) とすると,
\(P(B|A)\) (パラメータ \(A\) が与えられたもとでの データ \(B\) の確率)は以下の \(q\) に相当する.
$$ q(\alpha,\beta,\sigma, Y^{(1)},X^{(1)},...,Y^{(N)},X^{(N)}) = \Pi_{i=1}^{N} l(\alpha,\beta,\sigma, Y^{(i)}, X^{(i)}) $$
これは \(y\) にとっては確率で, パラメータ(\(\alpha\),\(\beta\),\(\sigma\))にとっては尤度と呼ばれる値でもある.

\(P(A)\) は事前確率と呼ばれていて, パラメータに関するの事前知識に基づいた偏りを示す.
何も知識がなければ定数にしておく(ほんとは確率(密度)だから足して(積分して) 1 にならないといけないんだけど…).

\(P(B)\) はデータが得られる確率で, \(A\) とは関係ない値.
今の状況ではデータ \(B\) は具体的な値なので, \(P(B)\) も何か具体的な値となる.
つまりこいつは 5 とか 8 とか 1000 とかの定数値となるはず.

まとめると,
$$ (事後分布) = (尤度) * \frac{(事前分布)}{(定数)} $$

事前分布が定数であれば事後分布は尤度と同じ(定数分だけ異なる)で, 目的関数は最尤推定とほぼ同様になるが, 最尤推定は最大の尤度をもつパラメータ 1 つをゲットするのに対して,
ベイズ推定ではパラメータの「分布」を考えるところが違うのであった.

事後分布の推定は, この分布に従う多数の粒子を発生させる方法(MCMC)を用いて数値的に行う.
コンピュータがいい感じでやってくれると思う.

ここからは, 求めるパラメータ(上の例では \(\alpha\), \(\beta\), \(\sigma\)) を \(\theta\) という一つの文字で代表させる.
$$ p(\theta|y) = \frac{p(y|\theta) p(\theta)}{p(y)} $$

\(p(\theta)\) を定数とすれば(= \(\theta\) に関して何の仮定もおかない),
$$ \log p(\theta|y) = \log p(y|\theta) + const(\theta) - const(y) $$
なので,

\begin{align*} \log p(\theta|y) &\propto& \log p(y|\theta) \end{align*}

となる.

一方, \(p(\theta) \sim normal(0,\tau)\) (正規分布) とすれば(= \(\theta\) の平均値は 0 っぽい, という仮定をおく. \(\tau\) が小さいほど \(\theta=0\) を確信してる, というイメージ),

\begin{align*} p(\theta) &=& \frac{1}{\sqrt{2\pi\tau^2}} \exp(-\frac{(\theta-0)^2}{2\tau^2}) \\ \log p(\theta) &\sim& - \frac{\theta^2}{2 \tau^2} \equiv - \gamma \theta^2 \end{align*}

よって,

\begin{align*} \log p(\theta|y) &\propto& \log p(y|\theta) + \log p(\theta) \\ &\propto& \log p(y|\theta) - \gamma \theta^2 \end{align*}

となる.
ここで, \(\gamma\) は分散に関係するパラメータ(計算前に自分で適当に決める値. \(\gamma\sim\frac{1}{\tau^2})\) である.

右辺第二項が追加された状態での回帰はリッジ(Ridge)回帰あるいはチノコフの正則化と呼ばれている.
右辺第二項にパラメータの制限を加えることを正則化というらしいのだが, これはパラメータの事前分布を上のように設定してベイズ的に推定してるのと同じ.

事前分布にラプラス分布
$$ p(y) = \frac{1}{2\sigma} \exp(-\frac{|y-m|}{\sigma}) $$
を使ったものがラッソ(Lasso)回帰と言われるものである(見た感じラプラス分布だと思うんだけど要確認です).

このあたりの超わかりやすい説明は, プログラミングのための確率統計 8.1.回帰分析と多変量解析から p277 に書いてあった.
正則化の説明をどっかでみたときに, 後ろに勝手に適当な関数を追加したりして, そんなことやっていいの? とか思ったのだが, ここを読んでやっと意味がわかったような気がした.

パラメータ \(\alpha\) は, 最尤推定ではすべてのデータに共通した一つの値として求まる..
ベイズ推定ではすべてのデータに共通した一つの分布として求まる.
階層ベイズでは, パラメータ \(\alpha\) はデータあるいはグループ毎に異なる複数の分布からなると考える.
ここでは, \(\alpha\) はグループに関するデータ \(z\) に依存すると考え, \(y\) について考えたのと同様に以下の 3 点を考えてみる.

\(y\) は正規分布と仮定.
$$ p(y) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp(-\frac{(y-m)^2}{2\sigma^2}) $$

\(y\) のパラメータ \(m\) (期待値)が, データ \(x\) と線形関係にあると仮定.
$$ m = \alpha x + \beta $$

データ \((X^{(i)},Y^{(i)})\) の確率は,
$$ l(\alpha,\beta,\sigma, Y^{(i)}, X^{(i)}) = p(Y^{(i)}) $$
同時確率(尤度) \(q\) は
$$ q(\alpha,\beta,\sigma) = \Pi_{i=1}^{N} l(\alpha,\beta,\sigma, Y^{(i)}, X^{(i)}) $$
対数尤度 \(L\) は
$$ L(\alpha,\beta,\sigma) = \sum_{i=1}^{N} \log l(\alpha,\beta,\sigma, Y^{(i)},X^{(i)}) $$
この部分の計算を行うと,
$$ -\sum_{i=1}^{N} (Y^{(i)}- (\alpha X^{(i)} + \beta))^2 $$
という形になる.
これは最小二乗法として知られている関数の形です.
つまり, 最小二乗法を使った直線の当てはめ = \(y\) に正規分布を仮定 してるのと同じ.

\(y\) は 0 or 1 であり, ベルヌーイ分布を想定するのが自然だろう.
$$ p(y) = t^y (1-t)^{(1-y)} $$

\(y\) のパラメータは \(t\) で, これは \(y=1\) である確率という意味を持つ.
確率なので, \(0 \leq t \leq 1\) です.
\(x\) は一般に \(-\infty < x < \infty\) なので, ロジスティック関数
$$ f(x) = \frac{1}{1 + \exp(-x)} $$
を使って,
$$ t = f(\alpha x + \beta) $$
とし, \(t\) を 0 から 1 の間に押し込めることにする.

データ \((X^{(i)},Y^{(i)})\) の確率は,
$$ l(\alpha,\beta, Y^{(i)}, X^{(i)}) = p(Y^{(i)}) $$
同時確率(尤度) \(q\) は
$$ q(\alpha,\beta) = \Pi_{i=1}^{N} l(\alpha,\beta, Y^{(i)}, X^{(i)}) $$
対数尤度 \(L\) は
$$ L(\alpha,\beta) = \sum_{i=1}^{N} \log l(\alpha,\beta, Y^{(i)},X^{(i)}) $$
この部分の計算はどんな感じかというと,

\begin{align*} L &=& \sum_{i=1}^{N} [ t \log(Y^{(i)}) + (1-t) \log(1-Y^{(i)}) ]\\ &=& \sum_{i=1}^{N} [ f(\alpha X^{(i)} + \beta) \log(Y^{(i)}) + (1-f(\alpha X^{(i)} + \beta)) \log(1-Y^{(i)}) ] \end{align*}

\(- A \log B\) の和 = 交差エントロピー(cross entropy) と言う.
L の最大化 = 交差エントロピーの最小化.

1-of-K 表記法を利用する.
この記法では,
$$ y = (y_1, y_2, .. y_k, .. y_K) $$
\(y_k\) は, クラス \(k\) であれば 1, それ以外は 0 である.
(\(y_k\) のうち一つだけが 1, それ以外は 0).
\(y\) について以下の式が成り立つ.
$$ \sum_{i=1}^{K} y_k = 1 $$

\(p(y)\) は以下のカテゴリカル分布と想定される.
$$ p(y) = \Pi_{k=1}^{K} t_{k}^{y_k} $$
\(t_k\) はクラス \(k\) である確率を示す.

\(t\) について以下の式が成り立つ.
$$ \sum_{i=1}^{K} t_{k} = 1 $$

\(t_k\) は確率なので \(0 \leq t_k \leq 1\) を満たし,
なおかつ
$$ \sum_{i=1}^{K} t_{k} = 1 $$
も満たさないといけない.

これを満たすものとして,
$$ t_k = \frac{\exp(\alpha_k x + \beta_k)}{\sum_{k=1}^{K} \exp(\alpha_k x + \beta_k)} $$
例えばこんなのを考える事ができる.
ソフトマックス関数と呼ばれてるらしい.

ロジスティック回帰とほぼおなじ.
対数尤度に 1 から \(K\) の和が加わる.
$$ L(\alpha,\beta) = \sum_{k=1}^{K} \sum_{i=1}^{N} \log l(\alpha_k,\beta_k, Y^{(i)},X^{(i)}) $$

重み \(w = \vec{w}\) とバイアス \(b\) をパラメータとして \(m = \vec{w} \vec{x} + b\)
とかる. 上の例と全く一緒である.

\begin{align*} \mu &=& \vec{w} \vec{x} + b \\ t &=& \frac{1}{1 + \exp(-\mu)} \end{align*}

$$ \vec{\mu} = \vec{w} \vec{x} + \vec{b} $$
\(\vec{\mu}\) は K 次元ベクトル.
\(\vec{w}\) は K 行 N 列の行列.
\(\vec{x}\) は N 次元ベクトル.
\(\vec{b}\) は K 次元ベクトル.

$$ t_k = \frac{\exp(\mu_k)}{\sum_{k=1}^{K} \exp(\mu_k)} $$