jiku log

データサイエンスの核心を掴む : 学びと発見の記録

「確率的機械学習 入門編I」を読む ~第4章 統計学 ③ベイズ統計学・その1~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第4章 統計学」のベイズ統計学における,ベータ・二項モデルに関する読書メモである。

4.6 ベイズ統計学

本節では,推定における不確実性を扱うことができるベイズ統計学について説明している。

「推定の不確実性」を考慮する応用シーンとして,能動学習や過剰適合の回避などが挙げられる。

ベイズ統計学における不確実性の表現

ベイズ統計学において,不確実性は事後分布(posterior distribution)によって表現する。


 \begin{align}
p(\boldsymbol{\theta} \mid \mathcal{D}) = \frac{ p(\boldsymbol{\theta}) p(\mathcal{D} \mid \boldsymbol{\theta})}{ p(\mathcal{D}) }
= \frac{p(\boldsymbol{\theta} ) p(\mathcal{D} \mid \boldsymbol{\theta}) }{ \int p(\boldsymbol{\theta'} ) p(\mathcal{D} \mid \boldsymbol{\theta'}) d\boldsymbol{\theta'}} \\ \\
\end{align}

ベイズ統計学における予測

パラメータの事後分布を計算した後に,入力が与えられた下での事後予測分布は,


 \begin{align}
p(\boldsymbol{y} \mid \boldsymbol{x}, \mathcal{D} ) = \int p(\boldsymbol{y} \mid \boldsymbol{x}, \boldsymbol{\theta}) p(\boldsymbol{\theta} \mid \mathcal{D}) d \boldsymbol{\theta} \\ \\
\end{align}
である。

●ブログ筆者註 :
上記の事後予測分布は,

  1.  p(\boldsymbol{\theta} \mid \mathcal{D})から \boldsymbol{\theta}を1つ生成する。
  2.  \boldsymbol{\theta}を入力として, p(\boldsymbol{y} \mid \boldsymbol{x}, \boldsymbol{\theta}) を構成する。
  3.  \boldsymbol{\theta}について足し合わせる(積分する)。

という手順で得ていると解釈できる。

このような足し合わせをベイズモデル平均と呼ぶ。

4.6.1 共役事前分布

本節では,共役事前分布(conjugate prior)について説明している。

共役事前分布とは,事前分布 p(\boldsymbol{\theta})と事後分布事前分布 p(\boldsymbol{\theta} \mid \mathcal{D})が同じ分布族 \mathcal{F}に属するような事前分布のことである。

4.6.2 ベータ・二項モデル

本節では,共役事前分布の例として,ベータ・二項モデルを紹介している。

ベータ・二項モデルでは,

  • 尤度 : ベルヌーイ分布  p(\mathcal{D} \mid \theta) =  \theta^{N_1} (1 - \theta)^{N_0}
  • 事前分布 : ベータ分布  p(\theta) \propto \mathrm{Beta}(\theta \mid \overset{\smallsmile}{\alpha}, \overset{\smallsmile}{\beta}) \propto \theta^{\overset{\smallsmile}{\alpha}-1} (1-\theta)^{\overset{\smallsmile}{\beta} - 1}

としている(尤度は二項分布でもよい)。

そのため事後分布は,


 \begin{align}
p(\theta \mid \mathcal{D}) 
&\propto \mathrm{Beta}(\theta \mid \overset{\smallsmile}{\alpha} + N_1, \overset{\smallsmile}{\beta} + N_0) \\ \\
&= \mathrm{Beta}(\theta \mid \overset{\smallfrown}{\alpha}, \overset{\smallfrown}{\beta}) \\ \\
\end{align}
となり,再びベータ分布となる。

下図に,事前分布がベータ分布の場合(左)と一様分布(右)の場合を示す。事前分布がベータ分布の場合は,尤度と事前分布のバランスを取るような事後分布が形成される。一方で事前分布が一様分布の場合は,尤度と事後分布が一致する。

事前分布がベータ分布の場合(左)と一様分布の場合(右)

事後分布のモード(MAP推定)

最も可能性が高いパラメータの値を与えるのがMAP推定である。


 \begin{align}
\hat{\theta}_{map} 
&= \operatorname*{argmax}_{\theta} p(\theta \mid \mathcal{D}) \\ \\
&= \operatorname*{argmax}_{\theta} \log p(\theta) + \log p(\mathcal{D} \mid \theta ) \\ \\
\end{align}
一様事前分布 p(\theta) \propto 1を用いると, \log p(\theta) = 0なので,MAP推定は最尤推定(MLE)と等しくなる。

事後平均

事後分布のモードは単一の点であるため,事後分布の特徴が十分に捉えきれていないことがある。モードよりも頑健な推定として,事後平均が挙げられる。


 \begin{align}
\bar{\theta} \equiv E[ \theta \mid \mathcal{D} ] \\ \\
\end{align}

事後予測分布

将来の値を予測したいと考えたとき,訓練データによって推定されたパラメータ \hat{\theta}をモデルに代入した p(y \mid \hat{\theta})を用いた予測をプラグイン近似と呼ぶ。ただし,プラグイン近似は過剰適合に陥りやすい。

これを避けるためのベイズ的な方法は, \thetaを周辺化する事後予測分布が挙げられる。例えばベルヌーイモデルでは,


 \begin{align}
p(y=1 \mid \mathcal{D}) = \int_{0}^{1} p(y=1 \mid \theta) p(\theta \mid \mathcal{D}) d\theta \\ \\
\end{align}
のようにして求める。


 M \gt 1回の硬貨投げの試行を考えると,ベルヌーイモデルの代わりに二項モデルを用いる。この場合,事後予測分布は,


 \begin{align}
p(y \mid \mathcal{D}, M) = \int_{0}^{1} \mathrm{Bin}(y \mid M, \theta) \mathrm{Beta}(\theta \mid \overset{\smallfrown}{\alpha}, \overset{\smallfrown}{\beta}) \\ \\
\end{align}
となるが,この積分を実行すると,ベータ・二項分布

 \begin{align}
Bb(x \mid M, \overset{\smallfrown}{\alpha}, \overset{\smallfrown}{\beta}) \equiv \binom{M}{x} \frac{ B(x+\overset{\smallfrown}{\alpha}, M-x+\overset{\smallfrown}{\beta}) }{ B(\overset{\smallfrown}{\alpha}, \overset{\smallfrown}{\beta}) } \\ \\
\end{align}
が得られる。

まとめと感想

今回は,「4 統計学」のベイズ統計学における,ベータ・二項モデルについてまとめた。

ベータ分布は共役事前分布の代表選手

本書では,ベータ・二項モデルを題材に,事後分布のモード(MAP推定),事後平均,事後分散,事後予測分布など,事後分布から導出される様々な統計量・分布を紹介していた。

事後分布は一般に計算しづらいと言われているが,ベータ・二項モデルは事後分布が解析的に求まるので,様々な性質を理解するうえで有益だと感じた。


本記事を最後まで読んでくださり,どうもありがとうございました。

2025年11月の学びの振り返り

はじめに

当ブログのサブタイトルは,「データサイエンスの核心を掴む : 学びと発見の記録」である。2025年11月の学びと発見を振り返ってみた。

インプット編

2025年11月は,2冊の本を読んだ。

「先輩データサイエンティストからの指南書」

https://gihyo.jp/assets/images/cover/2025/9784297151003.jpg

私は製造業のデータサイエンスチームを率いる立場として,メンバーの育成やプロジェクトの推進に日々取り組んでいる。製造業におけるデータ分析は,単発のモデル構築で終わるのではなく,長期運用・システム連携・品質保証までつなげることが重要である。そのため,チームとして成果を出すには,アルゴリズムの知識だけでなく,エンジニアリングの基礎体力・再現性のある開発環境・品質を保つ仕組みが欠かせないと考える。

浅野純季,木村真也,田中冬馬,武藤克大,栁泉穂 著「先輩データサイエンティストからの指南書」は、まさにその実務基盤を体系的に整理した書籍である。第1章ではデータサイエンスを価値に結びつけるためのエンジニアリング思考を,第2章から第6章では環境構築,コード品質,データ品質,実験管理,そしてプロトタイプ開発といったテーマを段階的に扱っている。

本書を通じて,メンバーが安心して開発・検証・共有できる環境を整えるためのヒントを得たいと考えている。また,属人的な取り組みから脱し,品質とスピードを両立するデータサイエンスチームの在り方について改めて考えるきっかけにしたいと考えている。

stern-bow.hatenablog.com

「Data-centric AI入門」

https://gihyo.jp/assets/images/cover/2025/thumb/TH800_9784297146634.jpg

近年AI・データ活用の成果は,アルゴリズムの高度化よりも,データ品質とデータ運用プロセスの整備によって左右されることが明確になりつつある。特に製造業においては,センサ情報,設備ログ,検査データなど,多種多様なデータが日常的に生成されているが,それらを適切に収集・整理・改善しなければ,AIは十分な価値を発揮できない。

片岡博雄 監修「Data-centric AI入門」は,データを中心としたAI開発の原則と実践を体系的に示した書籍である。データ品質の重要性,画像・テキスト・ロボット分野における具体例など,AIを業務成果へ結びつけるための考え方が整理されている。本書の内容は,AI活用を単発の技術検証で終わらせず,継続的な改善と価値創出につなげるための組織的基盤構築に資するものであると期待している。

本書を通じて,製造現場でのデータ利活用の知見と照らし合わせながら,データ管理方針,AIプロジェクト推進体制,品質保証プロセスに関する示唆を得て,実務への適用可能性を検討していきたい。

stern-bow.hatenablog.com

アウトプット編

ブログの記事

ブログの記事は,30個作成した。今月の記事は,「先輩データサイエンティストからの指南書」,「Data-centric AI入門」および「確率的機械学習 入門編I」の読書メモだった。


今月は,月間アクセス数が9,000件を超えた。またブログ開設からの合計アクセス数が90,000件を超えた。

振り返り

統計検定1級と弊ブログへのアクセス

今年は,2025/11/16に統計検定1級の試験があった。この日の前日にアクセス数が伸びており,弊ブログを受験勉強に使っている方がいるのではないかと推察している。

弊ブログは,最近では最近読んだ書籍の読書メモを記事にすることがほとんどである。一方で,読者の中には統計検定1級合格に向けたヒントを得たい方が多いようである。
過去に学んだ数理統計学のテキストをたまに読み返したり,実務で学んだ手法を使おうとしたときに,これまで理解できていなかったことに理解が及んでいることがある。
このような内容をテーマに記事を作成して,統計検定を受験される方の参考になるようにしていきたい。

実務寄りの本の読書

今月は,数理統計や機械学習の記事というより,実務寄りの本を読んだ。普段の読書メモは,読み進めていて足りない式展開を補足したり,理解がしづらかった部分を補ったりしているが,今月読んだ本はいずれも数式が少なめだったため,各章の要約と,実務に活用していくうえで考察したことなどを書いた。
これまでもビジネス寄りの本を読むことは多かったが,今回のように各章の要約と実務への応用を考えると,これまでよりも理解が深まったように感じた。今後は,実務寄りの本を読む際には,このような読み方を意識していきたい。


本記事を最後まで読んでくださり,どうもありがとうございました。

平均・分散に関するデルタ法 (2025年改訂版) #統計検定

はじめに

この記事は,統計・機械学習の数理 Advent Calendar 2025の1日目の記事である。

統計検定1級の統計数理・理工学や統計検定準1級において,デルタ法はわりとお馴染みの話題である(1級よりも準1級に出やすい印象がある)。デルタ法には,平均に関するものと分散に関するものがあるので,整理してみた。

※この記事は,過去に作成した記事である平均・分散に関するデルタ法 #統計検定 - jiku log に加筆修正したものである。

この記事が対象とする読者

  • 統計検定1級・準1級の受験者
  • 数理統計学に興味がある方

本記事で得られること(ベネフィット)

  • 確率変数の変数変換と分布収束の関係であるデルタ法の概要について説明している。
  • 確率変数の変数変換によって,分散と平均がどのように変換されるか,ということについて説明している。
  • サンプルが増えることによる挙動の変化を,数値実験によって説明している。

デルタ法の概要

定義と雰囲気

デルタ法に関する定義は,久保川達也 著「現代数理統計学の基礎」 (以下,久保川本と記載)の第5.3節に説明されている。

雰囲気だけ説明すると,デルタ法とは「ある確率変数が分布収束するときに,その確率変数を連続微分な関数で変換しても,変数変換後の確率変数も同様に分布収束する」という定理である*1。ポイントとしては,分布収束がかかわってくることである。

統計検定とのかかわり

統計検定において,デルタ法は重要な出題分野の1つである。

統計検定準1級の参考書である,日本統計学会 編「統計検定準1級対応 統計学実践ワークブック」では,いくつかの練習問題とともにデルタ法が説明されている。

また統計検定1級の統計数理・理工学において,デルタ法に関連する問題は,2012年~2022年の10年の間で,

  • 2018年 統計数理 (問1・問2)
  • 2012年 統計数理 (問3)

において出題されている。

デルタ法の説明

おさらい : 中心極限定理

デルタ法の説明をする前に,関連する定理についておさらいする。

分布収束するような確率変数の列に関する定理として代表的なものは,中心極限定理がある。中心極限定理は,「平均と分散を持つ確率変数であれば,その確率変数自体が正規分布にしたがわなくても,標本平均がしたがう確率分布は正規分布に近づく」という定理である。

より詳細には,中心極限定理は,確率変数 X_i (i=1,...,n)の分布関数 Fについて, E[X_i] = \mu, V[X_i ] =\sigma^2 が存在すると仮定したうえで,標本平均  \bar{X_n} = \frac{1}{n} \sum X_iを標準化したものが,標準正規分布に収束するというものである。


 \begin{align}
\frac{\sqrt{n}(\bar{X_n} - \mu)}{\sigma} \rightarrow _d \mathcal{N}(0, 1) \\ \\
\end{align}

またこの式から,


 \begin{align}
\sqrt{n}(\bar{X_n} - \mu) \rightarrow _d \mathcal{N}(0, \sigma^2) \\ \\
\end{align}
が得られる(なお \to _dという記号は,この後説明する分布収束を示す)。

おさらい : 分布収束

中心極限定理は,分布収束の一例である。分布収束は,久保川本では次のように説明されている。

確率変数の列 \{ U_n \} _{n=1,2,...}が確率変数 U分布収束するとは,


 \begin{align}
\lim _{n \to \infty} P(U_n \leq x) = P(U \leq x) = F_U(x) \\ \
\end{align}
が, F_U(x)の連続点で成り立つことをいい, U_n \to_d Uで表す。

「確率変数の列 \{ U_n \} _{n=1,2,...}」の具体例として,標本平均 \bar{X}_nの列が挙げられる。
標本平均において nを大きくするということは,平均の計算に用いるサンプルを増やすことに対応するが,先に説明した中心極限定理により,標本平均がしたがう分布は正規分布に収束する。

中心極限定理とデルタ法

デルタ法は,変数変換と分布収束に関連する定理である。確率変数自体が正規分布にしたがっている必要はないが,デルタ法は中心極限定理との組合せで出てくることが多いので,久保川本では,次のように説明されている。

確率変数の列 \{ U_n \} _{n=1,2,...}と連続微分可能な関数 g(\cdot)について,


 \begin{align}
\sqrt{n}(U_n - \mu) \rightarrow _d \mathcal{N}(0, \sigma^2) \\ \\
\end{align}
が成り立つときには,デルタ法から


 \begin{align}
\sqrt{n}(g(U_n) - g(\mu)) \rightarrow _d \mathcal{N}(0, \sigma^2 \{ g'(\mu) \}^2  ) \\ \\
\end{align}
となる。

デルタ法の証明の流れは以下のようになる。

  •  g(U_n) \mu周りでテイラー展開し, g(U_n) = g(\mu) + g'(\mu^*) ( U_n - \mu)とする。
    • ただし \mu^*は, \lvert \mu^* - \mu \rvert \lt \lvert U_n - \mu \rvert となる点である。
  • スラツキーの定理・連続写像定理から, g'(\mu^*) \to g(\mu)となる。
  • 再びスラツキーの定理から,デルタ法が導かれる。

平均・分散に関するデルタ法の覚え方

統計検定では,デルタ法を用いて変数変換後の平均・分散を評価する問題が出題されることがある。ここでは,デルタ法の覚え方について説明する。

以下では,

  • 確率変数列 \{X_n\}_{n=1,2,...}(たとえば標本平均の列)が構成されており,確率変数 \sqrt{n}(X_n - \mu) \to _d \mathcal{N}(0, \sigma^2)のように分布収束する。
  •  X_nには,平均 \mu_nと分散 \sigma^2_nが存在する。

とする。

分散に関するデルタ法の覚え方 (1次近似)

分散に関するデルタ法は,ある確率変数列 X_nを連続微分可能な関数 Y=g(X_n)で変換したときの分散を求めるための式である。直感的な覚え方は,「テイラー展開の1次近似を用いる」である。

まず Y=g(X_n)を, \mu_nの周りで1次の項までテイラー展開する。


 \begin{align}
Y = g(X_n) &\approx g(\mu_n) + g'(\mu_n)(X_n - \mu_n) \\ \\
\Rightarrow g(X_n) -  g(\mu_n) &\approx  g'(\mu_n)(X_n - \mu_n) \\ \\
\end{align}
なお左辺の「差」は,右辺の線形化によってどれくらい誤差が広がるかを示している。


両辺をそれぞれ2乗して期待値を取ると,


 \begin{align}
E [ (g(X_n) -  g(\mu_n))^2 ] &\approx  E [ \{ g'(\mu)(X_n - \mu_n) \}^2 ] = \{ g'(\mu_n) \}^2  E [ (X_n - \mu_n)^2 ] \\ \\
\end{align}
となる。

よって,


 \begin{align}
V [ g(X_n) ] &\approx  \sigma_n^2 \{ g'(\mu_n) \}^2   \\ \\
\end{align}
となり,関数 Y=g(X_n)で変換したときの分散が求められた。

平均に関するデルタ法の覚え方 (2次近似)

平均に関するデルタ法は,ある確率変数列 X_nを連続微分可能な関数 Y=g(X_n)で変換したときの平均を求めるための式である。直感的な覚え方は,「テイラー展開の2次近似を用いる」である。

 Y=g(X_n)を, \mu_nの周りで2次の項までテイラー展開する。


 \begin{align}
Y = g(X_n) \approx g(\mu_n) + g'(\mu_n)(X_n - \mu_n) + \frac{1}{2} g''(\mu_n) (X_n - \mu_n)^2 \\ \\
\end{align}

両辺の期待値を取ると, E [ X_n - \mu_n ] =0 より,


 \begin{align}
E [ g(X_n) ]  
&\approx E [ g(\mu_n) ]  + \frac{1}{2} E [ g''(\mu_n) (X_n - \mu_n)^2 ] \\ \\
&= g(\mu_n) + \frac{1}{2} \sigma_n^2 g''(\mu_n) \\ \\
\end{align}

となり,関数 Y=g(X_n)で変換したときの平均が求められた。

数値実験

デルタ法は,変数変換後の分布に関する漸近的な挙動を示した定理である。本節では,サンプル数を増やしていったときにおける,変数変換後の挙動を確認する。

問題設定

確率変数 X_i正規分布にしたがい, X_i \sim \mathcal{N}(\mu, \sigma^2)であるとする。またデルタ法で用いる変数変換は Y = g(X) = e^Xという変換であるとする。

サンプルを n個発生させたときの標本平均を \bar{X}_nとし、この標本平均にデルタ法を適用すると,


 \begin{align}
\sqrt{n}(g(\bar{X}_n) - g(\mu)) \rightarrow _d \mathcal{N}(0, \sigma^2 \{ g'(\mu) \}^2  ) \\ \\
\end{align}
となる。これから, g(\bar{X}_n)がしたがう確率分布は,

 \begin{align}
g(\bar{X}_n) \sim  \mathcal{N}(g(\mu), \sigma^2 \{ g'(\mu) \}^2 /n ) \\ \\
\end{align}
と近似できる。

サンプルサイズ nが小さい範囲では, g(\bar{X}_n)の分散とデルタ法による分散 \sigma^2 \{ g'(\mu) \}^2 /nの値が離れているが,サンプルサイズ nが大きくなると,これらの値は近付いていくと考えられる。


今回は以下の手順で,2種類の分散を比較した。

  • Step 1.  \mathcal{N}(\mu, \sigma^2)にしたがう変数を n個生成して,標本平均 \bar{X}_nを算出する。これを繰り返し,標本平均 \bar{X}_n n_{trial}個作成する。
  • Step 2. 標本平均 \bar{X}_nを関数 g(X)で変換した値を生成し( n_{trial}個生成される),標本分散 V[g(X)]を計算する。
  • Step 3. デルタ法による分散 \sigma^2 \{ g'(\mu) \}^2 /nを算出する。
  • Step 4. Step 2. で計算した標本分散 V[g(X)]平方根(標本標準偏差)と,Step 3. で計算したデルタ法による分散 \sigma^2 \{ g'(\mu) \}^2 /n平方根の比を計算する。

サンプル数は n=10, 50, 200, 1000, 10000とした。

実験結果1 : ヒストグラムの比較

まずStep 1. で作成した n_{trial}個の標本平均 \bar{X}_nヒストグラムと,Step 2. で作成した g(\bar{X}_n)ヒストグラムを並べて描いた。

標本平均のヒストグラム(上段)と,g(標本平均)のヒストグラム

サンプルサイズが小さい n=10のときは, g(\bar{X}_n)ヒストグラムは対数正規分布に近いことが確認できた。一方でサンプルサイズが大きくなるにつれて, g(\bar{X}_n)ヒストグラム正規分布に近付くことが確認できた。

考察

対数正規分布は,確率変数 X正規分布にしたがっているときに, \exp(X)がしたがう分布である。
サンプルサイズが小さいときは,標本平均 \bar{X}_nの値のばらつきが大きく,対数正規分布の定義域の広い範囲をカバーするため,ヒストグラムの形状も対数正規分布に近くなったと考えられる。
一方でサンプルサイズが大きいときは,大数の法則により標本平均 \bar{X}_nの値が母平均に近づきばらつきが減るため,ヒストグラムは母平均周りの正規分布のような形になったと考えられる。

実験結果2 : 標準偏差の比

次に,サンプルサイズと,「(デルタ法による標準偏差) / (g(標本平均)の標本標準偏差)」の比率を下図に示す。

(デルタ法による標準偏差)/(g(標本平均)の標本標準偏差)の比率

サンプルサイズが大きくなるにつれて,比率が1に近付くことが確認できた。

まとめ

統計検定1級や準1級でよく出てくるデルタ法について,平均・分散に関するデルタ法を紹介した。
最後に本記事の内容を,チートシートとしてまとめる。

平均・分散のデルタ法

宣伝

統計・機械学習の数理 Advent Calendar 2025 について

こちらのAdvent Calenderは,統計・機械学習に関する話題を扱っている。興味深いテーマが公開されていくので,ぜひご覧ください。
adventar.org

弊ブログについて

弊ブログは,数理統計学機械学習,データサイエンスについて学んだことを記事にまとめており,

といったコンテンツを公開しているので,他の記事もご覧ください。

この記事が面白いと感じていただけたら,こちらのバナーを押していってください↓

本記事を最後まで読んでくださり,どうもありがとうございました。

*1:「デルタ法」という名前だと何やら計算手法っぽい印象を受けるが,久保川本では「定理」と説明されている

「確率的機械学習 入門編I」を読む ~第4章 統計学 ②経験リスク最小化~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第4章 統計学」における経験リスク最小化および正則化に関する読書メモである。

4.3 経験リスク最小化(ERM)

前節では,最尤推定は「負の対数尤度」の最小化であると説明していた。本節では,この対数損失項 l(\boldsymbol{y}_n, \boldsymbol{\theta}; \boldsymbol{x}_n) = - \log p(\boldsymbol{y}_n \mid \boldsymbol{x}_n, \boldsymbol{\theta})を他の損失関数で置き換えることによって,MLEを一般化する。


 \begin{align}
\mathcal{L}(\boldsymbol{\theta}) = \frac{1}{N} \sum_{n=1}^N l(\boldsymbol{y}_n, \boldsymbol{\theta}; \boldsymbol{x}_n) \\ \\
\end{align}

これは,経験分布の下での損失の期待値なので,経験損失最小化(empirical risk minimization, ERM)として知られる。

二値分類問題のためのいろいろな損失

4.3.1 例:誤分類率の最小化

本項では,分類問題のときに使われることがある0-1損失について説明している。

0-1損失では,正しく分類できているのであれば0,分類を誤った場合は1となるような損失である。訓練データ上での経験的な誤分類率になっている。

4.3.2 代理損失

本項では,損失の損失の最小化が難しいときに,代わりに用いられる代理損失(surrogate loss function)について説明している。

0-1損失は,横軸の値が0より大きければ0,0より小さければ1となるような,不連続な変化をする損失関数である。上図に示す通り滑らかではないため,最適化が難しい。

そのかわりに,対数損失 l_{ll}(\tilde{y}, \eta) = \log (1+e^{-\tilde{y} \eta})を用いると,損失は連続な関数となる(図の赤線)。

4.4 その他の推定手法

4.4.1 モーメント法

本項では,パラメータの推定方法であるモーメント法(method of moments, MOM)について説明している。

MLEを解くには,負の対数尤度の最小値を求めるために, \nabla^2_{\boldsymbol{\theta}} \mathrm{NLL} (\boldsymbol{\theta}) = \boldsymbol{0}を解く必要がある(ただし, \mathrm{NLL} (\boldsymbol{\theta})は負の対数尤度)。ただし,この計算は困難な場合がある。

モーメント法は,パラメータを求めるための単純な方法である。理論的な k次モーメントは \mu_k = E[ Y^k ]で,経験的なモーメントは


 \begin{align}
\hat{\mu}_k = \frac{1}{N} \sum_{n=1}^N y_n^k \\ \\
\end{align}
である。モーメント法では,パラメータの数 K個分の式 \mu_k = \hat{\mu}_kを解けばよい。


たとえば正規分布の場合は,


 \begin{align}
\mu_1 &= mu = \bar{y} = \frac{1}{N} \sum_{n=1}^N y_n \\ \\
\mu_2 &= \sigma^2 + \mu^2 = s^2 = \frac{1}{N} \sum_{n=1}^N y_n^2 \\ \\
\end{align}
という連立方程式となる。また正規分布の場合,MOM推定とMLEは一致する。


一方で,MOM推定が妥当でない場合もある。本書P112には,一様分布の場合は,モーメント法によって推定されたパラメータの範囲が,データの実現値とずれるため,MOM推定が妥当ではないという説明がなされている。

4.4.2 オンライン(再帰的)推定

本項では,データセットが逐次的に得られる場合であるときの学習手法であるオンライン学習(online learning)について説明している。なお,データセットが全て手元にある場合の推定をバッチ学習(batch learning)という。


ガウス分布において,サンプル数が t-1のときと tのときにおけるバッチ推定値は,


 \begin{align}
\hat{\mu}_{t-1} = \frac{1}{t-1} \sum_{n=1}^{t-1} y_n, \quad \hat{\mu}_{t} = \frac{1}{t} \sum_{n=1}^{t} y_n \\ \\
\end{align}
となる。

オンライン学習における1ステップの更新則は,


 \begin{align}
\boldsymbol{\theta}_t = f(\hat{\boldsymbol{\theta}}_{t-1}, \boldsymbol{y}_t) \\ \\
\end{align}
,すなわち, tステップ目における推定値は, t-1ステップ目における推定値と tステップ目におけるデータから推定される,というものである。

これを用いると,ガウス分布の平均を再帰的に推定するための更新則は,


 \begin{align}
\hat{\mu}_{t} &= \frac{1}{t} \sum_{n=1}^{t} y_n = \frac{1}{t} ((t-1) \hat{\mu}_{t-1} + y_t ) \\ \\
&=\hat{\mu}_{t-1} + \frac{1}{t} (y_t - \hat{\mu}_{t-1}) \\ \\
\end{align}
となる。この形は移動平均として知られている。


この形を発展させて,


 \begin{align}
\hat{\mu}_{t} = \beta \mu_{t-1} + (1 - \beta) y_t \\ \\
\end{align}
のように,より最近の事例に対して大きな重みを与える(例えば \beta=0.99など)ような移動平均の法制方法を,指数移動平均(exponential moving average)と呼ぶ。

4.5 正則化

MLEやERMでは,訓練データ上の損失を最小化するパラメータを選ぶが,将来のデータに対する損失が小さいとは限らない。この現象を過剰適合(overfitting)とよぶ。
本節では,過剰適合を避けるための手法である正則化(regularization)について説明している。


正則化は,NLL(または経験損失)に罰則項を加えて,以下の形の目的関数を最適化することである。


 \begin{align}
\mathcal{L}(\boldsymbol{\theta}; \lambda) = \left[ \frac{1}{N} \sum_{n=1}^N l(\boldsymbol{y}_n, \boldsymbol{\theta}; \boldsymbol{x}_n)  \right] + \lambda C(\boldsymbol{\theta}) \\ \\
\end{align}

ここで, \lambda \geq 0正則化パラメータ, C(\boldsymbol{\theta})は複雑度罰則である。


特に, \lambda = 1および C(\boldsymbol{\theta})を「負の事前分布の対数」 -\log p(\boldsymbol{\theta})とすると,これは事後分布の最大化と等価になる。


 \begin{align}
\hat{\boldsymbol{\theta}} = \operatorname*{argmax} _{\boldsymbol{\theta}} \log p(\boldsymbol{\theta} \mid \mathcal{D})
= \operatorname*{argmax} _{\boldsymbol{\theta}} [ \log p(\mathcal{D} \mid \boldsymbol{\theta}) + \log p(\boldsymbol{\theta}) ] \\ \\
\end{align}
これは,最大事後確率推定(maximum a posteriori estimation)またはMAP推定と呼ばれる。

4.5.1 例:ベルヌーイ分布のMAP 推定

本項では,ベルヌーイ分布に対するMAP推定について説明している。

硬貨投げの例を考えた時,硬貨の表面(おもてめん)を1回だけしか観測していないと,MLEは \theta_{mle}=1となる。

このような, \theta=0, 1といった極端な値を避けるために,事前分布として p(\theta) = \mathrm{Beta} (\theta \mid a, b)のようにベータ分布を設けると,MAP推定値は,


 \begin{align}
\theta_{map} = \frac{N_1 + a - 1}{N_1 + N_0 + a + b -2 } \\ \\
\end{align}
となる(ただし, N_0, N_1はそれぞれ裏面・表面が出た回数)。MAP推定値から, a, bを調整することで,パラメータの値として極端な値が出ることを避けることができる(ゼロ頻度問題の回避)。

4.5.2 例:多変量ガウス分布のMAP 推定

本項では,多変量ガウス分布のMAP推定について説明している。

多変量ガウス分布の共分散行列 \boldsymbol{\Sigma}のMLEは経験共分散行列であるが,高次元において経験共分散行列は特異になることがある。
そのため, \boldsymbol{\Sigma}の事前分布として逆ウィシャート分布を用いることでMAP推定を行なう。

●ブログ筆者註 :
逆ウィシャート分布(Inverse-Wishart distribution)は,以下のように表される
(出典 : Inverse-Wishart distribution - Wikipedia)。


 \begin{align}
f_{\boldsymbol{\Sigma}}(\boldsymbol{\Sigma}; \boldsymbol{\Psi}, \nu) = 
\frac{ \lvert \boldsymbol{\Psi} \rvert ^{\nu/2} }{2^{\nu p / 2 } \Gamma_p(\nu / 2)} \lvert \boldsymbol{\Sigma} \rvert ^{-(\nu + p + 1)/2} \exp \left( \frac{1}{2} \mathrm{tr} \left( \boldsymbol{\Psi} \boldsymbol{\Sigma}^{-1}   \right)    \right) \\ \\
\end{align}

ただし,

  •  p :  \boldsymbol{\Sigma} \in \mathbb{R}^{p \times p}の次元数
  •  \nu \gt p-1 : 自由度のパラメータ
  •  \boldsymbol{\Psi} \in \mathbb{R}^{p \times p} : パラメータ(正定値行列)
  •  \Gamma_p : 多変量ガンマ関数

である。

これを用いるとMAP推定値は,


 \begin{align}
\hat{\boldsymbol{\Sigma}}_{map} = \lambda \boldsymbol{\Sigma}_0 + (1 - \lambda) \hat{\boldsymbol{\Sigma}}_{mle} \\ \\
\end{align}
となる。なお \boldsymbol{\Sigma}_0 = \mathrm{diag} (\hat{\boldsymbol{\Sigma}}_{mle})とすることが多い。

4.5.3 例:重み減衰

本項では,多項式回帰における過剰適合の抑制する手法について説明している。

多項式回帰では,平均ゼロのガウス事前分布 p(\boldsymbol{w})を設定する。このときMAP推定は,


 \begin{align}
\hat{\boldsymbol{w}}_{map} = \operatorname*{argmin}_{\boldsymbol{w}} \mathrm{NLL}(\boldsymbol{w}) + \lambda \lVert \boldsymbol{w} \rVert _2^2 \\ \\
\end{align}
となる。

この罰則項はL2正則化または重み減衰と呼ばれる。
また線形回帰の場合,L2正則化を用いる方法を特にリッジ回帰という。

4.5.4 検証セットを用いた正則化パラメーターの選択

本項では,正則化パラメータ \lambdaの選び方として,データを訓練集合と検証集合(validation set)を作成し,正則化パラメータ \lambdaの値を変えながら,検証集合への当てはまりの良さで正則化パラメータ \lambdaを選ぶ方法を説明している。

4.5.5 交差検証

4.5.4項の方法は,訓練データのサイズが小さいと,モデルの学習に使えるデータが減るため,モデルパラメータの推定の信頼性が低下するおそれがある。本項ではその対策として,交差検証(cross validation ,CV)を説明している。

4.5.6 早期停止

本項では単純な正則化法として,早期停止(early stopping)について説明している。早期停止は,特に複雑なモデルに対して,実用上はとても効果的であることが多い。

早期停止は反復的な最適化アルゴリズムについて,検証データへの当てはまりが悪くなりはじめた段階で学習を止める。

4.5.7 データの増加

本項では,モデルの複雑度を一定に保ちながらデータ量を増やしたときの挙動について説明している。

モデルの次数を固定した場合における,訓練データサイズと誤差の関係

上図において,問題設定は以下の通りである。

  • 真のモデルは2次多項式に, \sigma^2 = 4の分散のノイズを加える。そのため,真のモデルにおいてもノイズが含まれる(ベイズ誤差(Bayes error)やノイズ下限(noise floor)と呼ぶ)。
  • 推定する際のモデルの次数は,1次,2次,10次,20次として,モデルへの当てはめを行なう。
  • 横軸は訓練データのサンプル数,縦軸は訓練データの誤差・テストデータの誤差とする。

上図からは,以下のようなことが分かる。

  • 1次式を用いた場合(左上図), Nを大きくしてもテスト誤差は小さくならない。すなわち,実際の現象と捉えるにはモデルが単純すぎるといえる(過小適合)。
  • 次数が大きい場合, Nを大きくすると誤差が減る。ただしこの際,次数が小さい方が,推定するべきパラメータ数が少ないので,減少のスピードが速い。
  • 訓練誤差は, Nが比較的小さい場合は, Nとともに増加する。これは,新しいデータを追加するたびに異なる入出力パターンが増えるためである。

まとめと感想

今回は,「4 統計学」における経験リスク最小化および正則化についてまとめた。

代理損失による計算効率化

損失関数として,0-1損失はとても分かりやすいが,最適化計算の多くは損失関数の微分を用いるので,0-1損失は最適化計算には向いていない。
そのため微分が可能な各種損失関数で置き換えるというテクニックが重要だと理解できた。一方で,微分可能な損失においてもそれぞれの性質が異なっているはずなので,どの損失関数を用いるべきかは具体的な問題を通じて理解していきたい。

正則化の理論的裏付け

正則化は,過剰適合を避けるための代表的な手法の1つである。経験的には有効に働くことは分かっているが,ベイズ統計における事後分布最大化の枠組みでも説明ができるということは,正則化が経験的に上手くいく方法というだけではなく,理論的な裏付けにもなっていて興味深かった。


最尤推定は,数理統計学におけるパラメータの推測でもおなじみの手法である。本書では,最尤推定

という2つの方法と同等である,という説明をしていた。

有名な手法についても,別の観点から再解釈すると,発展性が出てきたり理解が深まったりするので,参考になった。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第4章 統計学 ①最尤推定~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第4章 統計学」における最尤推定に関する読書メモである。

4 統計学

本章の名称は「統計学」(Statistics)であるが,主な話題はパラメータの推定である。

4.1 はじめに

本節では,本章のテーマである確率モデルのパラメータをデータから学習することについて説明している。

データ \mathcal{D}から確率モデルのパラメータ \boldsymbol{\theta}を推定するプロセスは,モデル適合(model fitting)または訓練(training)と呼ばれるが,以下の最適化問題に帰着される。


 \begin{align}
\hat{\boldsymbol{\theta}} = \operatorname*{argmin}_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) \\ \\ 
\end{align}
という最適化問題に帰着される。ただし \mathcal{L}(\boldsymbol{\theta})は何らかの損失関数または目的関数である。

4.2 最尤推定(MLE)

本節では,パラメータ推定に対する最も一般的なアプローチである最尤推定(maximum likelihood estimation, MLE)について説明している。

4.2.1 定義

本項では,MLEの定義について説明している。

MLEの定義(MLEにおいて解く最適化問題)は


 \begin{align}
\hat{\boldsymbol{\theta}}_{mle} \equiv \operatorname*{argmax}_{\boldsymbol{\theta}} p(\mathcal{D} \mid \boldsymbol{\theta}) \\ \\ 
\end{align}
である。すなわち尤度をパラメータについて最大化するという問題になる。

ただし通常取り扱うのは,尤度の対数を取った対数尤度であり, l(\boldsymbol{\theta}) \equiv \log p(\mathcal{D} \mid \boldsymbol{\theta})で表される。

また,たいていの最適化アルゴリズムは,コスト関数を最小化するように設計されているため,負の対数尤度を最小化する問題として定義する。

4.2.2 MLE の正当化

本節では,MLEを正当化するいくつかの方法について説明している。すなわち,尤度最大化問題を,別の最適化問題から導くことについて説明している。

ベイズ事後分布の単純な1点近似として捉える

MLEを正当化する方法として,MLEを「一様事前分布の下でベイズ事後分布の一点近似」として捉える方法である。

事後分布最大化によって得られるパラメータ \hat{\boldsymbol{\theta}}_{map}は,事後分布のモードであり,


 \begin{align}
\hat{\boldsymbol{\theta}}_{map} 
&= \operatorname*{argmax}_{\boldsymbol{\theta}} p(\boldsymbol{\theta} \mid \mathcal{D})  \\ \\
&= \operatorname*{argmax}_{\boldsymbol{\theta}} p(\mathcal{D} \mid \boldsymbol{\theta}) + \log p(\boldsymbol{\theta}) \\ \\
\end{align}
で与えられる。

事後分布をデルタ関数 p(\boldsymbol{\theta} \mid \mathcal{D}) = \delta(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{map})とし,事前分布を一様分布とすると, \hat{\boldsymbol{\theta}}_{map} = \hat{\boldsymbol{\theta}}_{mle} となる。

予測分布と経験分布を近づける

MLEを正当化するもう一つの方法は,予測分布 p(\boldsymbol{y} \mid \hat{\boldsymbol{\theta}}_{mle})が以下の経験分布


 \begin{align}
p_{\mathcal{D}}(\boldsymbol{y}) \equiv \frac{1}{N} \sum_{n=1}^N \delta (\boldsymbol{y} - \boldsymbol{y}_n) = \frac{1}{N} \sum_{\boldsymbol{y}} \delta (\boldsymbol{y} - \boldsymbol{y}_n) \\ \\
\end{align}
に対して近くにする,という方法である。


真のパラメータ \boldsymbol{\theta}を持つ確率分布 q(\boldsymbol{\theta}) = p(\boldsymbol{y} \mid \boldsymbol{\theta})と上記の経験分布との距離は,KLダイバージェンス D_{KL} (p \mid \mid q)で測る。
KLダイバージェンスは,「定数+負の対数尤度」という形に変形できるため,KLダイバージェンス最小化は,最尤推定で行なっている「負の対数尤度の最小化」と同じになる。

●ブログ筆者註 :
本書P104ではデルタ関数の演算が出てくるが,ここでは本書P60に登場する「より分けの性質」を用いている。


 \begin{align}
 D_{KL} (p \mid \mid q) = \sum_{\boldsymbol{y}} (p_{\mathcal{D}}(\boldsymbol{y}) \log p_{\mathcal{D}}(\boldsymbol{y}) - p_{\mathcal{D}}(\boldsymbol{y}) \log q(y) ) \\ \\
\end{align}
であるが,右辺第2項は

 \begin{align}
 \sum_{\boldsymbol{y}} \left(   \frac{1}{N} \sum_{n=1}^N \delta (\boldsymbol{y} - \boldsymbol{y}_n) \right) \log p(\boldsymbol{y} \mid \boldsymbol{\theta}) \\ \\
\end{align}
であり,より分けの性質

 \begin{align}
\int_{-\infty}^{\infty} f(y) \delta(x - y) dy = f(x) \\ \\
\end{align}
から, \log p(\boldsymbol{y}_n \mid \boldsymbol{\theta})が出現する。

4.2.3 例:ベルヌーイ分布のMLE

本項では,最尤推定の具体例として,ベルヌーイ分布の最尤推定を説明している。

ベルヌーイ分布の最尤推定量は, Y=0, Y=1となった回数をそれぞれ N_0, N_1とすると, \hat{\theta}_{mle} = N_1/(N_0 + N_1)となり,これは Y=1となった回数なので,直感的な結果と言える。

4.2.4 例:カテゴリカル分布のMLE

本項では,最尤推定の具体例として,カテゴリカル分布の最尤推定を説明している。

カテゴリカル分布は, K面サイコロを振って出る目の確率,のように説明されることが多い。カテゴリカル分布においては,各面が出る確率 \theta_kの総和が1になるという制約があるため,制約付き最適化問題に用いるラグランジュ乗数法を用いる。

最尤推定量は,試行回数の総和を N Y=kとなった回数を N_kとすると, \hat{\theta}_k = N_K / Nとなり,経験的な比率そのものとなる。

4.2.5 例:単変量ガウス分布のMLE

本項では,最尤推定の具体例として,単変量ガウス分布最尤推定を説明している。

単変量ガウス分布なので,推定する対象のパラメータは \boldsymbol{\theta} = (\mu, \sigma^2)すなわち平均と分散である。最尤推定量は,


 \begin{align}
\hat{\mu}_{mle} = \bar{y}, \quad \hat{\sigma}^2_{mle} = \frac{1}{N} \sum_{n=1}^N (y_n - \hat{\mu}_{mle})^2 \\ \\
\end{align}
となる。

分散の推定量として代表的なものは不偏推定量であり,不偏推定量は次のように N-1で割ったものになる。


 \begin{align}
\hat{\sigma}^2_{unb} = \frac{1}{N-1} \sum_{n=1}^N (y_n - \hat{\mu}_{mle})^2 \\ \\
\end{align}

4.2.6 例:多変量ガウス分布のMLE

本項では,最尤推定の具体例として,多変量ガウス分布最尤推定を説明している。

多変量ガウス分布最尤推定の結果は,


 \begin{align}
\hat{\boldsymbol{\mu}} = \bar{\boldsymbol{y}}, \quad \hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (\boldsymbol{y}_n - \bar{\boldsymbol{y}}) (\boldsymbol{y}_n - \bar{\boldsymbol{y}})^T \\ \\
\end{align}
となる。

多変量ガウス分布最尤推定で注意するべきことは, N \lt D,すなわちサンプル数 Nが次元数 Dよりも小さいときに,数値計算が不安定になったり,特異行列になったりするということである。このような時の対処法は,4.5.2項において,正則化を用いる方法が説明されている。

●ブログ筆者註 :
共分散行列のMLEは経験共分散行列になることは,本書P108の式(4.51)に示されている。
経験共分散行列の計算自体は難しい問題ではないが, N \lt Dのときにはランクが落ちるため, \hat{\Sigma}逆行列計算が求められなくなる。
多変量ガウス分布には共分散行列の逆行列が現れるため,例えば推定された共分散行列を用いて予測をする場合などに,正しく予測ができなくなるおそれがある。

4.2.7 例:線形回帰のMLE

本項では,最尤推定の具体例として,線形回帰の最尤推定を説明している。

線形回帰モデルは


 \begin{align}
p(y \mid \boldsymbol{x}; \boldsymbol{\theta}) = \mathcal{N}(y \mid \boldsymbol{w}^T \boldsymbol{x}, \sigma^2) \\ \\
\end{align}
であるが, \sigma^2が固定されているとすると,推定対象は \boldsymbol{w}になる。この最尤推定量は,OLS(Ordinary least suquare)推定量と呼ばれ,

 \begin{align}
\hat{\boldsymbol{w}}_{mle} = (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y} \\ \\
\end{align}
となる。

まとめと感想

今回は,「4 統計学」における最尤推定についてまとめた。

最尤推定の解釈

最尤推定は,数理統計学におけるパラメータの推測でもおなじみの手法である。本書では,最尤推定

という2つの方法と同等である,という説明をしていた。

有名な手法についても,別の観点から再解釈すると,発展性が出てきたり理解が深まったりするので,参考になった。

本記事を最後まで読んでくださり,どうもありがとうございました。

「確率的機械学習 入門編I」を読む ~第3章 多変量の確率モデル③指数型分布族~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第3章 多変量の確率モデル」における指数型分布族・混合モデル・確率的グラフィカルモデルに関する読書メモである。

3.4 指数型分布族

本節では,主要な確率分布を多く含む指数型分布族(exponential family)について説明している。

3.4.1 定義

本節では,指数型分布族の定義や,各種用語を与えている。

指数型分布族は,以下の式で表される分布 p(\boldsymbol{y} \mid \boldsymbol{\eta})である。


 \begin{align}
p(\boldsymbol{y} \mid \boldsymbol{\eta}) 
&\equiv \frac{1}{Z(\boldsymbol{\eta})} h(\boldsymbol{y}) \exp [ \boldsymbol{\eta}^T \mathcal{T}(\boldsymbol{y}) ]  \\ \\
&= h(\boldsymbol{y}) \exp [ \boldsymbol{\eta}^T \mathcal{T}(\boldsymbol{y}) - A(\boldsymbol{\eta}) ]  \\ \\
\end{align}

変数の説明は以下の通りである。

  •  \boldsymbol{\eta} : 指数型分布族のパラメータで,特に自然パラメータ(natural parameter)または正準パラメータ(canonical parameter)と呼ばれる。
  •  h(\boldsymbol{y}) : スケーリング定数,または基底測度
  •  \mathcal{T}(\boldsymbol{y}) : 十分統計量(sufficient statistics)。
  •  Z(\boldsymbol{\eta}) : 正規化定数で分配関数と呼ばれる。
  •  A(\boldsymbol{\eta}) = \log Z(\boldsymbol{\eta}) : 対数分配関数と呼ばれる。

●ブログ筆者註 :
指数分布族の定義は上の通りであるが,いくつかの特徴的がある。

  • 自然パラメータ \boldsymbol{\eta}がデータ \boldsymbol{y}と結合して出てくる部分は,自然パラメータと十分統計量内積で現れる。
  • かつその内積が, \expの中身として現れる。そのため対数尤度においては,内積そのものが現れる。

3.4.2 具体例

本項では,指数型分布族の例として,ベルヌーイ分布を紹介している。

ベルヌーイ分布は以下のように表される。


 \begin{align}
\mathrm{Ber} (y \mid \mu) 
&= \mu^y (1 - \mu)^{1 - y} \\ \\
&= \exp \left[ y \log \left(  \frac{\mu}{1 - \mu}  \right)  + \log(1 - \mu) \right] \\ \\
\end{align}

3.4.3 対数分配関数とキュムラント母関数

本項では,対数分配関数とキュムラント母関数の関係について説明している。

キュムラント母関数(※ブログ筆者補足)

キュムラント母関数のおさらいをする。

確率変数 Xについて,モーメント母関数 M_X(t) = E[ e^{tX} ]で定義される。キュムラント母関数 K_X(t)モーメント母関数の対数として定義される*1


 \begin{align}
K_X(t) = \log(M_X(t)) \\ \\
\end{align}


 n次のキュムラントとは,キュムラント母関数を級数展開したとき


 \begin{align}
K_X(t) = \sum_{n=1}^{\infty} \frac{t^n}{n!} \kappa_n \\ \\
\end{align}
に現れる係数 \kappa_nのことであり,以下を満たす。

 \begin{align}
\kappa_n = \left. \frac{\partial ^n}{ \partial t^n} K_X(t) \right| _{t=0} = \left. \frac{\partial ^n}{ \partial t^n} \log(M_X(t)) \right| _{t=0}  \\ \\
\end{align}


なお本書P91には,1次と2次のキュムラントがそれぞれ平均・分散であるとの説明があるが,これは久保川「現代数理統計学の基礎」のP21に説明がある。

キュムラントと対数分配関数の関係(※ブログ筆者補足)

本書P91には,「指数型分布族の重要な特徴の1つは,対数分配関数の高階導関数が,対応する次数のキュムラントになっている」との説明があった。このことを確認する。


まず分配関数 Z(\boldsymbol{\eta}) = \exp(A(\boldsymbol{\eta}))は正規化定数であるため,


 \begin{align}
\int h(\boldsymbol{y}) \exp [ \boldsymbol{\eta}^T \mathcal{T}(\boldsymbol{y}) ] d \boldsymbol{y} = \exp(A(\boldsymbol{\eta})) \\ \\
\end{align}
となる。この式より, \boldsymbol{\eta} \boldsymbol{\eta}+ \boldsymbol{s}を代入すると,

 \begin{align}
\exp(A(\boldsymbol{\eta} + \boldsymbol{s})) = \int h(\boldsymbol{y}) \exp [ (\boldsymbol{\eta} + \boldsymbol{s})^T \mathcal{T}(\boldsymbol{y}) ] d \boldsymbol{y} \\ \\
\end{align}
となる。


次に十分統計量 \mathcal{T}におけるモーメント母関数・キュムラント母関数を考えると,


 \begin{align}
M_{\mathcal{T}} ( \boldsymbol{s} ) 
&= \int p(\boldsymbol{y} \mid \boldsymbol{\eta}) \exp( \boldsymbol{s}^T \mathcal{T}(\boldsymbol{y})) d \boldsymbol{y} \\ \\
&= \int h(\boldsymbol{y}) \exp [ (\boldsymbol{\eta} + \boldsymbol{s})^T \mathcal{T}(\boldsymbol{y}) - A( \boldsymbol{\eta} ) ] d \boldsymbol{y} \\ \\
&= \exp(- A( \boldsymbol{\eta} )) \int h(\boldsymbol{y}) \exp [ (\boldsymbol{\eta} + \boldsymbol{s})^T \mathcal{T}(\boldsymbol{y}) ] d \boldsymbol{y} \\ \\
&= \exp(A( \boldsymbol{\eta} + \boldsymbol{s}) - A( \boldsymbol{\eta} ) )\\ \\

\therefore \: K_{\mathcal{T}} ( \boldsymbol{s} ) &= \log (M_{\mathcal{T}} ( \boldsymbol{s} ) ) = A( \boldsymbol{\eta} + \boldsymbol{s}) - A( \boldsymbol{\eta} ) \tag{1} \\ \\

\end{align}
となる(なお,確率変数 \boldsymbol{y}のモーメント母関数・キュムラント母関数ではないことに注意する)。


 n次のキュムラントは,


 \begin{align}
\left. \frac{\partial ^n}{ \partial \boldsymbol{s}^n }  K_{\mathcal{T}} ( \boldsymbol{s} ) \right| _{\boldsymbol{s} = 0} \\ \\
\end{align}
で表される。

(1)式の右辺を \boldsymbol{s}微分することを考える。形式的に \boldsymbol{t} = \boldsymbol{\eta} + \boldsymbol{s}と置くと,


 \begin{align}
\frac{\partial}{ \partial \boldsymbol{s} } \left(  A( \boldsymbol{\eta} + \boldsymbol{s}) - A( \boldsymbol{\eta} )  \right)
&= \frac{\partial \boldsymbol{t} }{ \partial \boldsymbol{s} }  \frac{\partial }{ \partial \boldsymbol{t} } A( \boldsymbol{t} ) 
= \frac{\partial }{ \partial \boldsymbol{t} } A( \boldsymbol{t} ) \\ \\
\end{align}
となる。

 \boldsymbol{s} = 0とすると, \boldsymbol{t} = \boldsymbol{\eta}より,


 \begin{align}
\left. \frac{\partial ^n}{ \partial \boldsymbol{s}^n }  K_{\mathcal{T}} ( \boldsymbol{s} ) \right| _{\boldsymbol{s} = 0} 
= \frac{\partial ^n }{ \partial \boldsymbol{\eta}^n } A( \boldsymbol{\eta} ) \\ \\
\end{align}
となり,対数分配関数 A( \boldsymbol{\eta} )の高階導関数が,対応する次数のキュムラントになっている。

このキュムラント母関数は,十分統計量に関するキュムラント母関数であるため,本書P91の(3.84)式・(3.85)式


 \begin{align}
\nabla A(\boldsymbol{\eta}) &= E[ \mathcal{T}(\boldsymbol{y}) ] \\ \\
\nabla^2 A(\boldsymbol{\eta}) &= \mathrm{Cov} [ \mathcal{T}(\boldsymbol{y}) ] \\ \\
\end{align}
のように,このキュムラント母関数の1次・2次導関数は,十分統計量の1次・2次モーメントになる。

3.4.4 最大エントロピー基準による指数型分布族の導出

本項では,関数 f_k(\boldsymbol{x})の期待値 F_K が与えられている状況で,確率分布 p(\boldsymbol{x})を推定する,ということを考える。

推定するには,何らかの基準が必要になるが,この基準として,

とすると,確率分布として指数型分布族が得られる。

たとえば f_1(x)=x, f_2(x)=x^2としてエントロピー最大化を行なうと,ガウス分布が得られる。

3.5 混合モデル

本項では,より複雑な確率分布の作り方として,単純な分布の凸結合をとる混合モデル(mixture model)について説明している。

混合モデルは,


 \begin{align}
p(\boldsymbol{y} \mid \boldsymbol{\theta}) &= \sum_{k=1}^K \pi_k p_k(\boldsymbol{y}) \\ \\
\sum_{k=1}^K \pi_k &= 1, \quad \pi_k \in [0, 1] \\ \\
\end{align}
で表される。

3.5.1 混合ガウスモデル

本項では,ガウス分布の混合モデルである混合ガウスモデル(Gaussian mixture model, GMM)


 \begin{align}
p(\boldsymbol{y} \mid \boldsymbol{\theta}) &= \sum_{k=1}^K \pi_k \mathcal{N}(\boldsymbol{y} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k ) \\ \\
\end{align}
について説明している。

3つの2次元ガウス分布の混合モデル

混合ガウスモデルは,クラスタリングに用いることができるが,各正規分布への所属が確率的に決められるので,ソフトクラスタリングと呼ばれる。

3.5.2 ベルヌーイ混合モデル

本節では,データが2値の場合における,混合ベルヌーイモデル(Bernoulli mixture model, BMM)を用いる。


 \begin{align}
p(\boldsymbol{y} \mid z=k, \boldsymbol{\theta}) = \prod_{d=1}^D \mathrm{Ber}(y_d \mid \mu_{dk}) 
=  \prod_{d=1}^D \mu_{dk}^{y_d} (1 - \mu_{dk})^{1 - y_d} \\ \\
\end{align}

3.6 確率的グラフィカルモデル

本節では,

  • 一部の変数がほかの変数に対して条件付き独立であることを仮定する
  • この条件付き独立の過程をグラフによって表現する

という確率的グラフィカルモデル(probabilistic graphical model, PGM)について説明している。

3.6.1 表現

確率的グラフィカルモデルは,グラフ構造を用いて条件付き独立性を表現する。グラフが有向非巡回グラフ(directed acyclic graph, DAG)のとき,このモデルはベイジアンネットワークと呼ばれる。

代表的な確率的グラフィカルモデルとして,マルコフモデルが挙げられる。

まとめと感想

今回は,「3 多変量の確率モデル」における指数型分布族・混合モデル・確率的グラフィカルモデルについてまとめた。

指数型分布族の各種性質

指数型分布族は,数理統計学でもおなじみの確率分布である。

本節では,

  • 対数分配関数とキュムラント母関数の関係
  • 最大エントロピー基準による指数型分布族の導出

など,様々な性質が紹介されていた。

本節では,また各種性質の導入にとどまっていたが,パラメータ推定をはじめ機械学習の様々なシーンにおいてその特性を発揮するものと考えられる。
例えば第12章において,一般化線形モデルの説明で指数型分布族が活用できるとのことなので,その際の説明を楽しみにしたい。

確率分布の複雑化

実データや実システムを扱おうとすると,単純な多変量確率モデルでは扱いにくいことが想定される。複雑な確率分布を表現する方法として混合モデルと確率的グラフィカルモデルが紹介されていた。

これらのモデルだけでも,かなりの応用が利きそうではあるが,混合モデルは確率分布の凸結合であるし,確率的グラフィカルモデルは確率分布の積で表現されるなど,まだまだ非線形化の余地はあると考えられる。そのような確率分布の代表例が生成モデルであると考えられるが,混合モデルや確率的グラフィカルモデルとの特徴の違いについて今後理解していきたい。


本記事を最後まで読んでくださり,どうもありがとうございました。

*1:この定義は,Wikipediaの定義(キュムラント母関数 - Wikipedia)から引用している。一方で,久保川「現代数理統計学の基礎」では,「キュムラント母関数は特性関数の対数」として定義している。ただし本書においては複素数が出てこないようなので,本記事では特性関数ではなくモーメント母関数の対数である,として説明を進めている。

「確率的機械学習 入門編I」を読む ~第3章 多変量の確率モデル②線形ガウスシステム~

はじめに

持橋大地・鈴木大慈 監訳「確率的機械学習 入門編I」は,世界的に評価の高いK.P.Murphy著 "Probabilistic Machine Learning (Book1)" の和訳であり,確率モデルに基づく機械学習,深層学習といった基礎が丁寧に整理されている。私は統計検定1級として数理統計の基礎は学んできたが,機械学習・深層学習は実務に応じて場当たり的に学んできた。実務での応用に向けて,機械学習・深層学習の基礎を体系的に学び,チームの技術力を底上げしたいと考えている。また読書メモに自身の理解をまとめることで,製造業に携わる若いエンジニアにとっても有益な知識を還元できればと考えている。

※なおボリュームが多い本なので,知っているところは端折りながら読み進めたい。


本記事は,「第3章 多変量の確率モデル」における,線形ガウスシステムに関する読書メモである。

3.3 線形ガウスシステム

本節では,

  • 入力 \boldsymbol{z} \in \mathbb{R}^Lが未観測のベクトル(潜在変数)
  • 出力 \boldsymbol{y} \in \mathbb{R}^Dノイズありの観測されたベクトル

であり,かつ \boldsymbol{z}および, \boldsymbol{z}に条件付けられた \boldsymbol{y}の確率分布がそれぞれ


 \begin{align}
p(\boldsymbol{z}) &= \mathcal{N}(\boldsymbol{z} \mid \boldsymbol{\mu}_z, \boldsymbol{\Sigma}_z) \\ \\
p(\boldsymbol{y} \mid \boldsymbol{z} ) &= \mathcal{N}(\boldsymbol{y} \mid \boldsymbol{W} \boldsymbol{z} + \boldsymbol{b}, \boldsymbol{\Sigma}_y) \\ \\
\end{align}
であるとする。これを線形ガウスシステム(linear Gaussian System)と呼ぶ。

3.3.1 ガウス分布に対するベイズの定理

本節では,線形ガウスシステムにおける p(\boldsymbol{z})および p(\boldsymbol{y} \mid \boldsymbol{z} )に対して,ベイズの定理を適用することで,潜在変数の事後分布 p(\boldsymbol{z} \mid \boldsymbol{y} )の式を与えている。

ベイズ則適用の流れ(※ブログ筆者補足)

ガウス分布においてベイズ則を適用するモチベーションは,観測できない潜在変数を,観測変数を用いて推定するということである。
ベイズ則適用の流れは以下の通りである。

  • Step 1. 潜在変数と,潜在変数が与えられた場合の観測変数をモデル化する。
    • 潜在変数の確率分布は, p(\boldsymbol{z})である。これには観測されたデータ \boldsymbol{y} の情報がないので,データがない状態における事前知識(先験情報)からモデル化している。
    • 観測変数の条件付き確率分布は, p(\boldsymbol{y} \mid \boldsymbol{z} )である。これは, \boldsymbol{z}が与えられたときに \boldsymbol{y} が発生するモデルを示している。もし線形の関係を仮定するのであれば線形ガウスシステムになる。
  • Step 2.  \boldsymbol{y}, \boldsymbol{z}の同時分布を求める。
    • 同時分布の式 p(\boldsymbol{y} \mid \boldsymbol{z}) p(\boldsymbol{z}) = p(\boldsymbol{y}, \boldsymbol{z})から,同時分布を求める。
  • Step 3.  p(\boldsymbol{z})の事後分布 p(\boldsymbol{z} \mid \boldsymbol{y})を求める。
    •  p(\boldsymbol{y}, \boldsymbol{z}) = p(\boldsymbol{y}) p(\boldsymbol{z} \mid \boldsymbol{y})のように,同時分布から条件付き分布を求める際に,変数を入れ替える。
    • ベイズ則は p(\boldsymbol{z} \mid \boldsymbol{y}) p(\boldsymbol{y}) = p(\boldsymbol{y} \mid \boldsymbol{z}) p(\boldsymbol{z}) に表されるが,計算の際にはいったん同時分布が必要になる。


同時分布は,本書P84に与えられているが,対応する変数は ( p(\boldsymbol{z})^T,  p(\boldsymbol{y})^T)^Tのように,  p(\boldsymbol{z})が上段に来るようになっている。

同時分布の共分散行列は,「3.3.2 導出」に与えられているように,2次の項を整理することで,共分散行列の逆行列(精度行列)を求める。
同時分布の平均ベクトルは, \boldsymbol{z}については p(\boldsymbol{z})の平均 \boldsymbol{\mu}_zをそのまま用いればよい。 \boldsymbol{y}については, \boldsymbol{z} = \boldsymbol{\mu}_zが与えられたときの \boldsymbol{y} = \boldsymbol{W} \boldsymbol{z} + \boldsymbol{b}を用いればよい。


最後の潜在変数の事後分布は,同時分布と条件付ける変数 \boldsymbol{y}を用いて,多変量正規分布の条件付き分布の思い出し方(増補改訂版) #統計検定 - jiku log に示した公式を用いれば導出できる。

  

3.3.3 例:未知のスカラー値の推定

本節では,線形ガウスシステムにおけるベイズの定理の適用例として,潜在変数・観測変数がともに1次元の場合を例に挙げている。

1次元ガウスシステムでは,

  • 事後分布の精度(分散の逆数)は,事前分布の精度と尤度の精度の和になる。
  • 事後分布の平均は,事前分布の平均と最尤推定量の凸結合になる。

という特徴がある。

また,尤度を固定した場合における,事前分布の分散の大きさと事後分布の形状を下図に示す。

1次元ガウス分布における事前分布・尤度・事後分布

左図のように,事前分布の分散が小さい場合,事後分布は事前分布と尤度の間に来るようになる。
一方で右図のように,事前分布の分散が大きい場合は,事後分布は尤度にかなり近いものになる。

3.3.4 例:未知のベクトルの推定

本節では,変数がベクトル量であるときの例を挙げている。

2次元ベクトルにおいて,観測データと事前分布,事後分布の例を下図に示す。

2次元カウス分布にしたがうベクトルzに対するベイズ推定

データのばらつきは,縦軸方向よりも横軸方向の方がばらついている。そのため事後分布も,横軸方向の分散が大きいものになっている。

まとめと感想

今回は,「3 多変量の確率モデル」の線形ガウスシステムについてまとめた。

線形ガウスシステムは多変量確率モデルの代表選手

線形ガウスシステムは,

  • 「潜在変数と観測変数というモデルから潜在変数を推定する」,というよくある問題設定である。
  • 潜在変数と観測変数は,最も基本的な依存関係である線形モデルである。

という特徴を持つことから多変量確率モデルの代表選手であると言える。

またガウス分布を扱っているので,確率分布を求めることは,平均と共分散行列を求めることに帰着されるので,目指すべき計算もわかりやすい。

ただし計算自体は,かなり手間がかかるものである。その理由は,多変量ガウス分布 \expの中身に,ブロック行列の逆行列が出てくるためである。本書ではあまりこれらの計算について触れられていないが,行列計算については別書を参考にしながら読み進めるのが良いと感じた。


本記事を最後まで読んでくださり,どうもありがとうございました。