上極限集合/下極限集合とボレル・カンテリの補題

上極限集合と下極限集合

定義

\(A_{1}, A_{2}, \dots\) を集合の列とする。このとき、上極限集合下極限集合をそれぞれ

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}A_{k}$$
$$\lim_{n\to\infty}\inf A_{n}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}$$

で定義する。

数式による解釈

定義式の各パーツについて、条件式を用いて書き直すと、

$$\bigcup_{k=n}^{\infty}A_{k}=\{x\mid\exists l\geq n \quad x\in A_{l}\}$$
$$\bigcap_{k=n}^{\infty}A_{k}=\{x\mid\forall l\geq n \quad x\in A_{l}\}$$

であるので、定義式は

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}\{x\mid\exists l\geq n \quad x\in A_{l}\}$$
$$=\{x\mid\forall m\quad x\in\{x\mid\exists l\geq m \quad x\in A_{l}\}\}$$
$$=\{x\mid\forall m\quad\exists l\geq m \quad x\in A_{l}\}$$
$$\lim_{n\to\infty}\inf A_{n}=\bigcup_{n=1}^{\infty}\{x\mid\forall l\geq n \quad x\in A_{l}\}$$
$$=\{x\mid\exists m\quad x\in\{x\mid\forall l\geq m \quad x\in A_{l}\}\}$$
$$=\{x\mid\exists m\quad\forall l\geq m \quad x\in A_{l}\}$$

と表すことができる。

言語による解釈

前節の議論により、上極限集合は「どんな\(m\)を取ってきても、それ以上の\(l\)が存在し、\(x\)が\(A_{l}\)に含まれる」、すなわち「\(m\)以上のところに、必ず\(x\)を含む集合\(A_{l}\)が存在する」という\(x\)の集合である。したがって、上極限集合の要素は無限個の\(A_{k}\)に含まれる。

また、下極限集合は「とある\(m\)以上のすべての\(l\)について、\(x\)が\(A_{l}\)に含まれる」という\(x\)の集合である。これは裏返して言えば「\(x\)を含まない集合は最大でも\(A_{m}\)まで」ということであり、したがって、\(x\)を含まない\(A_{k}\)の数は有限個ということになる。

集合の列\(A_{1},A_{2},\dots\)について、上極限集合の元\(x\)について\(x\in A_{k}\)となるとき、\(A_{k}\)を\(T\)、そうでないときは\(F\)と置き換えることにする。このとき、集合の列は

$$TFFTFTT\dots FFT\dots FT\dots$$

のようになる。すなわち、集合の列において、あるところから急にTが出なくなるということはない。

下極限集合についても同様の操作を行うと、

$$TFFTFTT\dots TTT\dots TT\dots$$

のようになる。すなわち、あるところから急にTしか出なくなる。

性質

前節の例で見た通り、上極限集合の元が満たすべき条件は、下極限集合の元のそれよりも緩い。したがって、

$$\lim_{n\to\infty}\sup A_{n}\supset\lim_{n\to\infty}\inf A_{n}$$

の関係が成り立つ。上極限集合と下極限集合が一致するとき、集合列は収束すると言い

$$\lim_{n\to\infty}A_{n}$$

と書くことがある。

集合列が収束する例(単調増大列と単調減少列)

集合の列\(A_{k}\in\mathfrak{B}, k=1,2,\dots,\)において、すべての\(k\)について\(A_{k}\subset A_{k+1}\)が成り立つとき、\(A_{k}\)は単調増大列であると言い、逆にすべての\(k\)について\(A_{k}\supset A_{k+1}\)が成り立つとき、\(A_{k}\)は単調減少列であると言う。

\(A_{k}\)が単調増大列のとき、以下が成り立つ。

$$すべての m (\geq 1)について、\bigcup_{k=m}^{\infty}A_{k}は等しい\tag{1}$$
$$\bigcap_{k=m}^{\infty}A_{k}=A_{m}\tag{2}$$

これらを用いると、上極限集合と下極限集合はそれぞれ

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}A_{k}=\bigcup_{n=1}^{\infty}A_{n}\quad(\because(1))$$
$$\lim_{n\to\infty}\inf A_{n}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}=\bigcup_{n=1}^{\infty}A_{n}\quad(\because(2))$$

と表すことができ、互いに等しくなる。したがって、集合列は収束し

$$\limsup_{n\to\infty}A_{n}=\liminf_{n\to\infty}A_{n}=\lim_{n\to\infty}A_{n}=\bigcup_{n=1}^{\infty}A_{n}\tag{3}$$

が成り立つ。

\(A_{k}\)が単調減少列のとき、\(A_{k}\)の補集合の列\(A_{k}^{c}\)について\(A_{k}^{c}\supset A_{k+1}^{c}\)が成り立ち単調増大列となるので、上の結果から

$$\limsup_{n\to\infty}A_{n}^{c}=\liminf_{n\to\infty}A_{n}^{c}=\lim_{n\to\infty}A_{n}^{c}=\bigcup_{n=1}^{\infty}A_{n}^{c}$$

補集合の列が収束するので、単調減少列も収束し

$$\limsup_{n\to\infty}A_{n}=\liminf_{n\to\infty}A_{n}= \lim_{n\to\infty}A_{n}$$
$$=\left(\lim_{n\to\infty}A_{n}^{c}\right)^{c}=\left(\bigcup_{n=1}^{\infty}A_{n}^{c}\right)^{c}=\bigcap_{n=1}^{\infty}A_{n}\tag{4}$$

まとめ

上極限集合:「無限個の\(A_{k}\)に含まれる元の集合」
下極限集合:「それを含まない\(A_{k}\)が有限個である元の集合」

ボレル・カンテリの補題

内容

\(\mathfrak{B}\)を可測集合族(確率を定義した事象の集合。元\(A\)は可測集合と呼ばれる)とする。

  1. 事象の列\(A_{n}\in\mathfrak{B}, n=1,2,\dots,\)に対して、\(\sum_{n=1}^{\infty}P(A_{n})<\infty\)ならば、\(P(\limsup_{n\to\infty}A_{n})=0\)となる。
  2. 事象\(A_{1},A_{2},\dots,A_{n},\dots\in\mathfrak{B}\)が独立であるとき、\(\sum_{n=1}^{\infty}P(A_{n})=\infty\)ならば、\(P(\limsup_{n\to\infty}A_{n})=1\)が成り立つ。

意義

\(A_{n}\)は「サイコロを振って、\(n\)回目に\(1\)が出る」という事象であるとし、事象列\(A_{1},A_{2},\dots,\)を考える。前章の解釈より、この事象列の上極限は「サイコロを無限回(\(n\to\infty\))振ると、\(1\)が無限回出る」という事象、下極限は「サイコロを無限回振ると、\(1\)が出ない(\(2,3,4,5,6\)が出る)回数は有限回である」という事象をそれぞれ表す。すると直感的に

$$P(\lim_{n\to\infty}\sup A_{n})=1$$
$$P(\lim_{n\to\infty}\inf A_{n})=0$$

となることが予想されるが、これが実際に満たされるための十分条件を与えるのが、ボレル・カンテリの補題である。

1.の証明

\(B_{n}=\bigcup_{k=n}^{\infty}A_{k}\)とおくと、\(B_{n}\)は単調減少列となる。上極限集合の定義式と\((4)\)より

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}B_{n}=\lim_{n\to\infty}B_{n}$$

確率の連続性により

$$P(\lim_{n\to\infty}\sup A_{n})=P(\lim_{n\to\infty}B_{n})=\lim_{n\to\infty}P(B_{n})=\lim_{n\to\infty}P\left(\bigcup_{k=n}^{\infty}A_{k}\right)$$

ここで、\(P(\bigcup_{k=n}^{\infty})\leq\sum_{k=n}^{\infty}P(A_{k})\)であり、\(\sum_{n=1}^{\infty}P(A_{n})<\infty\)のとき\(\lim_{n\to\infty}\sum_{k=n}^{\infty}=0\)となることに注意して変形すると

$$P(\lim_{n\to\infty}\sup A_{n})=\lim_{n\to\infty}P\left(\bigcup_{k=n}^{\infty}A_{k}\right)\leq\lim_{n\to\infty}\sum_{k=n}^{\infty}P(A_{k})=0$$

が成り立つ。

2.の証明

$$(\lim_{n\to\infty}\sup A_{n})^{c}=\lim_{n\to\infty}\inf A_{n}^{c}$$

すなわち

$$\left(\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}A_{k}\right)^{c}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c}$$

であることに注意すると、\(P(\limsup_{n\to\infty}A_{n})=1\)は\(P(\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c})=0\)に等しいため、こちらを証明する。まず、

$$P\left(\bigcap_{k=n}^{\infty}A_{k}^{c}\right)=\lim_{m\to\infty}P\left(\bigcap_{k=n}^{m}A_{k}^{c}\right)=\lim_{m\to\infty}\prod_{k=n}^{m}P(A_{k}^{c})=\lim_{m\to\infty}\prod_{k=n}^{m}\left(1-P(A_{k})\right)$$

が成り立つ。ここで、\(x\geq0\)において常に\(1-x\leq e^{-x}\)より

$$\lim_{m\to\infty}\prod_{k=n}^{m}\left(1-P(A_{k})\right)\leq \lim_{m\to\infty}\prod_{k=n}^{m}\exp\left(-P(A_{k})\right)$$
$$=\lim_{m\to\infty}\exp\left(-\sum_{k=n}^{m}P(A_{k})\right)=0\quad(\because\quad 仮定)$$

がすべての\(n\geq1\)について成り立ち、すなわち\(P(\bigcap_{k=n}^{\infty}A_{k}^{c})=0\)である。ここで、\(B_{n}=\bigcap_{k=n}^{\infty}A_{k}^{c}\)とおくと\(B_{n}\)は単調増大列であり、\((3)\)から

$$\lim_{n\to\infty}B_{n}=\bigcup_{n=1}^{\infty}B_{n}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c}$$

となる。したがって、

$$P\left(\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c}\right)=P\left(\lim_{n\to\infty}B_{n}\right)=\lim_{n\to\infty}P\left(B_{n}\right)$$
$$=\lim_{n\to\infty}P\left(\bigcap_{k=n}^{\infty}A_{k}^{c}\right)=0$$

「意義」の節で述べたようなサイコロ投げの場合、事象\(A_{1},A_{2},\dots,A_{n},\dots\in\mathfrak{B}\)は独立であり、\(\sum_{n=1}^{\infty}P(A_{n})=\sum_{n=1}^{\infty}\frac{1}{6}=\infty\)より、直感通り\(P(\limsup_{n\to\infty}A_{n})=1\)となる。

しかし、\(P(A_{n})=\frac{1}{n^{2}}\)となるような事象を考えた場合、

$$\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\frac{\pi^{2}}{6}<\infty\quad(\because\quad バーゼル問題)$$

より、\(P(\limsup_{n\to\infty}A_{n})=0\)となる。 すなわち、このような事象は無限回は起こり得ない。

参考文献

  1. 久保川達也「現代数理統計学の基礎(共立講座 数学の魅力11)」共立出版(2017)
  2. 佐藤坦「はじめての確率論 測度から確率へ」共立出版(1994)
  3. 極限集合の解釈(イプシロンデルタ風に) – 岡竜之介のブログ(http://agajo.hatenablog.com/entry/2016/10/15/111717

複素数の指数関数・対数関数・べき乗と、その微分公式

以下、特に断りなく、複素数 \(z\) について

$$z=x+iy\quad(x, y は実数)$$

とおき、\(z\) の実部を \(x\) 、虚部を \(y\) で表す。

指数関数

定義

$$f(z)=e^z=e^{x+iy}=e^xe^{iy}$$
$$=e^x(\cos{y}+i\sin{y})\tag{1}$$

性質

定義式の形状より、指数関数は \(y\) 方向(虚軸方向)について周期関数となっている。すなわち

$$e^z=e^{z+2\pi i}$$

微分公式

$$f^{\prime}(z)=f(z)=e^z\tag{2}$$

微分は複素平面全体において定義され、指数関数は正則である。

対数関数

定義

実関数において、対数関数は指数関数の1価の(1対1対応する)逆関数として定義できたが、上記の通り、複素数における指数関数は周期関数となるため、対数関数の複素数への拡張は多価関数(1対多対応)となる。
定義式は、\(z\) が極形式で \(z=re^{i\theta}\) と表せることを用いて

$$\log{z}=\log{r}+i\theta\quad(z\neq 0)\tag{3}$$

となる。ここで \(r, \theta\) はそれぞれ、\(z\) の原点からの距離と偏角を表している。

定義式において、左辺の \(log\) は複素数に拡張した対数関数であるが、右辺の \(log\) は正の実数に対する対数関数であるため、意味が異なるので注意する。

また、偏角 \(\theta\) には \(2\pi\) の任意性があり

$$\theta=\arg{z}={\rm Arg}\ z+2n\pi\quad(n:整数)$$

と表せる。ここでは \(-\pi\leq {\rm Arg}\ z<\pi\) とし、整数 \(n\) は任意の値を取れることから、\(\log{z}\) は無限多価関数となる。

上式において \(n=0\) としたものは対数関数の主値と呼ばれ、

$${\rm Log}\ z=\log{r}+i{\rm Arg}\ z\tag{3*}$$

と表記する。

性質

偏角の主値の定義を \(-\pi\leq {\rm Arg}\ z<\pi\) としているとき、\(\log{z}\) は負の実軸と原点で連続ではない。

証明

実部が正で、虚部が十分に小さい複素数 \(z_{1}, z_{2}\) は、十分に小さい \(\delta >0, \epsilon >0\) を用いて

$$z_{1}=x+i\delta=xe^{i\epsilon}$$
$$z_{2}=x-i\delta=xe^{-i\epsilon}$$

と書ける( \(x>0\) )。
ここで、\(e^{i\theta}=\cos{\theta}+i\sin{\theta}\) であり、\(\theta\) が十分小さい時、\(\cos{\theta}\simeq 1,\ \sin{\theta}\simeq\theta\) であることから

$$xe^{i\epsilon}\simeq x(1+i\epsilon)=x+ix\epsilon$$

であり、\(x\epsilon=\delta\) と置くことで上式を得た。

よって \(z_{1}, z_{2}\) について対数関数の主値を取り、 \(\epsilon\to +0(\delta\to + 0)\)とすると

$${\rm Log}\ z_{1}=\log{x}+i\epsilon\to\log{x}$$
$${\rm Log}\ z_{2}=\log{x}-i\epsilon\to\log{x}$$

であり、正の実軸では連続である。

同様に、実部が負の複素数 \(z_{3}, z_{4}\) は

$$z_{3}=-x+i\delta=xe^{i(\pi-\epsilon)}$$
$$z_{4}=-x-i\delta=xe^{i(-\pi+\epsilon)}$$

と書ける( \(x>0\) )。対数関数の主値を取り、 \(\epsilon\to +0\)とすると

$${\rm Log}\ z_{3}=\log{x}+i(\pi-\epsilon)\to\log{x}+\pi i$$
$${\rm Log}\ z_{4}=\log{x}-i(\pi-\epsilon)\to\log{x}-\pi i$$

となり、負の実軸に近づくとき、極限値が異なる。
よって、対数関数は負の実軸において不連続である。

上記の性質は偏角の主値 \({\rm Arg}\ z\) の取り方に依存するため、注意が必要である。

その他の一般的な性質

複素数の対数関数は多価関数であるため、実数の対数関数がにおいて見られた性質が必ずしも成立しないことに注意する。

  • \(e^{\log{z}}=z\)
  • \(\log{e^z}\neq z\)
  • \(\log{z_{1}z_{2}}=\log{z_{1}}+\log{z_{2}}\)
  • \({\rm Log}\ {z_{1}z_{2}}\neq{\rm Log}\ {z_{1}}+{\rm Log}\ {z_{2}}\)
  • \(\log{z}\neq n\log{z}\)

微分公式

$$({\rm Log}\ z)^{\prime}=\frac{1}{z}\tag{4}$$
$$(\log{z})^{\prime}=\frac{1}{z}\tag{5}$$

対数関数の主値 \({\rm Log}\ z\) は、領域 \(|z|>0,\ -\pi\leq{\rm Arg}\ z<\pi\) において1価かつ正則である。

証明

\({\rm Log}\ z\) が1価であることは明らか。

領域 \(|z|>0,\ -\pi\leq{\rm Arg}\ z<\pi\) において \((3*)\) 式より、\(u(r,\theta)=\log{r},\ v(r,\theta)={\rm Arg}\ z=\theta\) とする。ここで、\(z\) の偏角 \(\theta\) は主値に固定した。
\(u,v\) を \(r,\theta\) で偏微分すると

$$\frac{\partial u}{\partial r}=\frac{1}{r},\ \frac{\partial u}{\partial \theta}=0$$
$$\frac{\partial v}{\partial r}=0,\ \frac{\partial v}{\partial \theta}=1$$

より、コーシー=リーマン方程式

$$\frac{\partial u}{\partial r}=\frac{1}{r}\frac{\partial v}{\partial \theta}$$
$$\frac{\partial v}{\partial r}=-\frac{1}{r}\frac{\partial u}{\partial \theta}$$

が成立するので、\({\rm Log}\ z\) は領域 \(|z|>0,\ -\pi\leq{\rm Arg}\ z<\pi\) において微分可能、すなわち正則である。

よって、複素関数の極限が複素平面上の経路によらないことに注意して、偏角を \(\theta_{0}\) に固定し、距離 \(r\) を変化させることで \(z_{0}=r_{0}e^{i\theta_{0}}\) における微分を考える。
\(z=re^{i\theta_{0}}\) とおくと \(\Delta z=z-z_{0}=(r-r_{0})e^{i\theta_{0}}=e^{i\theta_{0}}\Delta r\) であるので

$$({\rm Log}\ z)^{\prime}=\lim_{\Delta r \to 0}\frac{1}{\Delta re^{i\theta_{0}}}[u(r_{0}+\Delta r,\theta_{0})+iv(r_{0}+\Delta r,\theta_{0}) -u(r_{0},\theta_{0})-iv(r_{0},\theta_{0})]$$
$$=(\frac{\partial u(r_{0},\theta_{0})}{\partial r}+i\frac{\partial v(r_{0},\theta_{0})}{\partial r})e^{-i\theta_{0}}$$
$$=(\frac{1}{r}+0)e^{-i\theta_{0}}=\frac{1}{z}$$

が成り立つ。

また、\(z=\exp{(\log{z})}\) より、\(u=\log{z}\) とおいて両辺を \(z\) で微分すると

$$1=\frac{d}{dz}\exp{(u)}$$
$$1=\frac{d}{du}\exp{(u)}\frac{du}{dz}$$
$$1=\exp{(u)}\cdot\frac{du}{dz}$$

すなわち

$$\frac{d}{dz}\log{z}=\frac{1}{\exp{(u)}}=\frac{1}{z}$$

べき乗

定義

$$z^c=\exp{(c\log{z})}\quad(c:複素数)$$

微分公式

$$(z^c)^{\prime}=cz^{c-1}$$
$$(c^z)^{\prime}=c^z\log{c}$$
$$(z^z)^{\prime}=z^z(1+\log{z})$$

証明

\(u=c\log{z}\) とおくと

$$\frac{d}{dz}z^c=\frac{d}{dz}\exp{(u)}$$
$$=\frac{d}{du}\exp{(u)}\frac{du}{dz}$$
$$=\exp{(u)}\cdot\frac{c}{z}=z^c\cdot\frac{c}{z}=cz^{c-1}$$

\(v=z\log{c}\) とおくと

$$\frac{d}{dz}c^z=\frac{d}{dz}\exp{(v)}$$
$$=\frac{d}{dv}\exp{(v)}\frac{dv}{dz}$$
$$=\exp{(v)}\cdot\log{c}=c^z\log{c}$$

\(w=z\log{z}\) とおくと

$$\frac{d}{dz}z^z=\frac{d}{dz}\exp{(w)}$$
$$=\frac{d}{dw}\exp{(w)}\frac{dw}{dz}$$
$$=\exp{(w)}(z\cdot\frac{1}{z}+1\cdot\log{z})=z^z(1+\log{z})$$

ルジャンドル多項式の定義と性質(3/3)

定義

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n,\quad(n=0,1,2,…)\tag{1}$$

で定義された関数 \(P_{n}(x)\) を、ルジャンドル(Legendre)多項式という。
また、ルジャンドル多項式は、以下のテイラー級数

$$\frac{1}{\sqrt{1-2xt+t^2}}=\sum_{n=0}^{\infty}P_{n}(x)t^n\tag{2}$$

の係数として定義することもできる。

性質

ルジャンドル多項式は、以下の基本的な性質を持つ。

  1. 以下の漸化式が \(n\geq1\) に対して成り立つ。
    $$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)\tag{3}$$
  2. $$P_{n}(-x)=(-1)^n P_{n}(x)$$
  3. $$P_{n}(1)=1,\quad P_{n}(-1)=(-1)^n$$
  4. \(P_{n}(x)=0\) は相異なる \(n\) 個の実数解を開区間 \((-1, 1)\) 内に持つ。
  5. $$\int_{-1}^{1}x^k P_{n}(x)d x = 0 \quad (n\geq1,\quad k=0,1,2,…,n-1)$$
  6. 高々 \(n-1\) 次の多項式 \(Q(x)\) に対して、
    $$\int_{-1}^{1}Q(x)P_{n}(x)d x = 0$$

証明

性質 1., 2., 3.

性質 4.

性質 5.

性質 4. における議論より、\(0\leq k\leq n-1\) において常に

$$f^{(k)}(\pm 1)=0$$

が成り立つ。
\(\int_{-1}^{1}x^n P_{n}(x)d x\) について、上記を踏まえて部分積分を繰り返すと、

$$\int_{-1}^{1}x^k \frac{d^n}{d x^n}(x^2-1)^n d x$$
$$=\left[x^k\frac{d^{n-1}}{d x^{n-1}}(x^2-1)^n\right]_{-1}^{1}-\int_{-1}^{1}k x^{k-1}\frac{d^n}{d x^n}(x^2-1)^n d x$$
$$=-k\int_{-1}^{1}x^{k-1}\frac{d^n}{d x^n}(x^2-1)^n d x$$
$$=k(k-1)\int_{-1}^{1}x^{k-2}\frac{d^n}{d x^n}(x^2-1)^n d x$$
$$…$$
$$=(-1)^k k(k-1)…2\cdot 1\int_{-1}^{1}\frac{d^{n-k}}{d x^{n-k}}(x^2-1)^n d x$$
$$=(-1)^k k(k-1)…2\cdot 1\left[\frac{d^{n-k-1}}{d x^{n-k-1}}(x^2-1)^n\right]_{-1}^{1}$$
$$=0\quad \because 0\leq n-k-1 \leq n-1$$

より、

$$\int_{-1}^{1}x^k P_{n}(x)d x=\frac{1}{2^n n!}\int_{-1}^{1}x^k \frac{d^n}{d x^n}(x^2-1)^n d x=0$$

が示された。

性質 6.

高々 \(n-1\) 次の多項式 \(Q(x)=\sum_{k=0}^{n-1}a_{k}x^{k}\) に対して、

$$\int_{-1}^{1}Q(x)P_{n}(x)d x=\int_{-1}^{1}(\sum_{k=0}^{n-1}a_{k}x^{k})P_{n}(x)d x$$
$$=\sum_{k=0}^{n-1}a_{k}\int_{-1}^{1}x^{k}P_{n}(x)d x=0$$

より示された。

ルジャンドル多項式の定義と性質(2/3)

定義

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n,\quad(n=0,1,2,…)\tag{1}$$

で定義された関数 \(P_{n}(x)\) を、ルジャンドル(Legendre)多項式という。
また、ルジャンドル多項式は、以下のテイラー級数

$$\frac{1}{\sqrt{1-2xt+t^2}}=\sum_{n=0}^{\infty}P_{n}(x)t^n\tag{2}$$

の係数として定義することもできる。

性質

ルジャンドル多項式は、以下の基本的な性質を持つ。

  1. 以下の漸化式が \(n\geq1\) に対して成り立つ。
    $$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)\tag{3}$$
  2. $$P_{n}(-x)=(-1)^n P_{n}(x)$$
  3. $$P_{n}(1)=1,\quad P_{n}(-1)=(-1)^n$$
  4. \(P_{n}(x)=0\) は相異なる \(n\) 個の実数解を開区間 \((-1, 1)\) 内に持つ。
  5. $$\int_{-1}^{1}x^k P_{n}(x)d x = 0 \quad (n\geq1,\quad k=0,1,2,…,n-1)$$
  6. 高々 \(n-1\) 次の多項式 \(Q(x)\) に対して、
    $$\int_{-1}^{1}Q(x)P_{n}(x)d x = 0$$

証明

性質 1., 2., 3.

性質 4.

ロル(Rolle)の定理を用いる。

有界閉区間 \([a, b]\) 上で定義された連続関数 \(f(x)\) が開区間 \((a, b)\) で微分可能であり、

$$f(a)=f(b)$$

を満たすとき、

$$f^{\prime}(c)=0$$

となる \(c\in(a,b)\) が存在する。

Rolle の定理

まず、

$$f(x)=(x^2-1)^{n}=(x-1)^{n}(x+1)^{n}$$

とおき、これを \(g(x)=(x-1)^{n}\) と \(h(x)=(x+1)^{n}\) の合成関数と見ると、\(g(\pm 1)=h(\pm 1)=0\) であるから、ライプニッツ(Leibniz)の公式

$$(g(x)h(x))^{(n)}=\sum_{k=0}^{n}{}_n \mathrm{ C }_k g^{(n)}(x) h^{(n-k)}(x)$$

を順次導関数に適用することによって、\(1\leq k\leq n-1\) において常に

$$f^{(k)}(\pm 1)=0$$

が成り立つ。
また、\(f(x)\) は偶関数であり、開区間 \((-1, 1)\) で \(0\) にならない。

したがって、Rolleの定理より

$$f^{(1)}(a_{1}^{(1)})=0$$

を満たす \(a_{1}^{(1)},\quad(-1<a_{1}^{(1)}<1)\) が存在する。

したがって、\(f^{(1)}(-1)=f^{(1)}(a_{1}^{(1)})=f^{(1)}(1)=0\) であるから、\((-1, a_{1}^{(1)}), (a_{1}^{(1)}, 1)\) にRolleの定理をそれぞれ適用して、

$$f^{(2)}(a_{1}^{(2)})=f^{(2)}(a_{2}^{(2)})=0$$

を満たす \(a_{1}^{(2)}, a_{2}^{(2)}\quad(-1<a_{1}^{(2)}<(a_{1}^{(1)})<a_{2}^{(2)}<1)\) が存在することがわかる。

以上の手順を繰り返すと、定義式\((1)\)においては

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n=0$$
$$\frac{d^n}{d x^n} (x^2-1)^n=0$$

を満たす相異なる \(n\) 個の解 \(a_{1}^{(n)}, a_{2}^{(n)}, …, a_{n}^{(n)}\) が開区間 (\(-1, 1)\) に存在することが示される。

性質 5., 6.

ルジャンドル多項式の定義と性質(1/3)

定義

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n,\quad(n=0,1,2,…)\tag{1}$$

で定義された関数 \(P_{n}(x)\) を、ルジャンドル(Legendre)多項式という。
また、ルジャンドル多項式は、以下のテイラー級数

$$\frac{1}{\sqrt{1-2xt+t^2}}=\sum_{n=0}^{\infty}P_{n}(x)t^n\tag{2}$$

の係数として定義することもできる。

性質

ルジャンドル多項式は、以下の基本的な性質を持つ。

  1. 以下の漸化式が \(n\geq1\) に対して成り立つ。
    $$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)\tag{3}$$
  2. $$P_{n}(-x)=(-1)^n P_{n}(x)$$
  3. $$P_{n}(1)=1,\quad P_{n}(-1)=(-1)^n$$
  4. \(P_{n}(x)=0\) は相異なる \(n\) 個の実数解を開区間 \((-1, 1)\) 内に持つ。
  5. $$\int_{-1}^{1}x^k P_{n}(x)d x = 0 \quad (n\geq1,\quad k=0,1,2,…,n-1)$$
  6. 高々 \(n-1\) 次の多項式 \(Q(x)\) に対して、
    $$\int_{-1}^{1}Q(x)P_{n}(x)d x = 0$$

証明

性質 1.

定義式 \((2)\) の両辺を \(t\) で微分して、

$$\frac{d}{d t}\frac{1}{\sqrt{1-2xt+t^2}}=\frac{d}{dt}\sum_{n=0}^{\infty}P_{n}(x)t^n$$
$$-\frac{1}{2}(-2x+2t)(1-2xt+t^2)^{-\frac{3}{2}}=\sum_{n=1}^{\infty}n P_{n}(x)t^{n-1}$$
$$(x-t)\frac{1}{\sqrt{1-2xt+t^2}}= (1-2xt+t^2) \sum_{n=0}^{\infty}(n+1)P_{n+1}(x)t^{n}$$

ここで、根号を定義式\((2)\)により置き換えると、

$$(x-t) \sum_{n=0}^{\infty}P_{n}(x)t^n = (1-2xt+t^2) \sum_{n=0}^{\infty}(n+1)P_{n+1}(x)t^{n}$$

が得られる。
この式の両辺において、\(t^n\) の係数について比較すると、

左辺では

$$x P_{n}(x)-P_{n-1}(x)$$

右辺では

$$(n+1)P_{n+1}(x)-2xn P_{n}(x)+(n-1)P_{n-1}(x)$$

より、

$$x P_{n}(x)-P_{n-1}(x)=(n+1)P_{n+1}(x)-2xn P_{n}(x)+(n-1)P_{n-1}(x)$$
$$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)$$

が得られる。
この \(P_{n}(x)\) に関する漸化式は、ボネの漸化式と呼ばれる。

性質 2.

定義式 \((1)\) において、\((x^2-1)^n\) を二項定理を用いて展開すると \(x\) の \(2n\) 次の多項式になり、かつ、\(x\) の偶数乗の項のみから成り立っている。
これを \(n\) 回微分すると、\(x\) の \(n\) 次の多項式となり、かつ、各項の次数は \(2\) ずつ下がっていく。
すなわち、\(P_{n}(x)\) は \(n\) が偶数の時には偶関数、\(n\) が奇数の時には奇関数となる。これは性質 2. と同義である。

性質 3.

\((3)\)式を変形し、\(n\geq 0\) において

$$(n+2)P_{n+2}(x)=(2n+3)x P_{n+1}(x)-(n+1)P_{n}(x)\tag{3*}$$

とする。

$$P_{0}(x)=\frac{1}{2^0 0!}\frac{d^0}{d x^0}(x^2-1)^0=1$$
$$P_{1}(x)=\frac{1}{2^1 1!}\frac{d^1}{d x^1}(x^2-1)^1=x$$

より、

$$P_{o}(1)=P_{1}(1)=1\tag{4}$$

である。
\((3^{*})\)式に \(x=1\) を代入し、\(P_{n}(1)=P_{n+1}(1)\) を仮定すると、

$$(n+2)P_{n+2}(1)=(2n+3)P_{n+1}(1)-(n+1) P_{n+1}(1)$$
$$(n+2)P_{n+2}(1)=(n+2)P_{n+1}(1)$$
$$P_{n+2}(1)=P_{n+1}(1)$$

となるので、\((4)\)式の結果を順次代入して、\(P_{n}(1)=1\) が得られる。

同様に、

$$P_{o}(-1)=1,\quad P_{1}(-1)=-1,\quad P_{0}(-1)=-P_{1}(-1)\tag{5}$$

なので、\((3^{*})\)式に \(x=-1\) を代入し、\(P_{n+1}(-1)=-P_{n}(-1)\) を仮定すると、

$$(n+2)P_{n+2}(-1)=(2n+3)\cdot(-1)\cdot(-P_{n}(-1))-(n+1)P_{n}(-1)$$
$$(n+2)P_{n+2}(-1)=(n+2)P_{n}(-1)$$
$$P_{n+2}(-1)=P_{n}(-1)$$

となるので、\((5)\)式の結果を順次代入して、\(P_{n}(1)=(-1)^n\) が得られる。

性質 4.

性質 5., 6.

Neural Trojanミニレビュー

この記事は、2020年1月12日に行われた「第二回サイバーセキュリティ系LT会 in 東京」で発表した、「Neural Trojan — mini review」の内容を解説する記事です。
主に、arXivに投稿されているNeural Trojanについての論文をまとめています。

LT資料

[slideshare id=219010240&doc=seclt-dist20200112-200112075317]

Neural Trojanとは?

「ニューラル トロージャン」と発音します。
簡単に言うと、ニューラルネットワークの中に組み込まれたトロイの木馬のことです。

Liu et al.では、

IP(知財)提供者によってニューラルネットワークの中に組み込まれた、悪意ある隠し機能

(拙訳)

と定義されました。
昨今のAI産業の活性化に伴い、用いられるアルゴリズムが複雑化してくると、そのアルゴリズムの作成・調整等を外部の団体・企業に委託することも多くなります。
その際に、委託先が、顧客の意図しない動作をアルゴリズムの中に組み込むという攻撃が起こりうるのです。

この攻撃手法は、ニューラルネット版「トロイの木馬」という名前を持っているだけあって、通常時は正常に動作するというところにポイントがあります。
すなわち、攻撃者は好きなタイミングでトリガーを発動することによってアルゴリズムに誤作動を起こすことができますが、防御側が事前にその存在を検知することは非常に困難なのです。

Neural Trojanによる攻撃

Liu et al. “Neural Trojans”

Liu et al.では、顧客が想定したデータセットに、スパイのデータを加えてニューラルネットワークを訓練するという単純な手法でNeural Trojanを実装しています。

例えば、防犯カメラのセキュリティ認証システムを考えてみましょう。
顧客は、入館を許可されたメンバーの顔写真を用意し、委託先にニューラルネットワークの訓練を依頼します。
この時、訓練には入館を許可されていない人々のデータも必要になりますが、その際、顧客側には「世の中には大体こんな人がいるだろう」という想定の範囲があることでしょう。
そして後日、訓練済みのネットワークが届き、顧客はその想定の範囲内のデータで性能を検証します。
その際の性能が基準を満たしていたら、そのアルゴリズムは無事導入されることでしょう。

ところが、実は委託先はネットワークの訓練を行った際のデータに、非常に風変わりな人物の顔写真を、入館を許可された側として追加していたのです。
しかし、顧客はこれに気付くことはできません。
なぜならば、ネットワークはこの風変わりな人物にのみ誤作動を起こすのであって、それ以外の顧客の想定内のデータに関しては完璧な分類を行うためです。
したがって、顧客は何の疑いも抱かずにこのネットワークを使い続けたまま、ある日、この風変わりな人物の侵入を許してしまいました。

この風変わりな人物のことを、Neural Trojanの文脈の中ではトリガーと呼びます。

Gu et al. “BadNets: Identifying Vulnerablities in the Machine Learning Model Supply Chain”

Liu et al.では、そもそもが風変わりな画像をトリガーとしていました。
これに対し、Gu et al.では、通常の画像に特徴的なパーツを追加することで、Neural Trojanを作動させました。

このネットワークでは、通常時は画像(例えば、交通標識)の全体を見て、その内容を判断します。
しかし、その画像内に特徴的なパーツ(黄色の四角形、爆弾のマーク、花のマークなど)が追加されていると、このネットワークはその部分にのみ注目して、パーツと関連した内容を出力します。

このアルゴリズムが自動運転の車に積み込まれていたら、どうなるでしょうか?
通常時は車載カメラが標識を正確に認識し、運転を行っています。
しかし、特定のマークが貼られた標識を認めた途端、実際の標識の種類とは関係なく、あらかじめ仕組まれた内容の標識であると判断します。
すなわち、「止まれ」であろうと「進入禁止」であろうと関係なく、「速度制限80km」だと認識させてしまうことができるのです。
これにより、街中の交通標識に黄色いポストイットを貼って回るという行為により、大規模なテロを行うことが可能になります。
大きな標識にポストイットが1つ貼ってあるくらいなら、誰も不審に思わないかもしれません。

また、Gu et al.では、ネットワークを転移学習した後のNeural Trojanの効果についても検証しました。

転移学習とは、とある課題のために訓練されたネットワークを、他の似た課題のために訓練しなおして利用することを指します。
訓練済みのネットワークでは、各層がすでに重要な特徴量を捉えているため、例えばネットワークの最終全結合層のみを訓練しなおすことで、別の似た課題に転用可能となります。
これにより、ゼロの状態からネットワークを構築するよりも相当簡単に、学習モデルを利用することができます。

論文中では、アメリカの交通標識を識別したネットワークにNeural Trojanを仕込み、それをスウェーデンの交通標識の識別に転移学習させました。
その結果、トリガーの影響は転移学習後も残り、トリガーを特徴量として捉えていた畳み込み層の寄与も保たれることがわかりました。

すなわち、「良さそうだ」と思って取ってきたネットワークモデルにNeural Trojanが仕込まれていた場合、たとえある程度の再学習を行ったとしても、安心はできないということなのです。

Clements et al. “Hardware Trojan Attacks on Neural Networks”

Neural Trojanを仕込むことができるのは、ネットワークモデルを作成・訓練するソフトウェア企業だけではありません。
例えば、そのネットワークを載せるハードウェアを開発している企業も攻撃者となり得ます。

Clements et al.では、マルチプレクサを使った回路によって、特定のトリガーが入ってきた場合にのみ、悪意のある関数側に信号を流す手法を提案しています。

Zou et al. “PoTrojan: powerful neuron-level trojan designs in deep learning models”

同様に、ハードウェアにNeural Trojanを仕込む手法として、Zou et al.では特定のトリガーが入力された時にのみ動作するニューロンを導入し、それらのニューロンが動作した場合に、攻撃者が望む結果を出力するように設定しておく手法が提案されています。
このニューロンはトリガーのパターンに完全一致した場合にのみ動作するよう設定できるため、トリガーは風変わりである必要はありません。
一般的な画像であっても、被写体の位置や角度が完全一致していないと動作しないように設定することができます)

Li et al. “Hu-Fu: Hardware and Software Collaborative Attack Framework against Neural Networks”

ハードウェアにNeural Trojanを仕込むLi et al.の手法は、より凝った実装であると言うことができるかもしれません。

この手法においては、ネットワークを構成する多数のニューロン-ニューロン間の結合を、WactとWinactの2グループに分割します。
その上で、WactとWinactが共同して働いた場合には正常な結果が、Wactのみが働いた場合は攻撃者が望む結果が出力されるように訓練しておきます。
そして、マルチプレクサを用いた回路設計により、トリガーが入ってきた時にWinactが停止する(トリガーなし→Winactの計算結果を下流に通す、トリガーあり→下流には0を通す)ように設計しておくと、トリガーによって出力を変えることができるのです。

画像分類モデル以外への適応

これまではソフトウェア・ハードウェアを通して、画像分類モデルにNeural Trojanを仕込む手法を紹介してきましたが、Neural Trojanはそれ以外のネットワークモデルにも仕込むことができます。

例えば、Dai et al.ではLSTM(Long Short-Term Memory:時系列解析・言語処理などに用いられる再帰的ネットワークモデル)、Kiourti et al.では深層強化学習にNeural Trojanを仕込む手法について提案しています。

Neural Trojanへの防御

Neural Trojanを検知することは基本的に困難なのですが、すでに多数の防御手法が提案されています。
ここでは、Neural Trojanの攻撃手法を

  1. 想定外の訓練データを追加する手法(Liu et al.)
  2. トリガーパーツを追加する手法(Gu et al.)

の2種類に分けて、それぞれに対する代表的な防御手法を紹介します。

Liu et al. “Neural Trojans”

想定外の訓練データを追加する攻撃手法に対する防御法は、同じ論文の中で提案されています。
解説のため、もう一度攻撃手法についてのスライドを表示しておきました。

筆者らは、以下の3つの防御手法を提案しています。

  1. 入力されるデータに異常検知を施し、想定外のデータを除く
  2. ネットワークを再訓練する
  3. 入力されるデータに変換を施す(Input Processing)

ここで、1.と2.はやや自明であるとして、3.のInput Processingについて解説します。

Input Processingでは、Neural Trojanが仕込まれた(かもしれない)ネットワークにデータを入力する前に、そのデータを簡単なオートエンコーダを通します。
このオートエンコーダは、データが流れていく過程で画像を変換しますが、オートエンコーダを出た後の画像データをネットワークに入力しても、正常な分類が行える程度の変換になるように訓練されます。

ここで、オートエンコーダの訓練も、想定の範囲内の画像にしか行われないことがポイントになります。
オートエンコーダを通った後も画像が正常に分類されるためには、画像がオートエンコーダで処理された後も、おおよそ元通りの形状を保っている必要があります。
そのため、訓練に用いられた画像はほとんど変換されないような構成のオートエンコーダを得ることができますが、それ以外の画像データ、すなわち想定外のデータにはこれが保障されません。
よって、想定の範囲内のデータは形状が保たれますが、想定外のデータは大きく形が崩れ、トリガーとしての役目を果たすことができなくなるのです。

Chou et al. “SentiNet: Detecting Physical Attacks Against Deep Learning Systems”

トリガーパーツを追加する攻撃手法は、ネットワークをトリガーパーツに注目させて、問題を引き起こすことを意図しています。
そのため、Chou et al.では、ネットワークがどこを見ているかを可視化するGrad-CAMという技術を応用した防御手法を提案しています。

Grad-CAMは深層学習の判断根拠を可視化する手法です。
「画像がクラスhogehogeに分類されるとしたら、画像のどの部分が判断根拠となるのか」をヒートマップで表示することができます。
Chou et al.で提案されたSentiNetという防御手法は、この技術を活用して、以下の手順でトリガーとなったパーツを特定します。

  1. GradCAMを用いてトリガーパーツの候補を抽出
  2. 他の画像に、トリガーパーツの候補を貼ってみる
  3. 同時に、同じ画像のトリガーパーツの候補を貼った位置を隠したもの(マスク)も用意する
  4. トリガーパーツの候補とマスクによって、判断結果が変えられるかを検証
  5. マスクでは結果が変わらないが、トリガーパーツの候補によって変えられる場合、その候補がトリガーパーツであると判断する

以下、1つずつ見ていきます。

1. GradCAMを用いてトリガーパーツの候補を抽出

本当はクラスAなのに、クラスBと誤分類されてしまった画像があるとします。
GradCAMを用いると、それぞれのクラスに対する判断根拠が得られます。
ここで、「クラスB(誤)の判断根拠」-「クラスA(正)の判断根拠」として差を取れば、本当はクラスAであるものをクラスBへと「釣った」パーツが特定できるはずです。
これをトリガーパーツの候補として保持し、他の誤分類された画像にも同様の作業を行えば、複数の候補が得られることになります。

2. 他の画像に、トリガーパーツの候補を貼る

以下、得られたトリガーパーツの候補を1つずつ検証します。
そのために、本来正常に分類が行えるはずの画像群に、これらの候補を1つずつ貼ってみて、分類結果が変えられるかを検証します。

3. 同じ位置にマスクを貼ったものも用意する

とある画像にトリガーパーツの候補を貼った場合、同じ画像の同じ位置をただ隠した画像も用意します。
この行為を、ここでは「マスクを貼った」と表現します。
これは、コントロール群として用いられる画像です。

4. トリガーパーツの候補とマスクの影響を検証

トリガーパーツの候補とマスクのそれぞれを貼った画像を、ネットワークに入力して結果を見ます。
候補が本物のトリガーパーツであるならば、マスク(コントロール)との間には結果に差があるはずです。

5. トリガーパーツの判定

基本的に、パーツを貼り付けることによって画像の分類結果を変えられたなら、そのパーツはトリガーである可能性が高いでしょう。
しかしそれは、ネットワークがそのパーツに注目したことによって結果が変化したのか、パーツを貼り付けたことによって、画像の重要な部分が隠れてしまったことによるものなのかの判断はできません。
コントロール(マスク)は、この後者の影響を検証するために用意されました。

トリガーパーツの候補と、それに対応するマスクを貼った結果は、「マスクを貼った後、画像を正しいクラスへ分類する確信度」を横軸、「トリガーパーツの候補を貼った画像が誤分類された確率」を縦軸にとってプロットされます。
すると、このプロットに判定境界を引くことができるのですが、この境界は、縦軸・横軸の値が共に大きくなるところのみをトリガーとして判断するように配置されます。
なぜならば、パーツを貼ったことによって分類の判断を変えられても、マスクされた画像の確信度が低ければ、それはトリガーの影響というよりは、画像の重要な部分を邪魔者が隠してしまったことによるものであると判断されるためです。

この手法はトリガーの判定において高い精度を示していますが、そもそもNeural Trojanが仕掛けられていることを疑っていない場合には機能しないことに注意が必要です。

まとめ

広告:Geospatial Hackers Program 2019

今までの流れとは全く関係ありませんが、LT会では広告枠があったので載せました。

総務省が主催している地理空間情報活用ハッカソン「Geospatial Hackers Program 2019」にチューターとして参加し、約1か月をかけて全国4会場を巡ります。
各会場の日程は以下の通りです。

  • 愛知会場:2020年02月01日(土)~2020年02月02日(日)
  • 富山会場:2020年02月08日(土)~2020年02月09日(日)
  • 東京会場:2020年02月15日(土)~2020年02月16日(日)
  • 沖縄会場:2020年02月22日(土)~2020年02月23日(日)

この中で、筆者のホームである富山会場(富山県魚津市)では、ハンズオンの講師も務めます。
富山会場は少し特殊で、唯一「Unityを用いた地理空間情報ゲーム開発」をテーマにしています。
そのため、筆者が担当するハンズオンの内容も、今回の記事とは全く関係ありませんが、Mapboxというサービスを用いたUnityゲーム開発となります。

これは、富山会場のハッカソンが、富山県魚津市のゲームクリエイター育成イベント「つくるUOZU」と連携して行われるためです。

また、富山会場が他と異なる点は、2日間のイベントが泊まり込みで行われることです。
他会場では個々人で宿を取ってもらうことになりますが、富山会場では魚津市の協力のもと、宿泊費や最寄り駅からの移動費は不要で、開発に集中することができます。

富山は北陸新幹線が開通して、東京から片道2時間程度と非常に近くなりましたので、

  • Unityに興味がある人・自信がある人
  • 地理空間情報活用に興味がある人

は、ぜひ参加してみてください!

もちろん、他会場もよろしくお願いします。
(この記事の投稿時点で、愛知・東京会場は満席に達しています)

参加登録はこちらから↓
https://ghp.connpass.com/

【Pythonで異常検知】Chapter 1. 1変数正規分布に基づく異常検知

概要

この章では、以下の手順にしたがって1変数データの異常検知をPythonで実践することを目標とする。

  1. 訓練用データを正規分布にフィッティングする
  2. 得られた正規分布からテスト用データの異常度を求める
  3. 異常度の閾値を設定し、それを上回るデータを異常と判定する

データの準備

異常検知を実践するために、今回はScikit-learnのbrest_cancerのデータセットを用いる。このデータセットは、乳房にできた腫瘍の様々な特徴量の値と、その腫瘍が悪性・良性のいずれであったかを記載したラベルから成る。

変数dataにデータセットを読み込むと、data[“data”]に特徴量の値、data[“target”]に特徴量と関連付けられたラベルが格納される。ラベルは0が悪性、1が良性を意味し、特徴量の値からこれらのラベルを予測することが課題となる。

このデータセットは30の特徴量を持つため、今回はその中からworst areaの特徴量だけを用いる。また、このデータセットを訓練用とテスト用に分割することになるが、訓練用には良性(正常)のデータのみを用いるため、良性212例を訓練用としてモデルを作成し、残りの良性145例と悪性212例に対するラベルの予測を行うこととする。

正規分布へのフィッティング

訓練用データの特徴量のばらつきを、ヒストグラムで表示すると以下のようになる。

最尤推定法による正規分布へのフィッティング

特徴量の分布は500付近で多く、そこから遠ざかるにつれて数が少なくなる山型の分布を取っている。今回の異常検知では、この特徴量の確率分布が正規分布にしたがうと考えて、この山型が最も適合する正規分布の形へフィッティングする。このとき、得られるのはフィッティングした正規分布の平均\(\mu\)と標準偏差\(\sigma\)である。なお、正規分布の導出とその性質についての解説は以下の記事を参照のこと。

エントロピーの最大化による正規分布の導出

ここでは、Scipyに実装されているnorm.fitを用いてフィッティングを行う。

上図のオレンジ色の曲線が得られた正規分布である。なお、上図では特徴量の分布も確率密度として表示し直している。フィッティングの結果\(\mu\simeq 549, \sigma\simeq 158\)が得られ、このことから、特徴量の値が500程度になることはよくあることだが、100や1,000になることはほとんどあり得ないことが読み取れる。

確率変数と確率密度関数

異常度の計算

この、特徴量の値の起こりやすさを数値的に表現したものが「異常度」である。1変数正規分布では、特徴量が値\(x\)を取ることについての異常度\(a(x)\)は次のように表される。

$$a(x)=(\frac{x-\mu}{\sigma})^2$$

異常検知における異常度の定義

この値は\(x\)が\(\mu\)に近いときに小さく、\(\mu\)から遠いときに大きい。すなわち、名前の通り異常さの度合いを表している。また、データ数(今回は症例数)が非常に大きい時は、この異常度\(a(x)\)自身が自由度1、スケール因子1の\(\chi^2\)分布に従う。

$$a(x)\sim \chi^2(1,1)$$

異常検知における異常度が、カイ二乗分布に従うことの証明

この\(\chi^2\)分布を表示したのが下図である。

閾値の設定

上図に表示した\(\chi^2\)分布は確率密度関数であるため、これを\((0,\infty)\)の区間で積分すると1になる。また、図によると異常度が0~4程度の範囲が大部分の面積を占めているため、異常度が4以上となる確率はほとんどないことがわかる。異常検知の最終ステップでは、「異常度が〇以上になることは、△%程度の確率でしか起こり得ない(すなわち、ほとんど起こり得ない=異常である)」という異常度の閾値を設定し、それを超える異常度を持つデータ(症例)を異常(悪性)であると判断する。

たとえば、起こり得る確率が5%未満となる異常度の値は次のように求められる。

\(th_5\simeq 3.84\)より、異常度が3.84を超える確率は5%程度である。正常例の5%程度しか取り得ない値の特徴量であればそれはもう異常であると判断してもよく、ここでは上図の赤色で塗りつぶした範囲の異常度をもつデータはすべて異常とするモデルを作成する。

予測結果
悪性良性
真の状態悪性17735
良性7138

predict関数はデータ、パラメータ(\(\mu\)と\(\sigma\))、閾値を引数に与えることで、ラベルの予測結果を返す。また、sklearn.metrics.confusion_matrixは予測結果と真の状態の混同行列を返す。

混同行列の見方とその指標

今回は悪性例212例のうち177例が正しく悪性と判断され、35例の悪性は見落とされた。

ここで、閾値を変えることによって結果を改善することができる場合がある。次の実行例では、異常度の閾値を10%未満と、異常の範囲を広く取っている。

予測結果
悪性良性
真の状態悪性18626
良性19126

その結果、正しく悪性と判定される割合が増加したが、今度は先程よりも多くの良性例が間違って悪性と判断されていることがわかる。

次は逆に、異常度の閾値を1%未満とし、異常の範囲を狭めてみる。

予測結果
悪性良性
真の状態悪性16151
良性1144

今度は間違って悪性と判定される割合が非常に少なくなったが、かなり多くの悪性例を見落としている。

閾値をどのように定めるかによって予測性能は大きく変化しうるが、どの閾値が最適であるかの評価は難しい。最終的には異常検知を行う理由こそが閾値を決定するものとなり、たとえば今回の乳癌の判定のように異常(悪性)を見落とすと大きな不利益を被る場合には、多少偽陽性が増加しても異常の範囲は広く取った方がよいと考えられる。

情報量の定義とエントロピー

情報理論はその名の通り、情報の数量的構造を論ずる学問である。情報を学問として扱うためには、それを量として表すことができる指標を定義する必要がある。この記事では、情報を量的に扱うための指標である情報量について、それが満足すべき特徴から定義を導出し、さらに、それを一般的な状況に適用するためのエントロピーという指標についても定義する。

情報量

考慮すべき特徴

情報を数量的に扱うためには「情報の量」を定義する必要があり、その定義には、一般的に理解されているような情報の性質が組み込まれなければならない。

具体的に議論を進めるために、以下のような状況を設定する。

東地区西地区
北通りAD
中通りBE
南通りCF

とある市の美術館で絵の盗難事件が起き、犯人である怪盗Xは現在、この市のどこかに潜伏していることがわかっている。この市は東と西の2つの地区に分けられ、それぞれの地区には北・中・南の通りが存在する。つまり、この市には上の表のA~Fまでの6つの地域があるのだが、怪盗Xはこのうちのどこかに潜伏している。

この後、この市の警察署にさまざまな情報が寄せられ、それらを考慮した結果、怪盗Xの逮捕に至るのであるが、寄せられたそれぞれの情報には有効性に違いがある。この「有効性」のことを情報量と呼ばれる指標で表現することにしよう。

もし、この市が6つではなく60の地域に分けられていた場合、それらのうち1つに怪盗Xがいることを特定できる情報の価値は、6つの地域から1つを特定する場合よりも大きくなるだろう。よって、 \(n\) 個のものから1つを特定するときの情報量を \(f(n)\) と書くことにする(情報量は \(n\) の関数であり、 \(n\) が増えるほど大きくなる)。

ところで、今回の逮捕の決め手になったのは次の2つの情報であった。

  • \(\alpha\) :「怪盗Xは東地区にいる」
  • \(\beta\) :「怪盗Xは中通りにいる」

これらの情報を組み合わせた結果、怪盗XがBの地域にいることが確定した。さて、先程の議論によれば、 \(\alpha\) の情報は2つの地区から1つ、 \(\beta\) の情報は3つの通りから1つを特定したのであるから、それぞれの情報量は \(f(2),f(3)\) と書ける。そして、これらの情報を合わせたことにより6つの地域から1つを特定できたのであるから、以下の関係式が成り立つはずである。

$$f(6)=f(2)+f(3)$$

これを一般化すると、情報量 \(f(n)\) は

$$f(mn)=f(m)+f(n) \tag{1}$$

という性質を持つ。

定義

上記の特徴を考慮すると、情報量の定義は自然に定めることができる。

\((1)\) 式より

$$f(x+\epsilon x)=f((1+\epsilon)x)=f(1+\epsilon)+f(x)$$

が成り立つ。これを変形すると

$$\frac{f(x+\epsilon)-f(x)}{\epsilon x}=\frac{1}{x}\frac{f(1+\epsilon)}{\epsilon}$$

となり、 \(\epsilon\to 0\) の極限を取ると、左辺は明らかに導関数 \(f'(x)\) となる。右辺の極限値を

$$\lim_{\epsilon\to 0}\frac{f(1+\epsilon)}{\epsilon}=c$$

と定めると

$$f'(x)=\frac{c}{x}$$

となる。すなわち

$$f(x)=\int f'(x)dx=c\log_{e}x+d$$

である。ここで、 \(d\) は積分定数である。この値は、 \(f(1)=0\) であることを考慮すると

$$f(1)=c\log_{e}1+d=d=0$$

より、0となることがわかる。

さて、次は \(c\) の値であるが、これは情報の基本単位をどのように定めるかによって変化しうる。最も一般的な決定方法としては、2つのものから1つを選ぶ「2者択一」の情報量を1と定める、すなわち、 \(f(2)=1\) として \(c\) を計算する。このとき、 \(c=\log_{2}e\) となるため、情報量は

$$f(x)=\log_{2}e\cdot\log_{e}x$$

$$=\log_{2}e\frac{\log_{2}x}{\log_{2}e}=\log_{2}x$$

と表せる。その他の方法としては、 \(f(x)\) の対数の底が \(10\) 、または自然対数 \(e\) となるように \(c\) を定めることが考えられる。 \(f(x)=\log_{2}x\) となるように定めた「2者択一」の情報量の単位をビットというのに対し、 \(f(x)=\log_{10}x\) となるように定めたときの単位をディットといい、実際に情報量の計算を行う際に便利な定め方である。また、 \(f(x)=\log_{e}x\) となるように定めたときの単位をニットといい、こちらは理論計算を行う際に便利である。このように、用途によって対数の底を定めることができるが、これは情報量を測る単位を変えたに過ぎないので、本質的には同じものである。なお、この記事で今後特に記載のない場合は、情報量の単位にビット(底に \(2\) )を用いることにする。

ここで再び怪盗Xの例に戻り、より一般化した状況を考えてみよう。すなわち、とある市は \(k\) 個の地区 \(K_1, K_2, \cdots, K_k\) に分けられ、それぞれが \(l\) 本の通り \(L_1, L_2, \cdots, L_l\) を持ち、全体として \(n=kl\) 個の地域が存在する状況を考える。このとき、怪盗Xが潜伏している確率がすべての地域において等しいとすると、怪盗Xが特定の地区 \(K_i\) にいる確率 \(p(K_i)\) は

$$p(K_i)=\frac{l}{n}$$

と表せる。この、地区を特定する情報の情報量を \(I\) とすると、この情報にあと「 \(l\) 本の通りから1つを特定する」情報量が加わると、「 \(n\) 個の地域から1つを特定する」ことができるので、以下の式が成り立つ。

$$\log n=I+\log l$$

これを変形して

$$I=\log n – \log l=\log\frac{n}{l}=-\log\frac{l}{n}=-\log p(K_i)$$

より、確率 \(p\) で生じる現象を特定する情報量は

$$-\log p$$

で表せることがわかる。以上をまとめると次の情報量の定義が成り立つ。

確率 \(p\) の事象が実際に生起したことを知らせる情報に含まれている情報量を

$$-\log_{2}p ビット$$

と定義する。

すなわち、とある現象が起こりにくいと強く信じられているほど、それが実際に起こったと知らされた時の情報量は大きくなる。情報量を「驚きの程度」と捉えれば、これは日常的な感覚に即した定義である。

エントロピー

概念

これまで怪盗Xの例に登場してきた情報は、すべて「完全に正確な」情報であった。すなわち、「東地区にいる」という情報が寄せられた場合、怪盗Xは必ず東地区にいたのであるが、現実的には早とちりや曖昧な情報が含まれ、情報から直接居場所を特定できることは少ない。

そして、そのような状況設定のもとで定義された情報量 \(-\log p\) は、確率 \(p\) で起きる現象が実際に起こったときの情報量のことであって、「確率 \(p\) で起きる現象が起こる(だろう)」という情報そのものの情報量を表しているのではない。そのため、こうした情報そのものの情報量を考える際には、エントロピーという新たな指標を導入しなければならない。

定義

エントロピーを定義するために、 \(n\) 個の事象 \(A_1,A_2,\cdots,A_n\) が、それぞれ \(p_1,p_2,\cdots,p_n\) の確率で生じる状況を考える。このとき

$$\sum_{i=1}^{n}p_i=1$$

である。ここで、もし \(A_1\) が生じたとすると、そのことを確認したことで得られる情報量は \(-\log p_1\) である。同様に、 \(A_2\) が生じたときは \(-\log p_2\) 、 \(A_i\) が生じたときは \(-\log p_i\) の情報量が得られる。したがって、この状況下で得られる情報量の期待値は、現象ごとに得られる情報量にその確率をかけたものの総和として求められ

$$-\sum_{i=1}^{n}p_i\log p_i$$

と書ける。この期待値のことをエントロピーと定義する。

\(n\) 個の事象 \(A_1,A_2,\cdots,A_n\) が、それぞれ \(p_1,p_2,\cdots,p_n\) の確率で生じる状況におけるエントロピー \(H\) を

$$H(p_1,p_2,\cdots,p_n)=-\sum_{i=1}^{n}p_i\log p_i$$

と定義する。

エントロピーは一般に、「不確定さの度合い」として考えることができる。すなわち、 \(n\) 個の事象のうち、どれが起こるか全く予想がつかない場合では、どの現象が生じてもかなり驚きがある(得られる情報量の期待値が高い)と考えられるのに対し、起こる現象に大体目星がついている状態では、およそ予想通りの現象がおこるため驚きも少ないだろうということである。

性質

不確定さの度合いを表すエントロピーには、以下のような性質がある。

その1

エントロピー \(H\) は非負

$$H \geq 0$$

であり、

$$H=0$$

が成立するのは、どれか1つの \(p_i\) が \(1\) であり、その他はすべて \(0\) のときに限られる。

<証明>

任意の \(p_i\) について

$$1 \geq p_i \geq 0$$

より

$$-p_i\log p_i \geq 0$$

である。したがって

$$H \geq 0$$

が導かれる。以上の議論より、 \(H=0\) となるのは、すべての \(i\) について

$$-p_i\log p_i=0$$

となるときであり、このとき、 \(p_i=0\) または \(p_i=1\) である。ここで

$$\sum_{i}p_i=1$$

より、 \(p_i\) は1つのみが \(1\) で、その他は \(0\) になることがわかる。

その2

\(n\) 個の事象が表すエントロピーの最大値 \(H(n)\) は

$$H(n)=\log n$$

で、これは、すべての事象が等しい確率

$$p_i=\frac{1}{n}$$

で起こるときの不確定度である。

<証明>

$$\sum_{i}p_i=1$$

の条件のもとで

$$H=-\sum_{i}p_i\log p_i$$

を最大にする \(p_i\) を求めることを考える。よって、ラグランジュの未定乗数を \(\lambda\) とし、

$$\Phi=H-\lambda\sum_{i}p_i=-\sum_{i}p_i\log p_i-\lambda\sum_{i}p_i$$

とおく。これを各 \(p_i\) で微分して、それらが \(0\) となるよう各 \(p_i\) を定めればよい。 \(\Phi\) を \(p_1\) で微分すると

$$\frac{\partial \Phi}{\partial p_1}=-\log p_i-log e-\lambda=0$$

※なお、ここで \(log_{2}p_i=\frac{\log_{e}p_i}{log_{e}2}\) より、 \(p_i(\log p_i)’=p_i\cdot\frac{1}{p_i}\cdot\frac{1}{log_{e}2}=\log_{2}e\) となることを用いた。

より

$$\log p_1=-\log e-\lambda$$

$$\log p_1=-\log e-\log 2^\lambda$$

$$\log p_1=\log \frac{1}{e}2^{-\lambda}$$

$$p_1=\frac{1}{e}2^{-\lambda}$$

となる。他の \(p_i\) についても同様で、すべての \(p_i\) について等しい確率値が得られる。よって、 \(\sum_{i}p_i=1\) から

$$p_1=p_2=\cdots=p_n=\frac{1}{n}$$

が導かれる。このときのエントロピーを求めると

$$H=-\sum_{i=1}^{n}\frac{1}{n}\log\frac{1}{n}=\log n$$

となる。このように、すべての確率が等しくなるのは、どの事象が起こるかについて目星がまったくついていない状況を表しており、この状況において不確定度が最大といえるのはもっともな話である。

<参考>

今回は事象が離散的な、すなわち \(n\) が有限個の場合のエントロピーの最大化について考えたが、連続型の確率変数 \(x\) に対する確率密度関数 \(p(x)\) によるエントロピー

$$H[p]\equiv\int_{-\infty}^{\infty}p(x)\{-\ln p(x)\}dx$$

が最大になるのは、 \(p(x)\) が正規分布

$$\mathcal{N}(x|\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\{-\frac{1}{2\sigma^2}(x-\mu)^2\}$$

となるときである。

エントロピーの最大化による正規分布の導出

情報量とエントロピーの関係

「情報が持つ情報量」の定義

性質その2で見た通り、すべての事象が生じる確率が等しい時にエントロピーは最大になる。しかし、ここに「現象 \(A_i\) が起こりやすそう」などという情報が加わることで \(p_i\) が上昇し、その他の確率は低下した結果、エントロピーは減少する。このように、情報はエントロピー、すなわち不確定度を減少させる働きを持ち、このときの現象の度合いによって、その情報が持つ情報量を定義することができる。

情報を得ることによって、状況のエントロピーが \(H\) から \(H’\) へ変わるとき、この情報の持つ情報量を

$$I=H-H’$$

とする。

情報量とエントロピーの計算について、以下のような例を示す。

その1

エントロピー \(H\) の状況において、とある情報が得られたことにより、どの現象が起こったかが確定したとき、不確定度は0になるため、エントロピーは

$$H’=0$$

に変化する。したがって、この情報が持つ情報量は

$$I=H-H’=H$$

となり、まさしく \(H\) 自身である。

その2

8頭立ての競馬を考える。ある予想屋が1着の馬を予想したとき、実際にその馬が1着になる確率は50%であるという。この予想屋が「4番の馬が1着」と予想したとき、この情報の情報量はどれほどだろうか。

1着になれる確率がすべての馬で等しいとすると、4番の馬が1着になる確率はもともと、 \(\frac{1}{8}=12.5\%\) である。よって、この状況におけるエントロピー \(H\) は

$$H=-(\frac{1}{8}\log{\frac{1}{8}}+\frac{1}{8}\log{\frac{1}{8}}+\cdots+\frac{1}{8}\log{\frac{1}{8}})$$

$$=-\frac{1}{8}\log{\frac{1}{8}}\times 8=3.0$$

ここに予想屋の「4番の馬が1着」という情報が加わると、1着になる確率は4番の馬だけ \(\frac{1}{2}=50\%\) で、その他の馬は\(\frac{0.5}{7}=\frac{1}{14}\simeq 7\%\) である。よって、この予想後のエントロピー \(H’\) は

$$H’=-(\frac{1}{2}\log{\frac{1}{2}}+\frac{1}{14}\log{\frac{1}{14}}+\frac{1}{14}\log{\frac{1}{14}}+\cdots+\frac{1}{14}\log{\frac{1}{14}})$$

$$=-(\frac{1}{2}\log{\frac{1}{2}}+\frac{1}{14}\log{\frac{1}{14}}\times 7)$$

$$=-\frac{1}{2}(\log{\frac{1}{2}}+\log{\frac{1}{14}})$$

$$\simeq -\frac{1}{2}(-1.0-3.8)=2.4$$

したがって、予想屋の情報が持つ情報量 \(I\) は

$$I=H-H’\simeq 3.0-2.4=0.6ビット$$

となる。

確率変数と確率密度関数

この記事では、確率論で用いられる「確率変数」や「確率密度関数」などの用語について解説する。

確率変数

定義

確率変数とは、確率論において、起こり得る事柄(事象)に割り当てられている数(通常、整数や実数など)を値として取る変数のことである。

確率変数を\(x\)で表す。

例えばサイコロについて考えると、\(x\)は\(x=1,2,3,4,5,6\)の6種類の値を取り得る。なお、事象\(x\)が生じる確率を\(P(x)\)と書くことにすると、

$$P(1)=P(2)=P(3)=P(4)=P(5)=P(6)=\frac{1}{6}$$

である。

また、\(x\)は「世界の人の身長」というような実数値を取ることもできる。この場合、\(x\)の範囲は世界で最も小さい人の身長以上、最も大きい人の身長以下となり、サイコロについて考えた時のように離散的な値ではなく、連続した値を取る。

確率密度関数

\(x\)が離散的な値を取る場合には特定の\(x\)が生じる確率\(P(x)\)を求めることができたが、\(x\)が連続的な値を取る場合には大きな問題が生じることがわかる。例えば、1人の人間を選んだ時に、その人の身長が(\(1nm\)の誤差もなく)ピッタリ\(170cm\)である確率は0であると言わざるを得ない。このように\(x\)が連続値を取る場合には、任意の値の\(x\)について、\(P(x)=0\)が成り立つ。

そのため、実用上は連続値の特定の1点を考えるのではなく、\(x\)がとある範囲の値を取るときの確率を考える。すなわち、「身長が\(170cm\sim 180cm\)となる確率」などである1)確率変数が連続値を取る場合の他の例としては、時計を見た時に、「秒針が(1ナノ秒の誤差もなく)ピッタリ45秒を指している確率」は0だが、「秒針が45秒から60秒の間を指している確率」は\(\frac{1}{4}\)であることなどを考えよ。。この確率\(P[170\leq x\leq 180]\)は、以下のように積分の形で表される。

$$P[170\leq x\leq 180]=\int_{170}^{180}p(x)dx$$

このとき、式に含まれる\(p(x)\)のことを確率密度関数という。

定義

確率変数\(x\)が、\(a\leq x\leq b\)範囲の値をとる確率を\(P[a\leq x\leq b]\)とすると、

$$P[a\leq x\leq b]=\int_{a}^{b}p(x)dx$$

が成り立つとき、\(p(x)\)を確率密度関数という。

性質

定義より、確率変数\(x\)が、\(x\)から\(x+dx\)までの領域の値を取るときの確率は\(p(x)dx\)となる。

また、確率変数\(x\)の定義域を\({\bf R}\)とすると、全体確率は1となることより、

$$\int_{{\bf R}}p(x)dx=1$$

を満たす。上式の条件を、規格化条件という。また、確率密度関数は常に正値となる。したがって、\(p(x)\)が確率密度関数であるための本質的な条件は以下の2点である。

  1. $$\int_{{\bf R}}p(x)dx=1 \tag{1}$$
  2. $$p(x)\geq 0 \tag{2}$$

多変数への応用

これまでは確率変数が1次元の場合を考えてきたが、\(M\)次元確率変数\({\bf x}\)についても

$$P[{\bf a\leq x\leq b}]=\int_{{\bf a}}^{{\bf b}}p({\bf x})d{\bf x}$$

を満たす確率密度関数\(p({\bf x})\)を定義でき、確率変数\({\bf x}\)が、\({\bf x}\)から\({\bf x}+d{\bf x}\)までの領域の値を取る確率は\(p({\bf x})d{\bf x}\)となる。このとき、

$$d{\bf x}=\prod_{i=1}^{M}dx_i$$

である。

\(M=2\)における状況を図示して説明すると、

\(d{\bf x}=dx_1dx_2\)は上図の面積に等しい。これに\(p({\bf x})\)を掛けると確率になることから、\(p({\bf x})\)は単位面積当たりの確率と考えることができる。これが確率「密度」関数と呼ばれる理由である。

\(M\geq 3\)の場合も、\(d{\bf x}=dx_1dx_2\cdots dx_M\)を体積に相当するものと考えることで同様の議論が成り立つ。

確率密度関数の周辺化

\(M\)次元確率変数\({\bf x}=(x_1,x_2,\cdots,x_M)\)に対する確率密度関数\(p({\bf x\})\)について、\(x_3\)から\(x_M\)までのみ積分を行った結果を\(f(x_1,x_2)\)とすると、

$$f(x_1,x_2)=\int_{{\bf R}}p({\bf x})dx_3 dx_4 \cdots dx_M$$

となる。このとき、\(p({\bf x})\geq 0\)であることから、\(f(x_1,x_2)\geq 0\)である。また、\(f(x_1,x_2)\)を\(x_1\)と\(x_2\)について積分すると

$$\int_{{\bf R}}f(x_1,x_2)dx_1dx_2=\int_{{\bf R}}(\int_{{\bf R}}p({\bf x})dx_3 dx_4 \cdots dx_M)dx_1dx_2$$

$$=\int_{{\bf R}}p({\bf x})dx_1dx_2dx_3 dx_4 \cdots dx_M$$

$$=\int_{{\bf R}}p({\bf x})d{\bf x}=1$$

より、\((1),(2)\)式を満たすため、\(f(x_1,x_2)\)は\(2\)次元の確率変数\((x_1,x_2)\)についての確率密度関数となる。このように、\(M\)次元確率変数に対する確率密度関数を不要な\(N\,(N<M)\)次元について積分のみすることで、必要な\(M-N\)次元確率変数に対する確率密度関数を求めることができる。この手法を周辺化といい、得られた確率密度関数(今回の場合は\(f\))のことを周辺分布という。また、もともとの分布のことは、周辺分布と対比的に同時分布または結合分布という。

References   [ + ]

1.確率変数が連続値を取る場合の他の例としては、時計を見た時に、「秒針が(1ナノ秒の誤差もなく)ピッタリ45秒を指している確率」は0だが、「秒針が45秒から60秒の間を指している確率」は\(\frac{1}{4}\)であることなどを考えよ。

最尤推定法による正規分布へのフィッティング

観測された複数のデータがとある分布に基づいていると仮定して、その分布の形状を決定するパラメータを求める際、最尤推定法という手法がよく用いられる。この記事では、観測された結果が正規分布に従うと仮定した際に、最尤推定法を用いて平均 \(\mu\) (または \({\boldsymbol \mu}\) )と分散 \(\sigma^2\) (または共分散行列 \({\boldsymbol \Sigma}\) )を求める方法について解説する。

前提

1次元正規分布

$$\mathcal{N}(x|\mu,\sigma^2) \equiv \frac{1}{(2\pi\sigma^2)^{\frac{1}{2}}}\exp\{-\frac{1}{2\sigma^2}(x-\mu)^2\} \tag{1}$$

多次元正規分布

$$\mathcal{N}({\bf x}|{\bf \mu},{\boldsymbol \Sigma}) \equiv \frac{|{\boldsymbol \Sigma}|^{-\frac{1}{2}}}{(2\pi)^\frac{M}{2}}\exp\{-\frac{1}{2}({\bf x}-{\boldsymbol \mu})^T{\boldsymbol \Sigma}^{-1}({\bf x}-{\boldsymbol \mu})\} \tag{2}$$

尤度関数

最尤推定法は、観測結果 \(\mathcal{D}=\{x^{(1)},x^{(2)},\cdots,x^{(N)}\}\) (多次元の場合は \(\mathcal{D}=\{{\bf x}^{(1)},{\bf x}^{(2)},\cdots,{\bf x}^{(N)}\}\) )が与えられたとき

$$p(\mathcal{D}|{\boldsymbol \theta}) \equiv \prod_{i=1}^{N}p(x^{(i)}|{\boldsymbol \theta})$$

(多次元の場合は

$$p(\mathcal{D}|{\boldsymbol \theta}) \equiv \prod_{i=1}^{N}p({\bf x}^{(i)}|{\boldsymbol \theta})$$ )

で定義される尤度関数 \(p(\mathcal{D}|{\boldsymbol \theta})\) を最大化するようなパラメータ \({\boldsymbol \theta}\) を求める方法である。このとき、尤度関数をそのまま用いるよりも自然対数をとった方が計算に便利なため

$$L({\boldsymbol \theta}|\mathcal{D}) \equiv \ln{p(\mathcal{D}|{\boldsymbol \theta})}$$

で定義される対数尤度関数 \(L({\boldsymbol \theta}|\mathcal{D})\) の最大化を目指すことが多い。

1次元正規分布の最尤推定

最尤推定法により求めるパラメータは、 \({\boldsymbol \theta}=(\mu,\sigma^2)\) であり、尤度関数は

$$p(\mathcal{D}|\mu,\sigma^2) = \prod_{i=1}^{N}\mathcal{N}(x^{(i)}|\mu,\sigma)$$

と表される。ここで、 \((1)\) 式を代入すると

$$p(\mathcal{D}|\mu,\sigma^2) = (2\pi\sigma^2)^{-\frac{N}{2}}\exp[-\frac{1}{2\sigma^2}\{(x^{(1)}-\mu)^2+(x^{(2)}-\mu)^2+\cdots+(x^{(N)}-\mu)^2\}]$$

$$= (2\pi\sigma^2)^{-\frac{N}{2}}\exp\{-\frac{1}{2\sigma^2}\sum_{i=1}^{N}(x^{(i)}-\mu)^2\}$$

となり、対数尤度関数は

$$L(\mu,\sigma^2|\mathcal{D}) = \ln{p(\mathcal{D}|\mu,\sigma^2)} = -\frac{N}{2}\ln{2\pi\sigma^2}-\frac{1}{2\sigma^2}\sum_{i=1}^{N}(x^{(i)}-\mu)^2$$

と表される。これを最大化する \(\mu,\sigma^2\) を求めるためには、 \(\mu,\sigma^2\) でそれぞれ偏微分してその結果を0と置けばよいが、ここでは計算を簡略化するため \(\sigma^2\) ではなく \(\sigma^{-2}\) で偏微分する。したがって

$$\frac{\partial L}{\partial \mu} = \frac{\partial}{\partial \mu}[-\frac{N}{2}\ln{2\pi\sigma^2}-\frac{1}{2\sigma^2}\sum_{i=1}^{N}\{(x^{(i)})^2-2x^{(i)}\mu+\mu^2\}]$$

$$= -\frac{1}{\sigma^2}\sum_{i=1}^{N}(\mu-x^{(i)}) = 0$$

より

$$N\mu = \sum_{i=1}^{N}x^{(i)}$$

$$\mu = \frac{1}{N}\sum_{i=1}^{N}x^{(i)}$$

が導かれ、また

$$\frac{\partial L}{\partial \sigma^{-2}} = \frac{\partial}{\partial \sigma^{-2}}\{-\frac{N}{2}(\ln{2\pi}-\ln{\sigma^{-2}})-\frac{1}{2\sigma^2}\sum_{i=1}^{N}(x^{(i)}-\mu)^2\}$$

$$= \frac{N\sigma^2}{2}-\frac{1}{2}\sum_{i=1}^{N}(x^{(i)}-\mu)^2 = 0$$

より

$$\frac{N\sigma^2}{2} = \frac{1}{2}\sum_{i=1}^{N}(x^{(i)}-\mu)^2$$

$$\sigma^2 = \frac{1}{N}\sum_{i=1}^{N}(x^{(i)}-\mu)^2$$

が導かれる。以上より、1次元正規分布の最尤推定では次の関係式が成り立つ。

観測結果 \(\mathcal{D}=\{x^{(1)},x^{(2)},\cdots,x^{(N)}\}\) が1次元正規分布に従うと仮定した場合、その標本平均 \(\hat{\mu}\) と標本分散 \(\hat{\sigma}\) は、最尤推定法により以下のように求められる。

$$\hat{\mu} = \frac{1}{N}\sum_{i=1}^{N}x^{(i)}$$

$$\hat{\sigma^2} = \frac{1}{N}\sum_{i=1}^{N}(x^{(i)}-\hat{\mu})^2$$

多次元正規分布の最尤推定

作成中。詳細は

井手剛「入門 機械学習による異常検知――Rによる実践ガイド――」(2015)コロナ社

を参照。