2021-06-01

一様分布の上限の推定量の比較

さて、今回は一様分布の上限の推定方法の比較を行います。

初めに簡単に問題設定から。
$x_1,x_2,...,x_n$ が独立に、範囲 $0$ ≦ $x_i$ ≦ $\theta$ の一様分布に従うような状況を考えます。
この時、最尤推定法を用いて、 $\theta$ の推定量 $\hat{\theta}$ を考えると、 $\hat{\theta}=max\{x_1,...x_n\}=x_{max}$ になるのでした。

ここまでの内容に関しては、既に以下の記事で触れていますので、初耳！という方は以下をご覧ください。
www.bananarian.net

今回は、 $\theta$ の推定量は本当に $x_{max}$ で良いのか。もっとよりよい推定量は存在しないのかについて検討をします。

最大値の分布の導出
最大値の期待値と分散の導出
一様分布の上限値の不偏推定量
一様分布の上限値の不偏推定量 part2
一様分布の上限値の不偏推定量まとめ

最大値の分布の導出

まず、 $x_{max}$ は一体どんな分布をしていて、期待値や分散はいくらなのでしょうか。まずはそれらを計算します。

$M=max(X_1,...X_n)$ とします。この時 $M$ の分布関数は定義から、以下のように書けますね。

$P(M≦m)=P(max(X_1,...X_n)≦m)$

ここで、 $max(X_1,...X_n)≦m$ とは、全ての $X_i$ が $m$ 以下であることと同値であり、更に $X_i$ は互いに独立であるため、

$P(M≦m) = P(X_1≦m)×P(X_2≦m)×…×P(X_n≦m)$
$= (\int^{m}_{0} \frac{1}{\theta} dx)^n$

が成り立ちますね。 $\frac{1}{\theta}$ は一様分布の密度関数です。

これを計算すると、 $P(M≦m)=(\frac{m}{\theta})^n$ となります。

よって、 $m$ について微分を施すと、 $M=x_{max}$ の密度関数 $f(m)$ が得られ、

$f(m) = \frac{n}{\theta} (\frac{m}{\theta})^{n-1}$ とわかります。

最大値の期待値と分散の導出

密度関数がわかると、勿論期待値や分散を計算することが出来ます。

$E[x_{max}] = \int^{\theta}_0 x \frac{n}{\theta} (\frac{x}{\theta})^{n-1} dx = \frac{n}{n+1} \theta$
$E[x_{max}^2] = \int^{\theta}_0 x^2 \frac{n}{\theta} (\frac{x}{\theta})^{n-1} dx = \frac{n}{n+2} \theta^2$

よって、分散は以下のようになります。
$V[x_{max}] = E[x_{max}^2] - E[x_{max}]^2 = \frac{n}{n+2} \theta^2 - \frac{n^2}{(n+1)^2} \theta^2 = \frac{n}{(n+2)(n+1)^2} \theta^2$

一様分布の上限値の不偏推定量

期待値を導出してみると、どうやら $x_{max}$ の期待値は $\theta$ ではないということに気付くかと思います。
最尤推定量は必ずしも不偏性を持つわけではないので、このような状況は多々あります。

そこで不偏性を持たせるために以下のような修正を施してみます。
$\hat{\theta}_2 = \frac{n+1}{n} x_{max}$

これは、先ほどの計算から明らかですが、以下のようになるため不偏推定量です。
$E[\hat{\theta}_2] = E[\frac{n+1}{n} x_{max}] = \frac{n+1}{n} E[x_{max}] = \theta$

分散に関しても先ほどの計算から簡単に計算することが出来ます。
$V[\hat{\theta}_2]= V[\frac{n+1}{n} x_{max}] = \frac{(n+1)^2}{n^2} V[x_{max}] = \frac{1}{n(n+2)} \theta^2$

一様分布の上限値の不偏推定量 part2

そもそも $x_{max}$ に縛られずとも、一様分布の期待値は $E[x_i] = \frac{\theta}{2}$ ですから、 $\theta$ の推定量として次のようなものを考えることも出来ますね。

$\hat{\theta}_3 = \frac{2}{n} \sum_{i=1}^n x_i$

この推定量の期待値と分散を計算すると以下のようになります。

$E[\hat{\theta}_3]=\theta$
$V[\hat{\theta}_3] = \frac{1}{12n} \theta^2$

一様分布の上限値の不偏推定量まとめ

以上3種類の推定方法を紹介しました。まとめると以下のようになります。

$\hat{\theta}=max\{x_1,...x_n\}=x_{max}$
$\hat{\theta}_2 = \frac{n+1}{n} x_{max}$
$\hat{\theta}_3 = \frac{2}{n} \sum_{i=1}^n x_i$

さて、色々と出てきましたが、じゃあどれを使うべきかについて考えてみます。

n→∞のケース

nを無限に飛ばした場合、まず最尤推定量である $\hat{\theta}$ は一致性を持つので、真のパラメータ $\theta$ へ収束します。
更に、 $\hat{\theta}_2$ も、nを無限に飛ばせば $\hat{\theta}$ と等しいため、真のパラメータへ収束します。
$\hat{\theta}_3$ も、大数法則から真のパラメータへ収束します。

では、分散はどうでしょう。
$V[\hat{\theta}] = \frac{n}{(n+2)(n+1)^2} \theta^2$
$V[\hat{\theta}_2]= \frac{1}{n(n+2)} \theta^2$
$V[\hat{\theta}_3] = \frac{1}{12n} \theta^2$

nを無限に飛ばした場合、 $\hat{\theta}$ と $\hat{\theta}_2$ は分母に $n^2$ があるので、分母が $n$ の $\hat{\theta}_3$ よりもはやい速度で値が収束することが直感的に分かるかと思います。

よって、nがある程度大きい場合は $\hat{\theta}$ と $\hat{\theta}_2$ を使うのが良さそうです。

nが小サンプルのケース

さて、データの状況しだいではnを無限とみなしてよいほどnが大きくないということもあるかと思います。

その場合、まず最尤推定量である $\hat{\theta}$ を使うとバイアスが生じます。
また、 $V[\hat{\theta}] - V[\hat{\theta}_2] = \frac{2n+1}{n(n+2)(n+1)^2}>0$ であることを考えると、
分散の観点からも $\hat{\theta}$ よりは、 $\hat{\theta}_2$ を使う方が良さそうです。

より良い推定量を考えるのだから、分散は小さい方が良いですよね。

では、 $\hat{\theta}_3$ はどうでしょう。
$V[\hat{\theta}_2]$ と $V[\hat{\theta}_3]$ を比較してみると、一概にどちらが大きいとはいえなさそうですが、
$n>10$ で $\hat{\theta}_2$ の分散が小さくなることが分かります。上記二つはどちらも不偏推定量ですので、

$n≦10$ では $\hat{\theta}_3$
$n>10$ では $\hat{\theta}_2$

を使うのが良さそうです。

まとめ

以上まとめると、次のようになるかと思います。

$\hat{\theta}=max\{x_1,...x_n\}=x_{max}$
n→∞と見なせる場合に使用可能。

$\hat{\theta}_2 = \frac{n+1}{n} x_{max}$
n>10でよりよい推定が可能。 $\hat{\theta}$ よりも分散が小さい。

$\hat{\theta}_3 = \frac{2}{n} \sum_{i=1}^n x_i$
n≦10で $\hat{\theta}_2$ よりも良い推定が可能。

今回、最尤推定量は惨敗でしたね。
最近はビッグデータ時代で兎にも角にも最尤法や最適化で推定を行いがちですが、小サンプル下ではもっと良い推定がたくさんありますよという話でした。

※コメント欄にて、どうやら小サンプルの所で計算間違いがありそうとのことでご指摘いただきました。今しっかり確認出来ない状況なのですが、後日確認し、修正させていただきます。

2019-12-03

統計検定1級対策問題集~一様分布編~

統計検定対策問題統計検定統計統計学基礎

さて、意外と鬼門な一様分布に関する想定問題(及び過去問)とその解答例を作成して見ました。

離散一様分布の期待値、分散、歪度
連続一様分布の期待値、分散、歪度、積率母関数
変数変換後の分布
max,minの分布
最尤推定量
逆関数法
二番目に小さいものの分布

離散一様分布の期待値、分散、歪度

問1 確率変数Xが1,2,...nの値を取る離散一様分布に従うと仮定する。この時期待値、分散、歪度を求めよ

期待値について
$E[X]=\sum_{x=1}^n \frac{x}{n}=\frac{n+1}{2}$

分散について
$E[X^2]=\sum_{x=1}^n \frac{x^2}{n}=\frac{(n+1)(2n+1)}{6}$

$Var[X]=E[X^2]-E[X]^2=\frac{(n+1)(n-1)}{12}$

歪度について
$E[X^3]=\sum_{x=1}^n \frac{X^3}{n}=\frac{n(n+1)^2}{4}$

$E[\frac{(X-\mu)^3}{\sigma^3}]=\frac{1}{\sigma^3}E[X^3-3\mu X^2 +3\mu^2 X -\mu^3]=\frac{1}{\sigma^3} × \frac{(n+1)^2}{4} × 0=0$

連続一様分布の期待値、分散、歪度、積率母関数

問２確率変数Xが区間 $[a,b]$ 上で連続一様分布に従うと仮定する。この時期待値、分散、歪度、積率母関数を求めよ。

期待値について
$E[X]=\int_{a}^b \frac{X}{b-a} dX=\frac{a+b}{2}$

分散について
$E[X^2 ]=\int_{a}^b \frac{X^2}{b-a} dX=\frac{a^2+ab+b^2}{3}$

$Var[X]=E[X^2]-E[X]^2=\frac{(a-b)^2}{12}$

積率母関数について
$E[exp(tX)]=\frac{1}{b-a} \int_{a}^b exp(tX)dX=\frac{exp(tb)-exp(ta)}{t(b-a)}$

歪度について
$E[X^3]=\int_{a}^b \frac{X^3}{b-a}dX=\frac{(b+a)(b^2+a^2)}{4} [tex: E[\frac{(X-\mu)^3}{\sigma^3}]=\frac{1}{\sigma^3}E[X^3-3\mu X^2 +3\mu^2 X -\mu^3]=0$

変数変換後の分布

問３ Xが区間 $[0,1]$ 上で連続一様分布に従うと仮定する。この時 $X^2$ の累積分布、密度関数、期待値、分散を求めよ

密度関数について
$Y=X^2$ とおく。 $f(x)=1,g(y)$ :Yの密度関数とする。

密度関数の定義より
$P(x < X < x+\Delta x)=P(y < Y < y + \Delta y)$

$⇔ f(x) \Delta x = g(y) \Delta y$

よって

$g(y) = f(x) \frac{dx}{dy}=\frac{1}{2} Y^{-\frac{1}{2}}$

累積分布について
$G(y)=\int_0^y \frac{1}{2}Y^{-\frac{1}{2}} dY =y^{\frac{1}{2}}$

期待値について
$E[Y ]=\int_{0}^1 \frac{1}{2} Y^{\frac{1}{2}} dY =\frac{1}{3}$

分散について
$E[Y^2 ] = \int_{0}^1 \frac{1}{2} Y^{\frac{3}{2}} dY = \frac{1}{5}$

$Var[Y ]=\frac{4}{45}$

max,minの分布

問４ U,Vは独立に区間 $[0,1]$ の一様分布に従うと仮定する。X=max(U,V),Y=min(U,V)とおく。X,Yの累積分布、密度関数、同時密度関数、相関係数を求めよ

累積分布について
(X ≤ x)⇔(U ≤ x)&(V ≤ x)

$G_1(x)=Pr(X ≤ x)=Pr(U ≤ x)Pr(V ≤ x)=\int_{0}^x 1 du \int_{0}^x 1 dv=x^2$

$G_2(y)$
$=Pr(Y ≤ y)=1- Pr(Y ≥ y)=1-Pr(y ≤ U)Pr(y ≤ V)$
$=1-\int_{y}^1 1 du \int_{y}^1 1 dv = 1-(1-y)^2$

密度関数について
$g_1(x)=\frac{d G_1(x)}{dx}=2x$

$g_2(y) = \frac{d G_2 (y)}{dy}=2(1-y)$

同時密度関数について
同時分布 $G(x,y)$ を考える。

$G(x,y) = Pr(X ≤ x, Y ≤ y)=Pr(X ≤ x)-Pr(X ≤ x, y ≤ Y)$
$= Pr(X ≤ x)- Pr( y ≤ U, V ≤ x) = x^2 - \int_{y}^x 1 du \int_{y}^x 1 dv = y(2x-y)$

よって同時密度関数 $g(x,y)$ は

$g(x,y)=\frac{G(x,y)}{\partial x \partial y}=2$

相関係数について
$E[X]=\frac{2}{3}$
$E[Y]=\frac{1}{3}$

ここでY≤Xなので

$E[XY]=\int_{0}^1 \int_{0}^x 2xy dy dx =\frac{1}{4}$

$E[X^2]=\frac{1}{2}$
$E[Y^2]=\frac{1}{6}$

よって相関係数 $\rho(x,y)$ は
$\rho(x,y)=\frac{1}{36}$

最尤推定量

これは、前に記事にしたので、以下の記事の一様分布の欄を参照。
www.bananarian.net

逆関数法

問５連続確率変数Zについて累積分布F(z)=P(Z≤z)について、U=F(z)は区間[0,1]の一様分布に従うことを示せ

これは、ほとんど自明ですが

まず、累積分布関数は単調増加関数なので
$F(z)=P(Z ≤ z)=P(F(Z) ≤ F(z))$

つまり $F(z)=u$ とおくと、

$P(F(Z) ≤ u)=u$

これは一様分布の累積分布関数である。

二番目に小さいものの分布

問６ $U_1,U_2,U_3$ が独立に区間[0,1]の一様分布に従うと仮定する。二番目に小さいものをXとおく。この時のXの密度関数を求めよ

統計検定１級の過去問で「最小のものを $X_1$ ,次に小さいものを $X_2$ ,最も大きいものを $X_3$ とした時の密度と期待値を求めよ」が出題されていますが、本記事の問４と問６を組み合わせると解くことが出来ます。

分布関数Gをまず考える。

$G(x)= Pr(X ≤ x)$
ここで $U_1 ≤ U_2 ≤ U_3$ を一旦仮定してやると(X ≤ x)が成り立つに当たって考えられる条件は次の２通り。

$U_1 ≤ x かつ　U_2 ≤ x かつ U_3 ≤ x$ ...①

$U_1 ≤ x かつ U_2 ≤ x かつ U_3 ≥ x$ ...②

ここで、①は $U_1からU_3$ の大小問わず成り立つが、②は大小の組み合わせに応じて3通り存在する。

よって $Pr(X ≤ x)$ は次のようになる。

$G(x)=Pr(X ≤ x)=Pr(U_1 ≤ x)Pr(U_2 ≤ x) Pr(U_3 ≤ x)+3\{(\int_{0}^x 1 du)^2 (\int_{x}^1 1 du)\}$
$= x^3+3 x^2 (1-x)$

以上より密度関数は

$g(x)=\frac{d G(x)}{dx}=6x(1-x)$

2019-09-21

機械学習のバイアスバリアンス分解を古典的な統計学から丁寧に解説する

機械学習統計統計学基礎統計学雑記統計検定ディープラーニング

今日は、バイアスバリアンス分解の話をしようと思います。

はじめに
古典的な統計学におけるMSEの分解
- パラメータの推定
- MSEの分解
機械学習におけるバイアスバリアンス分解
まとめ
double descent

はじめに

「バイアスバリアンス分解」と言う言葉を聞いて、皆さんは何を思い浮かべるでしょうか。

恐らく古典的な統計学だけを勉強した人であれば、パラメータの推定量に関するMSEの分解を思い浮かべるかもしれません。

一方機械学習から勉強した人であれば、トレーニングデータで評価した予測値と予測変数(目的変数)とのMSEを思い浮かべることでしょう。

実は、後者のことを一般にバイアスバリアンス分解と呼び、前者はMSEをバイアスとバリアンスに分解できるだけの話で特に名前は与えられていない(と思われる)のですが、当然、両方MSEですし、やってることは同じなので、両方勉強している人であれば、そんなもの分けて議論することがナンセンスだともなることでしょう。

しかし、恐らく機械学習に馴染みのない方からすると、機械学習サイドでよく言われるところのバイアスバリアンス分解はわかりにくく、

躓きポイントじゃないかな〜と私は考えています。

そこで今回の記事では、まずはじめに古典的な統計学で語られるところのMSEの分解とその理解について説明した上で、それとのモチベーションの違いに触れながら機械学習で語られるところのバイアスバリアンス分解について説明しようと思います。

そして、最後にオマケとして、最近話題のdouble descentの話にも触れようと思います。

古典的な統計学におけるMSEの分解

パラメータの推定

まずはじめに、古典的な統計学の目的を確認しておきましょう。
統計学では、あるデータXの性質を'分布'という観点から考察したいというモチベーションから出発します。
ここであるデータXの性質とは「期待値や分散の値」だけでなく、「Xに影響を与える変数は何か」なども含みます。

とにかく、データXの分布(当然この分布は仮定する)のパラメータを推定し、それを解釈、利用することでXの性質を述べたり、その述べた性質が妥当か検証したり、パラメータの値を使って予測したりといったことを行いたいわけです。

なので、とにもかくにも仮定した分布のパラメータの値を、持っているデータを使ってうまく当てたいわけです。

ここでXに対して、何かしらのパラメータ $\theta$ によって決まる分布を仮定します。

パラメータ $\theta$ の値は神様しかわからないわけですが、持っているデータXを用いて得られる推定量 $g(X)$ で何とか誤差が小さくなるようにパラメータの値を推定したい！と言う話になってきます。

そこで、次のような評価方法Lを考えることになります。

$L=(\theta-g(X))^2$

このLは神様しか知らない真の値と推定に用いる値との二乗誤差です。

Xは確率変数を仮定しており、確率的にばらつくはずなので、このLをXの期待値で評価すると

$MSE=E[(\theta-g(X))^2]$

これが、MSE (Mean Squared Error) です。

と言うわけで古典的な統計学でMSEを用いる基本的なモチベーションとしては、

あるデータXの性質をパラメータに置き換えて述べたい
しかし、パラメータの値はわからないのでうまい推定量で推定したい
うまい推定量とは何かを検討するにあたり、MSEで評価したい

ということになります。

MSEの分解

さて、MSEですが、推定量のバイアスと分散に分解することができます。

ちなみにここでいう推定量のバイアスとは何を指すかというと、
推定したい真の値と、推定量の期待値の差です。
つまり $\theta-E[g(X)]$ がバイアスです。

それを示していきましょう。

$MSE$

$=E[(\theta-g(X))^2]$

$=E[(\theta-E[g(X)]+E[g(X)]-g(X))^2]$

$=E[(\theta-E[g(X)])^2+(g(X)-E[g(X)])^2+2(\theta-E[g(X)])(E[g(X)]-g(X))]$

ここで、最後の項は期待値をとるとゼロになりますね。
$E[2(\theta-E[g(X)])(E[g(X)]-g(X))]=0$
Xの期待値を考えていることに注目してください。

よって、
$MSE$

$=E[(\theta-E[g(X)])^2]+E[(g(X)-E[g(X)])^2]$

$=(\theta-E[g(X)])^2+V[g(X)]$

１項目はバイアスの二乗ですね。
２項目は推定量の分散です。

こうしてMSEをバイアスと分散に分解することができました。

ここまで見れば、例えば不偏推定量(バイアスが０の推定量)を考えるモチベーションもわかるかと思います。
MSEを小さくしたいが、MSEは正確には分散とバイアスの二つの要因から成る。
バイアスが0の推定量のクラスであればとにかく分散を小さくすることだけに注力すればいいので、考えやすいですね。

当然、モチベーション次第では、若干バイアスが生じてもそれ以上に分散を小さくできるような推定量があれば、不偏である必要もないですね。
そんなことを考えることができるわけです。

機械学習におけるバイアスバリアンス分解

予測と分類

一般に機械学習と言われる手法では、上記の話とはモチベーションが少しだけ違います。
パラメータには重きを置きません。

とにかく目的変数を予測(分類)したい！！とにかく正確に予測したい！！というモチベーションで話が進みます。

いやいや、先ほどの統計学の説明でも予測はモチベーションの一つとしてあったじゃないかというのは当然の指摘ですが、先ほどの古典的な統計学のモチベーションのメインはとにかくパラメータにあります。その成果を使って予測を行うこともあるというだけです。

古典的な統計学の文脈で良しとされる手法は、予測変数のMSEを出来うる限り小さくしているわけではないのです。

一方、機械学習でよしとされる手法は、パラメータのMSEを最小化しているわけではないのです。

若干ミスリードの感があるので補足しておくと、これは、統計学は予測に興味がないといっている訳ではありません。統計学は予測に"も"興味はあるが、他のことにも興味があるため、パラメータ目線で話を進めることが多かったというだけです。そのため機械学習を予測に特化した統計学として解釈することも当然可能です。

では、機械学習の話に戻ります。
例えばデータXを使ってデータYの値を予測したいといった状況を考えることにしましょう。

この時、得られたデータ $X^{(T)}$ を使って $g(X)$ というYを予測するための関数(予測器)を考えて、新しいXが与えられた時の'未知のY'を予測したいというのが機械学習のモチベーションです。

ここで、添字 $(T)$ を与えてモデルを作る時に使っている $X^{(T)}$ と、予測したいYに対応するXを区別していることに注意してください。

既に持っているデータは、答えであるYがもう分かっているのだから、わざわざ予測する必要なんてありません。
まだわからないものを予測してこその予測器です。

オーバーフィット

既に持っているデータ $X^{(T)}$ とそれに対応する答え $Y^{(T)}$ の組をテストデータと呼びます。

実はテストデータの二乗誤差

$\sum_{i=1}^{N}(Y_i^{(T)}-g(X_i^{(T)}))^2$

を0にする方法は簡単で、パラメータの次元を増やす(モデルを複雑にする)ことで誤差0の理想のモデル(!?)が完成します。

しかし、これには問題があって

このようにテストデータに対する誤差最小化モデルは、新しいデータに対する予測性能がありません。

これは当然で、テストデータの二乗誤差を最小化するモデルを作っただけで、新しいデータに対する二乗誤差最小化モデルを考えていないからです。

このようなテストデータにだけ過剰にフィットしている状況をオーバーフィットと呼んだりします。

よくこんな図とともに説明されます。
f:id:bananarian:20190917234218j:plain *1
点線がテストデータで、横軸のcomplexityは複雑度ですが、とりあえずパラメータの次元だと思ってください。縦軸もRiskですが、今回の話で言うところの二乗誤差だと思ってもらえれば問題ないです。

実線が新しいデータで見たものです。複雑度が増すと一定のラインを超えると、Riskが猛烈な勢いで増えています。

じゃあどうすればオーバーフィッティングを回避できるのかというと、評価の基準をテストデータに関する二乗誤差ではなく、新しいデータに関する二乗誤差にしてやれば良い訳です。

バイアスバリアンス分解

それでは、MSEを考えましょう。改めて文字を整理します。

予測器を作るにあたって用いるテストデータセットをDとして次のように書くことにします。

$D= \{(X_1^{(T)},Y_1^{(T)}),...,(X_N^{(T)},Y_N^{(T)})\}$

このテストデータセットDを用いて作った予測器gを次のように書くことにします。

$g(x;X^{(T)},Y^{(T)})=g(x;D)$

これは、gは新しい予測のための変数xを入れることでxに対応する予測値を出力する関数であるが、
作るにあたって得られたテストデータDの値を用いている(つまり、Dを条件付けたもとで得ている)ことを明示的に示すために、このような書き方をしています。

さて、テストデータとは区別して、データXが与えられた時に生じるであろう値Yをg(x)を使って予測したいので、次のような誤差Lを持って評価すれば良いことがわかりますね。

$L=(Y-g(X;D))^2$

ただしここで、Yは確率変数です。そのため、期待値をとって次のように考えてやる必要があります。

$E_{y|x}[(Y-g(X;D))^2]$

演算子 $E_{y|x}$ はxを条件づけた時のyの期待値を求めていることを表します。ちゃんと積分の形に直すと次のようなことを言っているのと同じです。

$\int (Y-g(X;D))^2 p(y|x)dy$

さて、実はもう一つ確率変数があります。Yを予測するにあたって、Xは固定なのですが、データセットDは固定していません。

モデルを考えるにあたって確率変数の組であるトレーニングセットを使って予測器を作っているため、この $g(X;D)$ は確率変数です。
要はXを固定しているだけで $X^{(T)}$ を固定した訳じゃないのです。よって、こいつの期待値も取る必要がありますね。

$MSE=E_{D|x}[E_{y|x}[(Y-g(X;D))^2]]$

これが機械学習で考えているところのMSEです。

式は一緒でも、考えているモチベーションが若干異なるため、このような形になる訳ですね。
さて、分解してみましょう。手順は先ほどと同様です。

$MSE$

$=E_{D|x}[E_{y|x}[(Y-g(X;D))^2]]$

$=E_{D|x}[ E_{y|x}[ (Y-E_{y|x}[Y]+E_{y|x}[Y]-g(X;D) )^2 ] ]$

$=E_{D|x}[ E_{y|x}[ (Y-E_{y|x}[Y])^2 + (g(X;D)-E_{y|x}[Y])^2 ] ]$

若干複雑ですが、中身をしっかりみてください。右の項はYの値によらずバイアスの二乗の形であることがわかりますね。

さて、続きの計算を行いましょう。
$=E_{D|x}[ V_{y|x}[Y]+(g(X;D)-E_{y|x}[Y])^2]$

$= V_{y|x}[Y]+E_{D|x}[(g(X;D)-E_{y|x}[Y])^2]$

さて、右側の項をさらに分解しましょう。分解方法は今までの手順と同じです。

$=V_{y|x}[Y]+E_{D|x}[(g(X;D)-E_{y|x}[Y])^2]$

$= V_{y|x}[Y] + E_{D|x}[(g(X;D)-E_{D|x}[g(X;D)]+E_{D|x}[g(X;D)] -E_{y|x}[Y])^2]$

$=V_{y|x}[Y] + E_{D|x}[(g(X;D)-E_{D|x}[g(X;D)])^2+(E_{D|x}[g(X;D)]-E_{y|x}[Y])^2]$

$=V_{y|x}[Y] +V_{D|x}[g(X;D)]+(E_{D|x}[g(X;D)]-E_{y|x}[Y])^2$

機械学習モデルにおいては、しばしばこの式の第三項目をバイアスと呼んだりします。
第二項目はバリアンス、そして第一項目はトレーニングデータセットと無関係なデータの分散です。

この第一項目についてはトレーニングデータではどうしようもないので、いじれません。

よって、予測においてはトレーニングデータとモデルを使っていじることが可能な第二、第三項が重要になります。

まとめ

結局、やっていることも、評価方法も同じであるということが分かっていただけたと思います。

分析において評価したいものが違っていて、それは手法の背後にあるモチベーションからこそなのです。

そして、評価したいものが違うがゆえに、期待値演算が二重になって若干複雑になっているのが機械学習(or 最新の統計学)におけるバイアスバリアンス分解なのです。

double descent

先ほど、オーバーフィッティングで次のような図を確認しました。
f:id:bananarian:20190917234218j:plain

しかし、2018年、実は次のような状況にあることが[Belkin, Mikhail, et al., 2018]によって確認されました。
f:id:bananarian:20190917234853j:plain *2

実は、もっとモデルの複雑度を増やすとまたRiskが減少し始めるらしい。これをdouble descentと呼んだりします。
これは、古典的な統計学の手法では、そもそもこの先がわかると言う状況がなかったので知られていなかったのですが、例えばニューラルネットワークなどを考えるとこのような状況が発生します。

この仕組みはまだ完全に判明していないので、理屈を説明することが出来ないのですが、

これは今注目を浴びているトピックで、テストデータでRMSEの谷を見るだけではどうやらダメらしいということで研究が進んでいる最中だったりします。

*1:Belkin, Mikhail, et al. "Reconciling modern machine learning and the bias-variance trade-off." arXiv preprint arXiv:1812.11118 (2018)より

*2:Belkin, Mikhail, et al. "Reconciling modern machine learning and the bias-variance trade-off." arXiv preprint arXiv:1812.11118 (2018)より

2019-07-13

統計検定1級対策解説~MCMC編1~

統計検定統計検定対策問題統計統計学基礎ベイズ統計学 MCMC R

統計検定1級の範囲表によると、いつのまにやらMCMCやベイズが追加されていますね。

しかし、公式教科書は2013年出版の物しか出ていないため、その辺の範囲には対応していません。そこで、一通り何個かの記事を通して解説を行い、関連問題も上げていこうと思います。

MCMCとは何か
マルコフ連鎖

MCMCとは何か

MCMCは、マルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo method)の略です。

ベイズ統計では、複雑な事後分布の期待値や分散を考える必要が出てきます。
その際、そうした期待値や分散は簡単に計算することが出来ない場合が多くあります。

そうした場合に対応するために、事後分布を定常分布とするマルコフ連鎖を構成することで、近似的な事後分布からのサンプルを得ることで、期待値や分散を計算しようというわけです。

そういった手法群を総称してMCMCと呼びます。

MCMCを考えるには、大雑把に言って次の3つを考える必要があるかと思います。

・マルコフ連鎖はどのような条件の下で収束するのか(定常分布が存在するのか)

・得たい事後分布を定常分布とするマルコフ連鎖をどうやって見つけるのか

・効率の良いアルゴリズムはどのようなものを考えればよいのか

今回は主に上二つに焦点を当てて説明していきます。

マルコフ連鎖

マルコフ連鎖と定常分布

一つ前の状態によって分布が変化する確率過程をマルコフ連鎖と呼びます。

例えば、次のような例を見てみます。

f:id:bananarian:20190712161745p:plain

この表は、状態A~Cから状態A~Cへの遷移行列です。
読み方は、例えば2行1列目には0.1が入っていますが、これは「状態Bから状態Aへ行く確率が0.1である」ことを表します。

ここで、経験的に初期状態で状態Aになる確率が0.3、状態Bになる確率が0.4,状態Cになる確率が0.3であったとしましょう。

この時、先ほどの遷移行列から考えれば

初期状態がAで、Bに遷移する確率は0.3×0.5=0.15であることがわかりますね。

同様にして初期状態BからBに遷移する確率は0.12

初期状態CからBへ遷移する確率は0.03

つまり、2時点目で状態Bである確率は0.3ということになります。

この計算を、全てのパターンで10回やってみることにします。ちょっと大変なので、Rに計算を投げてやると次の通り。

> XX=matrix(c(0.2,0.1,0.6,0.5,0.3,0.1,0.3,0.6,0.3),3,3)
> initV=c(0.3,0.4,0.3)
> res=initV
> 
> for(i in 1:10){
+ res=res%*%XX
+ print(round(res,digit=3))
+ }
     [,1] [,2] [,3]
[1,] 0.28  0.3 0.42
      [,1]  [,2] [,3]
[1,] 0.338 0.272 0.39
      [,1] [,2]  [,3]
[1,] 0.329 0.29 0.382
      [,1]  [,2]  [,3]
[1,] 0.324 0.289 0.387
      [,1]  [,2]  [,3]
[1,] 0.326 0.287 0.387
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386
      [,1]  [,2]  [,3]
[1,] 0.326 0.288 0.386

結果を見てみるとわかる通り、あるところから、確率が0.326 0.288 0.386で動かなくなっています。

この状態を、マルコフ連鎖が収束したと呼び、この時の定常分布は0.326 0.288 0.386です。

定常分布への収束条件

マルコフ連鎖は必ず定常分布へ収束するわけではありません。
次の3条件を満たしている必要があります。

・既約性
ある状態からある状態へ、有限回の遷移でたどり着けるという条件が既約性です。

・正再帰性
任意の要素が、有限期間内で何度も発生しうることを正再帰性と呼びます。
つまり、正再帰性を満たさないとは、70回目の遷移以降、状態Aは発生しなくなるなどと言った状況です。

・非周期性
遷移に周期が存在しないことを非周期性と呼びます。

先ほど用意した遷移行列は、この3条件をしっかり満たしているため、定常分布へ収束したというわけです。

詳細つり合い条件

定常分布への収束条件はわかりました。
しかし、私たちがやりたいことは、ある特定の定常分布へ収束するようなマルコフ連鎖を構成することです。

つまり、ある特定の定常分布 $p(x)$ へ収束するような遷移行列 $C$ を考えたいということになります。

そこで登場するのが詳細つり合い条件です。
ある状態iから状態jへ遷移する確率を $C(j|i)$ とおくと、詳細つり合い条件は次のようになります。

$C(i|j)p(x=j)=C(j|i)p(x=i)$

この条件が全ての $i,j$ で成り立つとき、この遷移行列を用いた、マルコフ連鎖は定常分布 $p(x)$ へ収束することが知られています。

つまり、私たちはこれで、欲しい事後分布を定常分布とするようなマルコフ連鎖を考えることが出来るようになったというわけです。

これらの条件を踏まえると、次に考えるべきは、

じゃあ、事後分布からのサンプルを取ってこれるような効率的なアルゴリズムって何？ってことですよね。

それについてはまた次回。

2019-07-12

統計検定1級対策解説~一様最強力検定編~

統計検定統計検定対策問題統計統計学基礎

今回は、一様最強力検定の話をします。

仮説検定
- 帰無仮説と対立仮説
- 誤差の確率
一様最強力検定
- 単純仮説
- 証明

仮説検定

とりあえず、ザックリと仮説検定の話をしていきます。

仮説検定とは、「ある仮説を否定(棄却)出来るかどうかを検討したい」だとか、「ある仮説に基づいて、何かしらの意思決定を行いたい(モデル選択等)」といった場合に利用する、判定方法です。

基本的な枠組みは次の通り。

帰無仮説と対立仮説

まず、棄却出来るかどうかを考えたい仮説として帰無仮説 $H_0$ を考えます。
例えば、ある正規分布に従っているデータを考えた時に、恐らく期待値は0ではないはずだと考えたとします。しかし、これは予想であって、客観的な根拠がありません。そこで、帰無仮説として期待値が0であるという仮説を設けます。この仮説を、仮説検定を使って棄却出来れば、「このデータの期待値は0ではない」ということが出来ますね。

つまり
帰無仮説 $H_0 : 期待値が0$
対立仮説 $H_1 : 期待値は0ではない$

というような2種類の仮説をまず、設定するということですね。

誤差の確率

検定をするにあたっては、帰無仮説が正しいのか、対立仮説が正しいのか分からない状態で意思決定を行います。つまり、帰無仮説が正しい状況と対立仮説が正しい状況の両方を考えた上で、妥当な意思決定を検討しなければなりません。

そこで気になるのは次のような数値でしょう。

・帰無仮説が正しい時に棄却される確率 (第一種の誤差)
・対立仮説が正しい時に受容される確率 (第二種の誤差)

両方の確率が0、つまり確率1で正しい仮説を選べるような検定なら最高ですが、少し考えるとわかりますが、そのような状況はほとんど発生しません。どうしてもミスをする確率は出て来てしまいます。

そこで、とりあえず第一種の誤差は小さい値に抑え込もうという方針で妥協します。

そこで利用するのが有意水準です。

つまり、有意水準5%に決めるとは、(帰無仮説内の任意のパラメータにおいて)第一種の誤差を5%以下に抑え込むことを指します。

数学的に書くと、

まず、考える1変量パラメータをもつ確率分布を $P_{\theta}$ とおく。
パラメータ空間全体を $\Theta$ 、帰無仮説におけるパラメータの集合を $\omega$ と置くと、対立仮説におけるパラメータの集合は $\Theta-\omega$ となりますが、

棄却域を $R$ とすると、第一種の誤差は $P_{\theta}(R),ただし\theta \in \omega$

よって、有意水準 $\alpha$ に決めるとは

$P_{\theta}(R)≦\alpha$ for all $\theta \in \omega$

当然、仮説検定の性能を考えるのであれば、

対立仮説が正しい場合に、きちんと帰無仮説を棄却する確率も知りたいですよね。これを検出力(Power)と呼びます。

つまり、

$P_{\theta}(R)$ ただし $\theta \in \Theta-\omega$

を検出力と呼ぶわけです。

一様最強力検定

ここまでが前段で、ここからが本題です。
で、有意水準を使って第一種の誤差は固定しました。

よって、第一種の誤差はとりあえず置いておいて、第二種の誤差をできうる限り小さくしたい(つまり、検出力を大きくしたい)ということに注力することにしましょう。

そこで真っ先に思い浮かぶのが、「対立仮説が正しい場合のパラメータ空間内で一様に検出力が最大になるような検定問題を考えることが出来ないか？」ということです。

つまり、検定の方法は色々考えられるわけですけど、パラメータがどんな値を取ったとしても(対立仮説が正しいもとで)、検出力が最大になる最強の検定があれば、それを使えばいいですよね？

結論から言うと、普通はそんな検定はありません。統計検定1級公式教科書でも、「一般に、そのような検定は存在しない」と書かれています。

ただし、「一般に」成り立たないだけで、特殊なケースでは一様最強力検定が存在します。今回はそれを示すことにします。

単純仮説

単純仮説とは、パラメータ空間全体 $\Theta$ が2種類の要素しか持たず、帰無仮説 $\omega=\{\theta_1\}$ 、対立仮説 $\Theta-\omega=\{\theta_2\}$ のような仮説を単純仮説と呼びます。

この時実は、棄却域として $R^*=\{x : \frac{p_{\theta_2}(x)}{p_{\theta_1}(x)}≧k\}$ を設定することで、これは最強力検定になります。

証明

ここで、任意の標本空間の領域を $R$ ,有意水準を $\alpha$ と置くことにします。

$P_{\theta_2}(R^*)-P_{\theta_2}(R)=\int_{\{R^* \cap R^c\}} p_{\theta_2}(x)dx - \int_{\{R^{*c} \cap R\}} p_{\theta_2}(x) dx$

これは、 $R^*=(R^* \cap R) + (R^* \cap R^c)$ なのですぐわかります。

ここで、 $R^{*} \cap R^c$ は $R^*$ の部分集合なので、 $p_{\theta_2}(x)≧k p_{\theta_1}(x)$ が成り立ちます。

よって、 $\int_{\{R^* \cap R^c\}} p_{\theta_2}(x)dx ≧k \int_{\{R^* \cap R^c\}} p_{\theta_1}(x)dx$

同様に考えると

$\int_{\{R^{*c} \cap R\}} p_{\theta_2}(x) dx ＜k \int_{\{R^{*c} \cap R\}} p_{\theta_1}(x) dx$

よって

$P_{\theta_2}(R^*)-P_{\theta_2}(R)≧k(P_{\theta_1}(R^*)-P_{\theta_1}(R))$

ここで、有意水準で $\alpha$ で抑えていることから

$P_{\theta_2}(R^*)-P_{\theta_2}(R)≧k(P_{\theta_1}(R^*)-P_{\theta_1}(R))≧0$

よって

$P_{\theta_2}(R^*)≧P_{\theta_2}(R)$

以上より最強力検定であることが示せた。

2019-05-15

統計検定1級対策解説~情報量基準編~

統計検定対策問題統計検定統計統計学基礎統計学雑記機械学習 R

今回はAICです。
AIC周りの話は、カルバックライブラー情報量とAICにどういう結びつきがあるのか概要を知っておくことが重要になってきます。

AICとは
何の役に立つんだって話
カルバックライブラー情報量
- カルバックライブラー情報量とAICの関係性

AICとは

AICは次のような式で表されます。

$AIC=-2 \sum_{i=1}^n log f(X_i,\hat{\theta})+2 dim(\theta)$

まず $\hat{\theta}$ は、n個のサンプルから得られるパラメータの最尤推定量です。そして、 $dim(\theta)$ は分布の自由度となっています。

このAICとはなんぞやというのを理解することがこの章で重要な話になってきます。

何の役に立つんだって話

このAIC、何の役に立つのかというと、オーバーフィットを回避するためのモデル選択に利用できます。

オーバーフィットとは

例えば、こんな感じのデータがあったとします。
f:id:bananarian:20190515123427p:plain

x=seq(0,1,0.1)
train=sin(2*pi*x)+rnorm(11,sd=0.3)

で、どういう関数を当てはめたら良いかよくわからないので、最小二乗法(又は最尤法)を使って、多項式を当てはめてみることにするとします。

多項式っていうのはこういうやつですね。

$f(x)=\sum_{m=1}^M \alpha_{m-1}X^{m-1}+\epsilon$

これはM-1次の多項式です。
Mはどうやって決めるのか問題がありますが、最尤法であれば、対数尤度最大、最小二乗法であれば二乗誤差を最小にするように、Mを適当に動かしてやると探すことが出来ますね。ちなみに、機械学習界隈ではよくＲＭＳＥが利用されますがどれを使おうとさしたる違いはありません。

で、普通に最小二乗法、又は最尤法を使うと次のような曲線が得られます。
f:id:bananarian:20190515124216p:plain

「おぉ～～全ての点を通ってて良い感じ～！」

と満足してはいけません。この曲線には問題があります。

同じ分布から新しいサンプルをこんな風に取ってきたとします。

xtest=seq(0,1,0.01)
test1=sin(2*pi*xtest)+rnorm(length(xtest),sd=0.3)

今の曲線をこの新しいサンプルにあてはめるとこんな感じ。
f:id:bananarian:20190515124809p:plain

変なところにすっ飛んでいってる所があり、予測が見当違いになってしまいます。

このように、既存のデータに極端にフィットしてしまい、新しいデータに対する予測性能が無くなってしまう現象をオーバーフィッティングと呼びます。

パラメータの次元(今回ならM)も最適化しようとした場合に、最尤法や最小二乗法を用いると、このような状況が発生します。それを回避するために利用される手法がAICです。

AICの式の意味

AICの式をもう一度確認してみます。

$AIC=-2 \sum_{i=1}^n log f(X_i,\hat{\theta})+2 dim(\theta)$

「このAICを最小化するようなパラメータを選択するのが良い」というのがAICが利用される場合の主張です。

第1項は最大対数尤度ですね。普通の最尤法ではこのAICの第一項を最小化(マイナスがついているので)するような値を選択しているわけですが、それではパラメータの次元を弄り始めるとさっきのような問題が発生します。

そこでAICでは第二項目としてパラメータの自由度を加えます。

つまり、パラメータの自由度を増やせば増やすほど罰則がかかるようになっているわけですね。

AICを使ってみる

では、自由度をいじりながらAICの値を確認してみます。簡単にプロットしてやると次のようになりました。

f:id:bananarian:20190515130147p:plain

最小点は6ですね。つまり、「10次多項式なんていらない、6次で十分だ」とAICで主張できるわけです。

ちなみに6次でさっきのプロットをしてみるとこんな感じ。
f:id:bananarian:20190515130703p:plain
f:id:bananarian:20190515130715p:plain

見た目的にも最小二乗法を使った場合よりも良いですね。

余談(機械学習でよく使われる手法の話)

ちなみに、機械学習や予測の文脈では次のような手法も使われます。

持っているデータを

・最小二乗法(最尤推定)用のトレーニングデータ
・次元Mを決めるためのテストデータ

にわけることで、先ほどの問題を解決するといった手法です。

こんな感じでテストデータとトレーニングデータがあったとします。

x=seq(0,1,0.1)
train=sin(2*pi*x)+rnorm(11,sd=0.3)
test=sin(2*pi*x)+rnorm(11,sd=0.3)

この２つのデータに対して、トレーニングデータで学習した最尤推定量を元に、パラメータの自由度をいじって、その際のRMSEを見ることを考えます。

すると次のような推移になります。
f:id:bananarian:20190515131537p:plain

テストデータを用いた場合は、途中でRMSEの減少が止まっています。
つまり、テストデータを用いた場合、10次元も6次元もいらない。3次多項式で十分だと主張できることになります。

先ほどのAICは、尤度を利用していることからも分かる通り、サンプルサイズが少ない場合十分に機能しない可能性があります。そのため、データセットをこのように分ける手法が取られることもあるわけです。

カルバックライブラー情報量

カルバックライブラー情報量は、２つの分布間の(ある種の)距離を測る考え方です。
次のように定式化します。

２つの密度関数f,gについて

$KL(f|g)=\int f(x)log(\frac{f(x)}{g(x)})dx$

２つの密度関数の対数の差の期待値を取っているのが見えますか？
もし、二つの密度が全く同じであれば0になります。

カルバックライブラー情報量とAICの関係性

ある分布 $f(y;\theta)$ から標本Xが得られたとします。
この標本を使って分布のパラメータ $\theta$ を推定することを考える場合、その推定量 $\hat{\theta}$ を使って尤度 $f(y;\hat{\theta})$ を考えた場合、真の分布 $f(y;\theta)$ に近くならなければならないはずです。

ここで、未知の標本(架空の物、得られていない) $Y$ を使って真の分布と推定量を突っ込んだ分布の距離をカルバックライブラー情報量で考えてやると次のようになります。

$KL(f(y; \theta)|f(y;\hat{\theta}))=\int f(y;\theta)log(f(y; \theta))dy -\int f(y; \theta) log(f(y; \hat{\theta)}) dy$

この値を小さくするためには、第1項目は弄りようが無いので、第2項目を考える必要が出てきますね。
つまり、第二項目の

$\int f(y; \theta) log(f(y; \hat{\theta)}) dy$

を最大にするような推定量が汎化性能(予測とかの意味で良い性質を持つこと)の意味で良い推定量であると言えます。
これを平均対数尤度と呼んだりします。

もし、ここで新しいデータとして今の分布と同様の分布から $Z$ が得られたとすると、

平均対数尤度の推定値は $\frac{1}{n}\sum log(f(z; \hat{\theta}))$ となるわけですが、 $\hat{\theta}$ はデータXの下での最尤推定量だったので、

$\frac{1}{n}\sum log(f(X; \hat{\theta)}) ≧ \frac{1}{n}\sum log(f(z; \hat{\theta)})$

が成り立ちます。つまり、よりよい推定量を考えるためには平均対数尤度を最大化しなければならないのに、それよりも値が大きいXの対数尤度を最大化して満足してしまっていることがわかります。

しかし、そうはいっても「もしZを取ってこれたら」という話をしているに過ぎないので、実際はデータXでなんとかしなければなりません。

そこで、うまく平均対数尤度に近づくように補正する項を持ってきてやろうといことで漸近的に出てくるのが、先ほどAICで出てきた $dim(\theta)$ です。

こうして、平均対数尤度を考える代わりに、既存のデータでどうにかできるようAICが生まれてきたというわけですね。

2019-05-10

統計検定1級対策問題集~フィッシャー情報量編~

統計検定対策問題統計検定統計統計学基礎

統計検定1級では、フィッシャー情報量を求める問題が頻出しています。そこで、フィッシャー情報量についてまとめました。

フィッシャー情報量とは
- フィッシャー情報量の第二の導出
フィッシャー情報量（行列）の具体例
- ベルヌーイ分布
- 正規分布

フィッシャー情報量とは

尤度関数 $p(x;\theta)$ を考えます。
この時、この尤度関数について対数を取った対数尤度をパラメータ(ベクトル)について１階微分し、二乗して期待値を取ったものをフィッシャー情報量(行列) $I(\theta)$ と呼びます。つまり

$I(\theta)= E[\{\frac{\partial}{\partial \theta} log( p(x;\theta) ) \}^2]$

フィッシャー情報量(行列)は様々な場面で利用されますが、もっとも有名なのはクラメルラオの下限です。

あるパラメータ(ベクトル)における不偏推定量の分散の下限はフィッシャー情報量(行列)の逆数(逆行列)に等しくなります。これをクラメルラオの下限と呼びます。

不偏推定量は、バイアス0の推定量でしたので、不偏推定量のクラスで推定量を考える場合、MSEを最小にするような推定量を考えるには分散を小さくすることだけを考えれば良いので、その下限がわかるのは非常に強力です。

フィッシャー情報量の第二の導出

フィッシャー情報量は定義通り１階微分の二乗の期待値でも求めることは出来ますが、次の計算も適当な正則条件の下で同値になります。

$I(\theta)=E[-\frac{\partial^2}{\partial \theta^2} log(p(x;\theta))]$

場合によっては二乗するより二回微分した方が簡単になることもあるので、この関係を知っておくことは重要です。一応簡単に証明しておきます。

$\frac{\partial^2}{\partial \theta^2} log(p(x;\theta))=\frac{\partial}{\partial \theta} \frac{1}{p}\frac{\partial p}{\partial \theta}$

$=-\frac{1}{p^2} (\frac{\partial p}{\partial \theta})^2 + \frac{1}{p} \frac{\partial^2 p}{\partial \theta^2}$

$=-(\frac{\partial logp}{\partial \theta})^2 +\frac{1}{p} \frac{\partial^2 p}{\partial \theta^2}$

ここで $E[ \frac{1}{p} \frac{\partial^2 p}{\partial \theta^2} ]=0$ …(微分と積分の交換が成り立てば)

よって $E[- \frac{\partial^2}{\partial \theta^2} log(p(x;\theta)) ]=I(\theta)$ であることがわかります。

だんだんこの公式で慣れてくると、

「んんん？？二乗したやつにマイナスつけるんだっけ？二回微分した奴にマイナスつけるんだっけ？？」と混乱してくるかもしれませんが、その時はクラメルラオの下限を思い出すと良いかと思います。

分散の下限になるってことは、値は正になるはずってのが感覚的に普通ですよね。だから二乗したやつにマイナスがつくことなんてあり得ません。

フィッシャー情報量（行列）の具体例

具体的な計算が無いとよくわからないと思うので、１変量の場合と多変量の場合をそれぞれ具体例で確認してみようと思います。

ベルヌーイ分布

独立同一にベルヌーイ分布に従うサンプルを考えた場合、その同時尤度は次のよう。

$p^{\sum x_i}(1-p)^{n-\sum x_i}$

そこで、対数を取ってやると

$log(p) \sum x_i + log(1-p) ( n-\sum x_i)$

これを $p$ について二回微分すると次のようになりますね。

$-\frac{\sum x_i}{p^2}-\frac{n-\sum x_i}{(1-p)^2}$

ここで、ベルヌーイ分布の期待値は $p$ なので、期待値を取ってマイナスをとるとフィッシャー情報量 $I(p)$ は

$I(p)=\frac{n}{p}+\frac{n}{1-p}=\frac{n}{p(1-p)}$

正規分布

パラメータが２つあるので、フィッシャー情報行列になります。

$\theta=(\mu,\sigma^2)$ とおく。

ここで、同時尤度は次のようになります。

$(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum (x_i-\mu)^2}{2\sigma^2})$

ここで、対数を取ると

$\frac{-n}{2} log(2\pi \sigma^2)-\frac{\sum (x_i-\mu)^2}{2\sigma^2}$

さて、ここで $\mu$ に関して２回微分してやると

$\frac{-n}{\sigma^2}$

これについてマイナスをつけて期待値を取ると

$\frac{n}{\sigma^2}$

また $\sigma^2$ に関して２階微分してやると

$\frac{n}{2 \sigma^4}-\frac{\sum (x_i-\mu)^2}{\sigma^6}$

これについてマイナスをつけて期待値を取ると

$\frac{n}{2 \sigma^4}$

$\mu$ , $\sigma^2$ で１回ずつ微分してやると

$\frac{-\sum (x_i-\mu)}{\sigma^4}$

これについてマイナスをつけて期待値を取ると

$0$

以上より、フィッシャー情報行列 $I(\theta)$ は次のようになります。

$I(\theta)=\begin{pmatrix} \frac{n}{\sigma^2} & 0 \\ 0 & \frac{n}{2 \sigma^4} \\ \end{pmatrix}$

ポアソン分布の場合や指数分布の場合も練習になるのでやってみると良いかと思います。

2019-05-03

十分統計量に関する小話

統計学雑記統計統計検定

ここ最近、統計検定関連の記事が続いていますね。

統計検定の範囲の中に「十分統計量」という単元があります。

前に２記事くらい練習問題記事を書きましたが、今回は十分統計量についてツラツラ思っていることを書いていきます。

完全に個人の感想です。

そもそも統計検定の範囲になっているが

統計検定に関する勉強をしている受験生の方々は、当然今の時代に、統計学を必要としている方々のはずなので、「この十分統計量って何？？？」と立ち止まってしまったり、「わざわざこんなこと考える必要ある？？」なんて感じる人も少なくないかと思います。

「何やら、十分統計量だなんだと難しい計算をわざわざゴリゴリやっていて、ラオブラックウェルの定理でMSEを小さくできる推定量を探せる！嬉しい！」と言われても(これが、統計検定１級の該当範囲)、

「うーん、そんな難しいことをせずとも、乱数を発生させて想定される推定量を比較してやれば良いんじゃないか」

と思う方も多いと思います。

そもそも、十分統計量の前に、不偏推定量についてもなんの意味があるんだろうと素朴な疑問を覚えた方も多いのではないでしょうか。

「わざわざUnbiassedな推定量を探して、そのクラスで分散小さい選手権を開催してる暇があったら、UnbiassedだろうがBiassedだろうが何でもいいから、想定される推定量を片っ端から集めてきて、乱数を100000個くらい発生させて、empiricalなバイアスと分散を推定して、使われる場面に最もフィットするであろう推定量を選べば良いではないか。」

なんて思いますよね。

また、「今の時代、ビッグデータの時代なのだから一致推定量のクラスで分散の小さいものを探したり、MSEが小さくなるようなものを探したほうがよっぽど建設的なのでは？」なんて疑問も持たれるかもしれません。

私も初めて不偏推定量について学んだとき思いましたし、前者については今もそう思っています。

今あげた疑問はコンピュータの計算技術が発展した現代においては(恐らく)真っ当な批判で、当然今はそういうシミュレーションが個人のPCレベルでも簡単に出来てしまうのだから、そうした環境にいる私たちは、小難しい計算を前に(面倒臭がって)そんな文句を言ってしまう訳です。

十分統計量が役に立っていた頃に思いを馳せる

コンピュータシミュレーションが簡単に行える今の時代においては、例えば大学院で指導教官に「僕は十分統計量の研究がしたい」なんて言っても、「うーん。まあ止やしないし、何も出てこないとまでは言わないけど、そこまで面白い話はないと思うよ〜。」と言われてしまうような分野です。

これは、「もうすでに一昔前に多くの研究者が議論し尽くした分野であり」「今や、十分統計量が役に立つ場面は少なくなってしまった」ためです。

そもそも、何で十分統計量や不偏推定量といった分野が、こうも厳密に議論される必要があったのでしょうか。
もう皆さん御察しの通りかと思いますが、これらの概念は、コンピュータで大規模なシミュレーションや計算を行うことが難しい時代に盛んに研究されていた概念で、要は陽に計算できる範疇で、理論上わかる範囲で如何に推定量に対して数学的な正当性を与えるか、妥当性を与えるかが研究されていた時代の考え方なのです。

だから「機械学習」「計算機統計学」「ニューラルネット」などと、どんどんコンピュータのスペックありきの手法が登場する現代で統計学を学ぼうとすると、どうにも十分統計量や不偏推定量のイメージが湧きにくい。この分野、なんか役に立つのか？？なんて思ってしまうわけですね。

当然、役立つ分野は今でもあります。特に私の所属する経済学系の研究では、不偏性を持つことは非常に重要な意味を持ちますし、理論的な正当性を完備十分性の観点から考えることも理論研究では必要かもしれません。

ただ、そうした分野に属していないとどうしても「何のための理論なのかよくわからん?」となってしまいます。

統計検定、又は数理統計学の勉強をし始めているという方が本を開いてまず、詰まるのはここだと思いますので、とりあえずは研究の道に進む方以外は、「シミュレーションや大サンプルで正当化するのが難しかった時代は、こんなに細々と理論を積み上げていったのだなあ」と昔話でも読むように、大鏡や源氏物語を読んで、昔の日本人の生活に思いを馳せるかのように理解すれば十分かなと思います。

当然、俺は研究の道に進むので、あるいは教員の道に進むので、しっかり理解しなければならないんだという方であれば、過去の論文を遡ったり、昔から読まれてきている名著を読んで、じっくり格闘する必要はあるかと思います。

十分統計量についてしっかり理解するには、必然測度論が必要になってくるので、まずはそこからでしょうか。

でも、当然勉強し始めの段階では、そういった事情、背景を知らない方がほとんどのはずですので、

「どうか、ここの分野で挫折しないよう、時間との兼ね合いですが、あまり深入りしすぎないようにした方がいいですよ」

と勉強中の皆さんにはお伝えしておきたいです(昔、何も知らずに深入りしすぎて、ひどく時間を溶かした過去が....)

2019-05-01

統計検定1級対策問題集~最尤推定量編~

統計検定対策問題統計検定統計学基礎統計

統計検定1級では、最尤推定量を求める問題が頻出しています。そこで、最尤推定量を求める問題についてまとめました。

一様分布
- パラメータ1つの場合
- パラメータが2つの場合
ベルヌーイ分布
ポアソン分布
正規分布
指数分布

一様分布

パラメータ1つの場合

最大値が未知パラメータの次のような一様分布を考えます。

$f(x)=\frac{1}{\theta}$

$0＜x≦\theta$

この時の、パラメータ $\theta$ の最尤推定量を求めます。

n個のサンプルが得られたとして、その同時尤度 $l(\theta)$ は

$l(\theta)=\frac{1}{\theta^n}$

この尤度 $l(\theta)$ を最大化する推定量 $\hat{\theta}$ を考えます。
ただし、暗黙の条件に次の条件があることに注意します。

$max\{x_1,\cdots,x_n\}≦\theta$

不等式制約があるので、正の変数 $\lambda$ を用いて、次のようなKKT条件を考える。

$L(\theta,\lambda)=\frac{1}{\theta^n}-\lambda (\theta-x_{\{max\}})$

$\frac{\partial L(\theta,\lambda)}{\partial \theta} =-\frac{n}{\theta^{n+1}}-\lambda=0$ …①

$\frac{\partial L(\theta,\lambda)}{\partial \lambda} =-(\theta-x_{\{max\}})≦0$ …②

$\lambda (\theta-x_{\{max\}})=0$ …③

$\lambda=0$ であるとすると、①の等式が成り立たない。

よって $\hat{\theta}=x_{\{max\}}$

ちなみに、この問題は実際の試験で出題されています。
また、最大値ではなく最小値が未知パラメータである場合も同様の方法で $\hat{\theta}=x_{\{min\}}$ と得られます。

パラメータが2つの場合

最大値も最小値も未知パラメータであるような一様分布として、次のようなものを考えます。

$f(x)=\frac{1}{\theta_2-\theta_1}$

$\theta_1＜x＜\theta_2$

この時、n個サンプルを得た時の同時尤度 $l(\theta_1,\theta_2)$ は次のようになります。

$l(\theta_1,\theta_2)=\frac{1}{(\theta_2-\theta_1)^n}$

更に、暗黙的に次の条件があることがわかります。

$\theta_1≦x_{\{min\}}$

$x_{\{max\}}≦\theta_2$

不等式制約なので、正の変数 $\lambda_1,\lambda_2$ を用いて次のようなKKT条件を考えます。

$L(\theta_1,\theta_2,\lambda_1,\lambda_2)=\frac{1}{(\theta_2-\theta_1)^n}-\lambda_1(x_{\{min\}}-\theta_1)-\lambda_2(\theta_2-x_{\{max\}})$

$\frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \theta_1}=\frac{n}{(\theta_2-\theta_1)^{n+1}} +\lambda_1=0$ …①

$\frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \theta_2}=\frac{-n}{(\theta_2-\theta_1)^{n+1}} -\lambda_2=0$ …②

$\frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \lambda_1}=(x_{\{min\}}-\theta_1)≦0$ …③

$\frac{\partial L(\theta_1,\theta_2,\lambda_1,\lambda_2)}{\partial \lambda_2}=(\theta_2-x_{\{max\}})≦0$ …④

$\lambda_1(x_{\{min\}}-\theta_1)=0$ …⑤

$\lambda_2(\theta_2-x_{\{max\}})=0$ …⑥

⑤、⑥について、もし $\lambda_1=0,\lambda_2=0$ であったとすると、①、②の等式は成り立たない。

よって、それぞれのパラメータの最尤推定量は

$\hat{\theta_1}=x_{\{min\}}$

$\hat{\theta_2}=x_{\{max\}}$

ベルヌーイ分布

ベルヌーイ分布の確率質量関数 $f(x)$ は次のよう。

$f(x)=p^x(1-p)^{1-x}$

ここで、n個のサンプルを考えた場合の同時尤度 $l(p)$ は

$l(p)=p^{\{\sum_{i=1}^n x_i\}} (1-p)^{\{n-\sum_{i=1}^n x_i\}}$

ここで、 $l(p)$ の対数を取ると次のよう。

$ln(p)=(\sum_{i=1}^n x_i) log(p) +(n-\sum_{i=1}^n x_i) log(1-p)$

ここで $p$ の最尤推定量 $\hat{p}$ は次の等式を満たす $p$ に等しい。

$\frac{\partial ln(p)}{\partial p}=\frac{\sum_{i=1}^n x_i}{p}-\frac{n-\sum_{i=1}^n x_i}{1-p}=0$

よって
$\hat{p}=\frac{\sum_{i=1}^n x_i}{n}$

ポアソン分布

ポアソン分布の確率質量関数 $f(x)$ は次の通り。

$f(x)=\frac{e^{-\lambda} \lambda^{x}}{x!}$

サンプルをn個取り出したとすると同時尤度 $l(\lambda)$ は次のよう。

$l(\lambda)=\frac{e^{-n\lambda} \lambda^{\sum_{i=1}^n x_i}}{\prod_{i=1}^n (x_i !)}$

ここで、対数を取ってやると

$ln(\lambda)=-log(\prod_{i=1}^n (x_i !))-n\lambda +(\sum_{i=1}^n x_i)log(\lambda)$

よって $\lambda$ の最尤推定量 $\hat{\lambda}$ は次の等式を満たす $\lambda$ に等しい。

$\frac{\partial ln(\lambda)}{\partial \lambda}=-n+\frac{\sum_{i=1}^n x_i}{\lambda}=0$

以上より

$\hat{\lambda}=\frac{\sum_{i=1}^n x_i}{n}$

正規分布

正規分布の確率密度関数は次の通り。

$f(x)=\frac{1}{\sqrt{2\pi \sigma^2}} exp(-\frac{((x-\mu)^2)}{2\sigma^2})$

同時尤度は
$l(\mu,\sigma)=(2 \pi \sigma^2)^{\frac{-n}{2}} exp(-\frac{\sum_{i=1}^n (x_i-\mu)^2}{2 \sigma^2})$

ここで対数を取ると

$ln(\mu,\sigma)=\frac{-n}{2} log(2 \pi \sigma^2)-\frac{\sum_{i=1}^n (x_i-\mu)^2}{2 \sigma^2}$

この時、 $\mu,\sigma$ の最尤推定量 $\hat{\mu},\hat{\sigma}$ は次の等式から得られる。

$\frac{\partial ln(\mu,\sigma)}{\partial \mu}=\frac{\sum_{i=1}^n (x_i-\mu)}{\sigma^2}=0$

$\frac{\partial ln(\mu,\sigma)}{\partial \sigma}=\frac{-n}{\sigma}+\frac{\sum_{i=1}^n (x_i-\mu)^2}{\sigma^3}=0$

よって、

$\hat{\mu}=\frac{\sum_{i=1}^n x_i}{n}$

$\hat{\sigma}=\sqrt{\frac{\sum_{i=1}^n (x_i-\hat{\mu})^2}{n}}$

指数分布

指数分布の確率密度関数は次の通り

$f(x)=\lambda exp(-\lambda x)$

この時同時尤度は

$l(\lambda)=\lambda^n exp(-\lambda \sum_{i=1}^n x_i)$

ここで、対数を取ると

$ln(\lambda)=n log(\lambda) -\lambda \sum_{i=1}^n x_i$

この時、 $\lambda$ の最尤推定量 $\hat{\lambda}$ は次の等式から得られる。

$\frac{\partial ln(\lambda)}{\partial \lambda}=\frac{n}{\lambda}-\sum_{i=1}^n x_i=0$

よって、

$\hat{\lambda}=\frac{n}{\sum_{i=1}^n x_i}$

2019-04-30

統計検定1級対策問題集~十分統計量編2~

統計検定対策問題統計検定統計学基礎統計負の二項分布

十分統計量に関する問題2記事目です。

ラオブラックウェルの定理
- ラオブラックウェルの定理とは
- ラオブラックウェルの定理証明
フィッシャーネイマンの分解定理
リンク

ラオブラックウェルの定理

ラオブラックウェルの定理とは

初めに完備十分統計量を考える上で重要になってくる「ラオブラックウェルの定理」の証明についてやっておきます。

実際の数理統計の本であれば、ラオブラックウェルの定理を示した後、完備性について解説し、完備十分統計量について話が移ります。

ただ、統計検定1級公式参考書では完備十分統計量までは触れていません。
ラオブラックウェルの定理までの説明で終わっているので、とりあえずこの定理の証明までは確認しておきます。

次のような定理をラオブラックウェルの定理と呼びます。
$T$ を $\theta$ の十分統計量とする。ここで、 $\theta$ のある推定量 $\delta(X)$ について、次のような推定量 $\delta_1(T)$ をラオブラックウェル推定量と呼ぶことにする。

$\delta_1(T)=E_{\theta}[\delta(X)|T]$

そして、ラオブラックウェル推定量が満たす次のような性質をラオブラックウェルの定理と呼ぶ。

$E_{\theta}[(\delta_1(T)-\theta)^2]≦E_{\theta}[(\delta(X)-\theta)^2]$

不等式の両サイドは平均二乗誤差になっています。
つまり、この不等式からわかることは

「ある推定量 $\delta$ を考えた時に、それよりも平均二乗誤差を小さくする(又は同等)推定量を、十分統計量を条件付けることで考えることが出来る」

ということです。単純ですが強力な定理です。

ラオブラックウェルの定理証明

まず、 $E_{\theta}[\delta_1(T)]=E_{\theta}[\delta(X)]$ であることを示します。

$E_{\theta}[\delta_1(T)]=\int_T \int_X \delta(X) dP(X|T) dP(T)$

$=\int_X \delta(X) dP(X)=E_{\theta}[\delta(X)]$

また、 $E_{\theta}[\delta_1(T)^2]≦E_{\theta}[\delta(X)^2]$ であることも示します。

$E_{\theta}[\delta_1(T)^2]=E_{T;\theta}[ E_{\theta}[\delta(X)|T]^2]$

$E_{\theta}[\delta(X)^2]=E_{T;\theta}[E_{\theta}[\delta(X)^2|T]]$

更にイェンゼン不等式を用いて、
$E_{T;\theta}[ E_{\theta}[\delta(X)|T]^2]≦E_{T;\theta}[E_{\theta}[\delta(X)^2|T]]$

以上より示せた。

最後にラオブラックウェルの定理を示します。
$E_{\theta}[(\delta_1(T)-\theta)^2]=E_{\theta}[\delta_1(T)^2]-2\theta E_{\theta}[\delta_1(T)]+\theta^2$

$E_{\theta}[(\delta(X)-\theta)^2]=E_{\theta}[\delta(X)^2]-2\theta E_{\theta}[\delta(X)]+\theta^2$

上二つの性質から
$E_{\theta}[(\delta_1(T)-\theta)^2]≦E_{\theta}[(\delta(X)-\theta)^2]$

フィッシャーネイマンの分解定理

負の二項分布

負の二項分布の確率質量関数 $f(x;p,r)$ は次のようになります。

$f(x;p,r)=\begin{eqnarray*} && {}_{r+x-1} C _x \\ \end{eqnarray*} p^r (1-p)^x$

$x_1,\cdots,x_n$ のサンプルを独立同一に得たとすると、同時分布は

$P(x_1,\cdots,x_n ;n,p,r) = \prod_{i=1}^n \{ \begin{eqnarray*} && {}_{r+x_i-1} C _{x_i} \\ \end{eqnarray*} p^r (1-p)^{x_i} \}$

$= \{ \prod_{i=1}^n \begin{eqnarray*} && {}_{r+x_i-1} C _{x_i} \\ \end{eqnarray*} \} p^{nr} (1-p)^{ \sum_{i=1}^n x_i }$

この時、 $T(X)=\sum_{i=1}^n x_i$ がパラメータ $p$ の十分統計量であることを示します。

フィッシャーネイマンの分解定理より、

$h(X)=\{\prod_{i=1}^n \begin{eqnarray*} && {}_{r+x_i-1} C _{x_i} \\ \end{eqnarray*} \}$

$g(T(X),p)=p^{nr}(1-p)^{T(X)}$

とみると、 $T(X)=\sum_{i=1}^n x_i$ がパラメータ $p$ の十分統計量であることがわかる。

ガンマ分布

ガンマ分布の確率密度関数 $f(x;\alpha,\beta)$ は次のようになります。

$f(x;\alpha,\beta)=\frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}$

$x_1,\cdots,x_n$ のサンプルを独立同一に得たとすると、

$T_{\alpha}(X)=\prod_{i=1}^n x_i$ は $\alpha$ の十分統計量
$T_{\beta}(X)=\sum_{i=1}^n x_i$ は $\beta$ の十分統計量です。これを示します。

同時分布は

$P(x_1,\cdots,x_n ;\alpha,\beta)=\frac{\beta^{n\alpha}}{\Gamma(\alpha)^n} \{\prod_{i=1}^n x_i\}^{\alpha-1} e^{-\beta\sum_{i=1}^n x_i}$

よって、

$h_{\alpha}(X)=e^{-\beta\sum_{i=1}^n x_i}$

$g(T_{\alpha}(X),\alpha)=\frac{\beta^{n\alpha}}{\Gamma(\alpha)^n} \{T_{\alpha}(X)\}^{\alpha-1}$

とみると、フィッシャーネイマンの分解定理より、 $T_{\alpha}(X)=\prod_{i=1}^n x_i$ は $\alpha$ の十分統計量。

また、 $h_{\beta}(X)= \frac{1}{\Gamma(\alpha)^n} \{\prod_{i=1}^n x_i\}^{\alpha-1}$

$g(T_{\beta}(X),\beta)=\beta^{n\alpha} e^{-\beta T_{\beta}(X)}$

と見ると、フィッシャーネイマンの分解定理より、 $T_{\beta}(X)=\sum_{i=1}^n x_i$ は $\beta$ の十分統計量

一様分布

下限が0,上限が未知パラメータ $\theta$ であるときの一様分布を考えます。この時、密度関数は次のようになります。

$f(x;\theta)=\frac{1}{\theta}$ $(0≦x≦\theta)$

これは、見方を変えると次のように書くことも出来ます。

$f(x;\theta)=\frac{1_{\{0≦x≦\theta\}}}{\theta}$

ここで、 $x_1,\cdots,x_n$ のサンプルを独立同一に得たとすると、 $\theta$ の十分統計量は $max\{x_1,\cdots , x_n\}$ となります。

これを示します。

同時分布は次のようになるので

$P(x_1,\cdots x_n;n,\theta)=\frac{1_{\{max\{x_1,\cdots , x_n\}≦\theta\}}}{\theta^n}$

フィッシャーネイマンの分解定理より、 $\theta$ の十分統計量は $max\{x_1,\cdots , x_n\}$ であることは明らか。

リンク

統計検定1級対策問題集~十分統計量編１~ - バナナでもわかる話