Uncertainty in Gradient Boosting via Ensembles[A4一枚まで備忘録]

A4用紙1枚にまとめるイメージでメモを残そうという取り組み。

今回は機械学習でとりわけ、回帰の問題を扱うことが仕事であった際に参照できるメモの一つとして残したい。
先日、Catboostのドキュメントに参考文献としてUncertainty in Gradient Boosting via Ensemblesの存在を知った。業務において、価格の予測とか、期間の予測とかマーケティングのデータでも良く扱うものなので、分類だけでなく回帰についてもバランス良く知見を貯めたいところ。

目次

Abstract
Introduction
Background and related work
Generating ensembles of GDBT models
Experiments on Synthetic Data
Experiments on classfication and regression datasets
Conclusion
Broader Impact

Abstract

機械学習のモデルは予測において点推定しかしないが、実務上は不確実性を定量化したいニーズがある。この研究では、予測の不確実性推定値を導出するためのアンサンブルベースの枠組みが紹介されている。予測における不確実性は「データに起因する複雑さとノイズの不確実性」「特徴空間において、情報がないことに起因する不確実性」に分解される。
手法としては確率的勾配ブースティングと確率的勾配ランジュバンブースティング(SGLB)が扱われている。SGLBは一つの勾配ブースティングモデルからのアンサンブルを可能にし、モデルの複雑さを減らすことができる。また、予測性能の向上について様々なデータセットを使って示されている。

Introduction

一般的な勾配ブースティングについて、回帰の際に点推定しか返さないことを挙げている。NGBoostについても言及している。NGBoostなら、勾配ブースト決定木(GBDT)を訓練して平均値を返せるとのこと。ただし、データの不確実性に起因する不確実性しか捕捉できないという問題がある。モデルに関する知識の欠落によって生じる不確実性に関して扱うことはできない。近年ではニューラルネットワークを用いて知識を起因とする不確実性の研究がなされているらしい。今回扱う、アンサンブルアプローチは全体の不確実性を「データの不確実性」と「知識の不確実性」に分解する。扱うモデルは確率的勾配ブースティングと、確率的勾配ランジュバンブースティング(SGLB)となっており、SGLBの利点として、漸近的に真のベイジアン事後分布からサンプリングされることが挙げられている。また、SGLBによって1つの勾配ブースティングのみをモデルを使ってアンサンブルをすることができるとされている。これにより計算の複雑さを減らすことができる。結果として、予測性能が高まり、誤りやドメインの外側のデータの入力の検知ができるようになったとされている。

Background and related work

  • ベイジアンのアプローチ
    ベイズ統計学を使っている。予測の多様性について考察することで、知識の不確実性の推定を行うことができるとしている。さらに、確率モデルのアンサンブルを行なっている。アンサンブルモデルの不一致さ、広がりの程度が知識の不確実性につながるとしている。事後分布の近似のための研究が色々となされ、ニューラルネットワークのモデルなども使われた。モデルの事後分布をもとに、予測の事後分布を手に入れる。予測の事後分布のエントロピーが、予測の不確実性を推定する。知識の不確実性は、全ての不確実性からデータの不確実性の期待値を差し引いた相互情報量で表される。確率的な回帰モデルのアンサンブルを考える時、予測事後分布のエントロピーと相互情報量を推定することができない。代わりに総変動の法則に従ってなら推定することができるが、1~2次のモーメントに基づいているだけなので、高次元での不確実性に関する見落としが起きうる。
  • 勾配ブースティング
    弱学習器は浅い決定木によって形成される。ただ、勾配ブースティングの回帰モデルは点での予測しかしない。予測の不確実性に関する研究はほとんどされていない。近年、NGBoost (Natural Gradient Boosting)が登場し、予測の不確かさを推定することができるようになった。NGBoostは入力変数xが与えられたもとでの、出力変数の条件付き事後分布の推定を行ってくれる。

Generating ensembles of GDBT models

NGBoostはデータの不確実性を捉えることができるものの、知識の不確実性については捉えることができない。
知識の不確実性を捉えるための方法としては、事後分布からサンプリングされた予測モデル(GBDT)のアンサンブルが挙げられる。
そもそもGBDT自身が木のアンサンブルであるが、そのGBDTのアンサンブルを考える。

アンサンブルをする方法としては、確率的勾配ブースティングのアプローチで生成された独立なモデルを扱うことがあげられている。
もとのモデルにノイズを注入して、出力をぶれさせることで予測の分布を得ようとしている。ただ、これは事後分布の推定を保証しないらしい。

ただし、確率的勾配ランジュバンブースティング(SGLB)アルゴリズムを用いることで真の事後分布からGBDTモデルをサンプリングすることができるらしい。
非凸の損失関数においての最適化の収束のために、確率的勾配ランジュバン力学が用いられているとのこと。ランジュバン力学ってなんだろうか。

確率的勾配ブースティングとの違いは2つある。

  • 1.ガウシアンノイズが勾配に明示的に加えられている。
    • 正規分布の平均は0で分散に逆拡散温度というものがある。ノイズをコントロールするものであることに間違いはない。確かに物理っぽい。STANのリープフロッグとかを思い出した。このノイズが最適値を見つけるのに役立つらしい。
  • 2.ブースティングモデルの更新式に正則化パラメータを設けている。
    • 定常分布に弱く収束していくマルコフ連鎖にパラメータが従うらしい。

対数尤度に正則化項が残るものの、比例するらしく、真の事後分布に従うと見なせるらしい。

確率的勾配ブースティングや確率的勾配ランジュバンブースティングは計算する上で重大なオーバーヘッドが生じるらしくイタレーションを減らすなど質を下げたり複雑さを増したりする必要があるとのこと。それを避けるために、仮想SGLBアンサンブルというものが提案されている。

SGLBアルゴリズムのイタレーションごとのパラメータは弱く定常に収束するマルコフ連鎖に従うことから、それらをモデルのアンサンブルと見なす。
しかしながらモデル同士が強く相関してしまうので、アンサンブルの質が下がる。ただし、各々のパラメータのK番目の集合のみを使うことで克服できるとある。
K番目のパラメータ集合のみのモデルで構築したアンサンブルを仮想アンサンブルと呼んでいる。

Experiments on Synthetic Data

ここではGBDTモデルのアンサンブルアルゴリズムを合成データを用いて実験している。
合成データとしては、以下の2つの1次元データセットが扱われている。

  • 合成データセット
    • 分散が増加していく絶対値x
    • 分散が減少していく絶対値x
  • 比較するモデル
    • SGLB
    • SGB
    • vSGLB(仮想SGLB):SGLBから生成
  • 結果
    • 分散が減少していく絶対値xのケースだと、データの不確実性に関して、SGBではサンプルの領域外で過剰に不確実性が増してしまう。
    • 分散が増加していく絶対値xのケースだと、知識の不確実性に関して、SGBではサンプルの領域外で不確実性が増してしまう。

サンプルの領域外のほうが、領域内よりも不確実性が大きいはずだが、実験データではそうではないらしい。

Experiments on classfication and regression datasets

GBDTモデルのアンサンブルアルゴリズムを色々なデータのいろいろなタスク(回帰、分類)で予測性能をみている。

  • タスク
    • 回帰タスク
    • 分類タスク
  • データセット
    • 18種類くらいある。(Wineとか、解約のデータとかAmazonのデータとかいろいろ)
  • 比較するモデル
    • SGB:このアンサンブルは10個の独立したモデルで、それぞれ1K本の木を持つ。
    • SGLB:このアンサンブルは10個の独立したモデルで、それぞれ1K本の木を持つ。
    • vSGLB(仮想SGLB):[500, 1000]の区間から50番目のモデルをアンサンブルに追加していっている。
  • どう実装しているか
    • CatBoostライブラリをベース
    • 分類では二値分類のクラスの確率分布を生成
    • 回帰では正規分布の平均と分散を生成
    • 負の対数尤度を最適化している
    • ハイパーパラメータはランダムサーチ
  • 何の指標を見ているか
    • 負の対数尤度
    • 誤分類率
    • RMSE
    • 不確実性の合計、知識の不確実性:Prediction-Rejection Ratio(PPR)というもので測る。どれくらい不確実性が誤差と順番と相関するかを測るらしい。100が最大で、ランダムだと0になる。
    • AUC:領域外のデータに関する評価に用いる。領域外データは合成して生成した。
    • 分類:負の対数尤度(NLL)
      • アンサンブル手法に軍配が上がるが、違いは軽微である。
      • 仮想アンサンブルが勝るときがある。
    • 分類:誤分類率
      • アンサンブル手法に軍配が上がるが、違いは軽微である。
    • 回帰:RMSE
      • アンサンブル手法に軍配が上がるが、違いは軽微である。
    • SGBやSGLBはパフォーマンスが良い。
      • vSGLBはわずかにパフォーマンスが悪い。自己相関があるなかでアンサンブルしているためパフォーマンス悪化が生じていると考えられる。
        • ただ単一モデルのものよりかは良い。
        • コストをかけずに精度を上げたいという観点ではvSGLBはよい選択肢となる。

    Conclusion

    この論文では、GBDTを用いたアンサンブルベースの不確実性推定モデルが扱われた。SGB、SGLB、仮想SGLBについて合成データや、分類や回帰のタスクにおいて比較した結果、アンサンブル手法が精度面で秀でているわけではないことが示された。
    しかし、知識の不確実性(モデルで十分に想定できていない不確実性)を扱う際にアンサンブル手法が役立つことがわかった。
    知識の不確実性を扱う場合に、計算コストの少ない手法として仮想のSGLBが提案され、精度において十分ではないが、低コストで知識の不確実性を扱う観点において有用であることが示された。

    Broader Impact

    • 不確実性の推定は、機械学習モデルをより安全に、より信頼できるものにするという観点で重要。
    • GBDT系の手法はニューラルネットワーク系の手法よりも計算コストが低くて良い。
    • 不確実性の推定が貧弱なものとならないように引き続き研究が求められる。

    追記

    不確実性について扱われた、ソースコードの載っている記事を集める。

    [Python]音声解析でトランペットの録音データを可視化してみる

    はじめに

    私はデータ分析以外の趣味としては、トランペットの演奏があります。個人でプイプイ吹いているだけのしょうもない趣味ですけどね。
    先日、楽器屋に新しいマウスピース(楽器を口に当てる金属部分、振動させる部分)を買いに行ったのですが、吹いているうちによくわからなくなったので、録音した音声波形などを可視化することで、いいなと思う音がどんな波形をしているのかを明らかにしたいと思います。

    進め方

    • 録音データを用意(iPhoneのムービー)
    • FFmpegで動画から音声のみを切り出し、サンプルレートを変え、モノラルにする
    • Pythonで基本周波数(F0)などを可視化

    データの準備

    今回はチューニングB♭(ベー)の音3種類(Yamaha 11C4、Bach 5C、Bach 3C)と、練習曲(ルパン三世のフレーズ)の音2種類(Yamaha 11C4、Bach 5C)の計5つの音源を扱います。

    マウスピースのイメージ

    • Yamaha 11C4:オールマイティー(これまで使っていたやつ)
    • Bach 5C:少し深くて、少し内径が小さい
    • Bach 3C:少し深くて、少し内径が小さい(5Cと3Cの違いは当て心地?)

    指標の説明

    • 音声波形:波動の振動の大きさを表す非負のスカラー量を縦軸に、横軸を時間とした場合のグラフを音声波形と呼ぶ。
    • ピッチマーク:単位波形を抽出する基準点でエポックマーカーと呼ばれる。音が鳴った際に1をとり音がない際は0を取る。
    • 基本周波数:信号を正弦波の合成で表したときの最も低い周波数成分の周波数で、音程に該当する。
    • 相関係数:ピッチマークの検出に使う指標でNCCF(Normalized Cross-Correlation Function)というものを使って計算している。自己相関のようなもので-1から1の間を取る。

    音声ファイルのあるディレクトリで以下のコマンドを実行します。(私の環境はMacです。)

    音声波形などの可視化〜チューニングB♭(ベー)〜

    ここでは単純にチューニングベーの音を比較します。まず音源ですが、以下で実際に聴けるようになっています。(耳を痛めないように、最も小さい音で聴いてから調整してください。)
    しかし、慣れないマウスピースで音程が不安定になっているのが悲しいですね。

    ・Yamaha 11C4のチューニングベー

    ・Bach 5Cのチューニングベー

    ・Bach 3Cのチューニングベー

    以前書いたブログと同様に、pyreaperというREAPERのラッパーとなっているライブラリを用います。REAPERはGoogleが開発していて、C++で書かれています。

    音声波形の比較


    やはり慣れている分、11c4が安定してそうな波形に見えますね。3Cは結構ブレてそうです。実際の音源でもわかりますね。

    ピッチの比較


    発音のあったところで1となる指標です。今回は伸ばしているだけなので、何の面白味もないです。

    基本周波数(F0)の比較


    音程をつかさどるF0に関しては3Cが安定してそうです。聴いたらブレブレなのに。

    相関係数の比較


    いつも使っている11c4がブレにくいようで、5Cが比較的に上下動が多いように見えます。

    音声波形などの可視化〜練習曲〜

    結局、マウスピースは5Cを購入して、その後スタジオで練習曲の音を比較しました。先ほどと同様に音源が以下で実際に聴けるようになっています。(耳を痛めないように、最も小さい音で聴いてから調整してください。)私は、実際に録音して、マウスピースでこんなに変わるんだと驚きを隠せないです。チューニングベーを比較するだけではわかりにくいんですかね。私としては、5Cの音の方が好きです。

    ・Yamaha 11C4での練習曲

    ・Bach 5Cでの練習曲

    音声波形の比較


    5Cのほうが上下の大きさが対称に近いような印象を受けます。いい音を捉える参考になるんですかね?

    ピッチの比較


    おそらく、11c4の音源の息継ぎのタイミングがやや長くなった結果、ピッチに差が出たのではないかと思います。

    基本周波数(F0)の比較


    比較が難しいですが、最後の伸ばしは5Cのほうが安定しているようです。

    相関係数の比較


    5Cのほうが上下動をしているように見えます。いいなと思う音の方が上下動しやすいのでしょうか。気になるところ。

    追記:スペクトログラムの比較

    どうやら音色を考えるうえで、縦軸に周波数、横軸に時間をとったスペクトログラムというものを可視化するのも良いらしいです。参考文献に従って行ってみます。

    5Cよりも11C4が赤い線が1つ多いように見えます。幅広く音が鳴っていてモサっとしてしまっているのでしょうか。私としては5Cの方が音が鋭くて好きです。

    おわりに

    音声波形で音色の比較がうまくできるかどうか試してみたのですが、良いと思われる音とそうでない音で相関係数や音声波形などのデータの可視化において違いが生じやすいのかなという印象でした。そして音色をつかさどるとされるスペクトログラムにおいても、違いがありそうでした。
    もっと音源を集めていっていい音の研究をしていこうと思います。あと、音楽の理論の方も細々と学んでいきます。

    参考文献

    IPythonデータサイエンスクックブック ―対話型コンピューティングと可視化のためのレシピ集
    pyreaperで音声データのピッチを掴むためのメモ書き(F0の抽出)
    Pythonで音声解析 – 音声データの周波数特性を調べる方法
    基本周波数についてのまとめ
    イラストで学ぶ 音声認識 (KS情報科学専門書)

    [数理統計学]統計的検定のまとめ


    通勤電車のなかで私が勉強する用のシリーズ第5弾です。今回は統計的検定についてまとめておこうと思います。

    【これまでのシリーズへのリンク】



    目次

      統計的仮説検定
      検出力
      一様最強力検定
      ネイマン-ピアソンの基本定理
      不偏検定
      尤度比検定
      尤度比検定(1標本のケース)
      尤度比検定(2標本のケース)
      \( \chi^2 \)適合度検定
      分割表の検定

      統計的仮説検定

      確率ベクトル(標本確率変数)である\( X = (X_1, \dots, X_n) \)が標本空間\( \mathfrak X \)上にある分布\( P_{\theta} , \theta \in \Theta \)に従うとして、空集合ではない\( \Theta_0 \)は\( \Theta \)の真部分集合であらかじめ決まっているものとする。その前提のもとで、Xの従う分布が
      帰無仮説
      $$ P_{\theta} , \theta \in \Theta_0 $$
      なのか、対立仮説
      $$ P_{\theta} , \theta \in \Theta $$
      なのかを主張する行為のことを統計的仮説検定と呼ぶ。なお、\( \Theta_0 \)が一つだけの場合は単純仮説、複数ある場合は複合仮説と呼ぶ。

      仮説検定を行う際は、先ほどあげた標本空間\( \mathfrak X \)を、\( \mathfrak X = C \cup C^C \)との二つの領域に分割する。ここで、\( C \)は棄却域を\( C^C \)は受容域を表す。その際、以下のような棄却域の定義関数を考える。

      $$
      \varphi (X) =
      \begin{cases}
      1 \ (x \in C ) \\
      0 \ (x \in C^C )
      \end{cases} \\
      0 \leq \varphi (x) \leq 1
      $$

      この定義関数\( \varphi (X) \)は検定関数と呼び、標本xを得たときに、確率\( \varphi (X) \)で帰無仮説(\( H_0 \))を棄却する際に用いる。先ほどの定義関数の場合、0-1の値しかとらないので、非確率化検定と呼ばれる。他方、0-1でない場合は確率化検定と呼ばれる。

      棄却域を決めると自ずと生じるものが、第一種の誤りと第二種の誤りです。

      第一種の誤り
      帰無仮説(\( H_0 \))が正しいにも関わらず、それを棄却してしまうことによる誤り。

      第二種の誤り
      対立仮説(\( H_1 \))が正しいにも関わらず、帰無仮説(\( H_0 \))を受容したことによる誤り。

      このどちらをも共に小さくすることはできません。
      図にすることで理解が捗ります。

      ここでは2種類の棄却域\( C_1, C2 \)を考えます。
      図からわかるように、
      $$ P(C_2 | H_0) < P(C_1 | H_0) \\
      P(C_2^C | H_1) < P(C_1^C | H_1)
      $$
      となり、\( C_2 \)は第一種の誤りが小さい代わりに、第二種の誤りが大きくなっているのがわかります。

      検出力

      対立仮説\( H_1 \)が正しい時に、\( H_1 \)を受容する確率を検出力と呼ぶ。検定関数\( \varphi (X) \)の良さを表す。

      $$\beta (\theta; \varphi ) = E_{\theta} \left [ \varphi (X) \right ] \ ( \theta \in \Theta_1 ) $$

      以下の\( \alpha \)のことを検定関数\( \varphi (X) \)の大きさと呼ぶ。
      $$ \alpha_0 = \sup_{\theta \in \Theta_0} E_{\theta}(\varphi(X))
      $$
      第一種の誤りをする確率に関しての、帰無仮説のもとに成り立つ\( \Theta_0 \)空間上での上限となる。
      ここで\( \alpha \)を
      $$0 \leq \alpha \leq 1$$
      としたとき、
      $$ \sup_{\theta \in \Theta_0} E_{\theta}(\varphi(X)) \leq \alpha
      $$
      となる検定を有意水準\( \alpha \)の検定と呼ぶ。
      ようは、多くの人が有意差あるとかないとか言っているのは、この第一種の誤りをする確率のことで、社会科学では5%などに置くことが多い。

      一様最強力検定

      \( \alpha : 0 \leq \alpha \leq 1 \)が与えられたとき、水準\( \alpha \)の検定の中で検出力(\( H_1 \)が正しいときに、\( H_1 \)を受容できる確率)を一様に最大にするものが存在する場合、そのような検定を一様最強力検定(Uniformly Most Powerful test:UMP検定)と呼ぶ。

      $$ E_{\theta} \left [ \varphi (X) \right ] \leq E_{\theta} \left [
      \varphi_0 (X) \right ], \ \forall \theta \in \Theta_1 $$

      \( \theta \)の値が大きいときは、それが対立仮説に従うと採択する確率が高まり、帰無仮説の領域だと、値が小さいにも関わらず、対立仮説に従うと言ってしまっているので、検出力が下がる。

      検出力関数の例

      $$ \pi_{\gamma} (\theta) = P \left [ H_0 \ is \ rejected | X_1, \dots, X_n \sim f(x | \theta) \right ] \\
      \pi_{\gamma} (\theta) \in [0, 1]
      $$

      という検出力関数について、Xの従う分布を以下のように仮定する。
      $$ X_1, \dots, X_n \sim N(\theta, 25) $$
      ここで、帰無仮説で
      $$ H_0 : \theta \leq 17 $$
      とする際の検出力について考える。

      標本平均の\(\bar X\)が
      $$ \bar X > 17 + 10 / \sqrt n $$
      に従うとすると、検出力関数は、
      $$ \pi_{\gamma} = P \left [ \bar X > 17 + 10/ \sqrt n | \theta \in \mathbb R \right ] \\
      = P \left [ \frac{\sqrt n (\bar X – \theta)}{5} > \frac{17+ 10/ \sqrt n – \theta}{5 / \sqrt n} \right ] \\
      = 1 – \Phi \left ( \frac{17+ 10/ \sqrt n – \theta}{5 / \sqrt n} \right ) $$
      となる。

      一応、可視化のためのコードも載せておきます。

      ネイマン-ピアソンの基本定理

      確率分布\(P_{\theta}\)の確率密度関数\( f(x; \theta) \)について、以下の統計検定を行う際の最強力検定に関する定理をネイマン-ピアソンの基本定理と呼ぶ。

      $$ H_0 : \theta = \theta_0 \\
      H_1 : \theta = \theta_1 $$

      ここでは有意水準を\( \alpha : 0 \leq \alpha \leq 1 \)とし、検定関数\( \varphi_0 (x) \)は以下の式で与えられる。
      $$ \varphi_0 (x) = \begin{cases}
      1 \ ( f(x ; \theta_1) > k f(x ; \theta_0) ) \\
      \gamma \ ( f(x ; \theta_1) = k f(x ; \theta_0) ) \\
      0 \ ( f(x ; \theta_1) < k f(x ; \theta_0) )
      \end{cases}$$
      ここで定数
      $$ \gamma : 0 \leq \gamma \leq 1 \\
      k : k \geq 0 $$
      は以下の式から定まるものとする。
      $$ E_{\theta_0} \left[ \varphi_0(X) \right ] = \alpha $$

      ネイマン-ピアソンの基本定理の証明について

      任意の統計検定の関数について、ネイマン-ピアソンの基本定理にしたがった統計検定の検出力を上回らないことを示せば良い。

      Xが連続であるとして、
      $$
      C = \{ x \in \mathfrak X | f(x; \theta_1) > k f(x; \theta_0) \} \\
      D = \{ x \in \mathfrak X | f(x; \theta_1) = k f(x; \theta_0) \} \\
      E = \{ x \in \mathfrak X | f(x; \theta_1) < k f(x; \theta_0) \}
      $$
      とおくと、これらは互いに素であるから、
      $$ \mathfrak X = C \cup D \cup E $$
      と表される。

      ここで統計検定の関数\( \varphi (X) \)を
      $$
      H_0 : \theta = \theta_0 \\
      H_1 : \theta = \theta_1
      $$
      の検定問題に対する有意水準\( \alpha \)の任意の検定とする。
      このとき、定義から\( \varphi (X) \)は
      $$ E_{\theta_0} [ \varphi (X)] \leq \alpha $$
      を満たす。
      \( \varphi_0 (X) \)の定義から、

      $$
      E_{\theta_1}[ \varphi_0 (X) ] – E_{\theta_1}[ \varphi (X) ] \\
      = \int_{\mathfrak X} \varphi_0 (X) f(x ; \theta_1) dx – \int_{\mathfrak X} \varphi (X) f(x ; \theta_1) dx \\
      = \int_{\mathfrak X} (\varphi_0 (X) – \varphi (X) ) f(x ; \theta_1) dx \\
      = \int_{C} (1 – \varphi (X) ) f(x ; \theta_1) dx + \int_{D} (\gamma – \varphi (X) ) f(x ; \theta_1) dx + \int_{E} ( – \varphi (X) ) f(x ; \theta_1) dx \\
      \geq \int_{C} (1 – \varphi (X) ) (k f(x;\theta_0))dx + \int_{D} (\gamma – \varphi (X) ) (k f(x;\theta_0))dx + \int_{E} ( – \varphi (X) ) (k f(x;\theta_0))dx \\
      = k \int_{C} (1 – \varphi (X) )f(x;\theta_0)dx + k \int_{D} (\gamma – \varphi (X) )f(x;\theta_0)dx + k \int_{E} ( – \varphi (X) )f(x;\theta_0)dx \\
      = k \int_{\mathfrak X} (\varphi_0 (X) – \varphi (X) ) f(x ; \theta_1) dx \\
      = k \{E_{\theta_0}[ \varphi_0 (X) ] – E_{\theta_0}[ \varphi (X) ]\} \\
      = k \{ \alpha – E_{\theta_0}(\varphi(X))\} \\
      \geq 0$$

      よって、
      $$ E_{\theta_1}[ \varphi_0 (X) ] \geq E_{\theta_1}[ \varphi (X) ] $$
      となり、\( \varphi_0 (X) \)が任意のものよりも大きな検出力を持つことから、有意水準\( \alpha \)の最強力検定であることが示された。

      ネイマン-ピアソンによる最強力検定の例(正規分布)

      平均が未知で分散が既知のケースにおいて、確率変数が\( X_1, \dots , X_n \sim N(\mu, \sigma_0^2) \)に従うとして、以下の検定問題を考え、有意水準\( \alpha \)の最強力検定を求める。
      
      $$
      H_0 : \mu = \mu_0 \\
      H_1 : \mu < \mu_0
      $$
      これは単純仮説と単純対立仮説からなる検定問題なので、以下のように書き換えることができる。
      $$
      H_0 : \mu = \mu_0 \\
      H_1 : \mu = \mu_1 \\
      (\mu_1 < \mu_0)
      $$

      なお、この例は片側検定問題となる。
      \( X_i \)の確率密度関数を\( f(x_i ; \mu) \)として、ネイマン-ピアソンの基本定理に従うと、有意水準\( \alpha \)の最強力検定は以下のようになる。

      $$ \varphi_0 (x) = \begin{cases}
      1 \ ( \prod_{i=1}^n f(x_i ; \mu_1) > k \prod_{i=1}^n f(x_i ; \mu_0) ) \\
      \gamma \ ( \prod_{i=1}^n f(x_i ; \mu_1) = k \prod_{i=1}^n f(x_i ; \mu_0)) \\
      0 \ ( \prod_{i=1}^n f(x_i ; \mu_1) < k \prod_{i=1}^n f(x_i ; \mu_0) )
      \end{cases}$$

      ここで、
      $$
      \gamma \in [0, 1] \\
      k \in [0, \infty ]
      $$

      $$ E_{\mu_0} [ \varphi_0(X) ] = \alpha $$
      を満たす定数で与えられるものとする。

      \( X = (X_1, \dots , X_n ) \)が連続型の場合、

      $$P_{\mu_j} \left \{ \prod_{i=1}^n f(x_i ; \mu_1) = k \prod_{i=1}^n f(x_i ; \mu_0) \right \} = 0 \\
      j = 0, 1$$

      であるから、このケースは無視することができる。
      仮定より、
      $$ \prod_{i=1}^n f(x_i ; \mu) = \prod_{i=1}^n \frac{1}{\sqrt {2\pi} \sigma_0} e^{-\frac{(x_i – \mu)^2}{2\sigma_0^2}} \\
      = \frac{1}{(\sqrt {2\pi})^n \sigma_0^n} e^{-\frac{\sum_{i=1}^n (x_i – \mu)^2 }{2\sigma_0^2}} \\
      = \frac{1}{(\sqrt {2\pi}\sigma_0)^n } e^{-\frac{n(\bar x – \mu)^2 + \sum_{i=1}^n(x_i – \bar x)^2}{2\sigma_0^2}}
      $$
      よって検定関数において、棄却するケースについて以下のように書き表すことができる。
      $$\prod_{i=1}^n f(x_i ; \mu_1) > k \prod_{i=1}^n f(x_i ; \mu_0) \\
      \Leftrightarrow \frac{\prod_{i=1}^n f(x_i ; \mu_1)}{\prod_{i=1}^n f(x_i ; \mu_0)} > k \\
      \Leftrightarrow \frac{\frac{1}{(\sqrt {2\pi}\sigma_0)^n } e^{-\frac{n(\bar x – \mu_1)^2 + \sum_{i=1}^n(x_i – \bar x)^2}{2\sigma_0^2}}}{\frac{1}{(\sqrt {2\pi}\sigma_0)^n } e^{-\frac{n(\bar x – \mu_0)^2 + \sum_{i=1}^n(x_i – \bar x)^2}{2\sigma_0^2}}} \\
      \Leftrightarrow e^{\frac{n}{2\sigma_0^2} ((\bar x – \mu_0)^2 – (\bar x – \mu_1)^2) } > k \\
      \Leftrightarrow (\bar x – \mu_0)^2 – (\bar x – \mu_1)^2 > \frac{2\sigma_0^2}{n} \log k = k’
      $$
      ここで、\( \mu_1 < \mu_0 \)より
      $$(\bar x – \mu_0)^2 – (\bar x – \mu_1)^2 < 0$$
      であるから、
      $$ \bar x < \frac{\mu_1^2 – \mu_0^2}{2(\mu_1 – \mu_0)} = k” $$
      となる。
      一方、\( H_0 \)が真のとき、\( \bar X \)は\( N \left ( \mu_0, \frac{\sigma_0^2}{n} \right ) \)に従うので、
      $$\alpha = E_{\mu_0}[ \varphi_0 (X) ] \\
      = P_{\mu_0} \left \{ \prod_{i=1}^n f(x_i ; \mu_1) > k \prod_{i=1}^n f(x_i ; \mu_0) \right \} \\
      = P_{\mu_0} \left \{ \bar X < k” \right \} \\
      = \int_{-\infty}^{k”}\frac{\sqrt{n}}{\sqrt{2\pi}\sigma_0}e^{-\frac{n(t-\mu_0)^2}{2\sigma_0^2}}dt
      $$
      ここで\( \frac{\sqrt{n}(t-\mu_0)}{\sigma_0} = u \)、\( k_0 = \frac{\sqrt{n}(k” – \mu_0)}{\sigma_0} \)とおいて変数変換を行うと、
      $$ \alpha = \int_{-\infty}^{k_0} \frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du \\
      = \Phi (k_0)
      $$
      となる。\(\Phi(\cdot)\)は\( N(0,1) \)の分布関数を表すので、標準正規分布表から\( \Phi (k_0) = \alpha \)をみたす\( k_0 \)の値を読み取ることで、以下のような今回の統計検定の棄却域が得られる。
      $$k^{*} = \mu_0 + \frac{\sigma_0}{\sqrt{n}}k_0 = \mu_0 – \frac{\sigma_0}{\sqrt{n}}z(\alpha) $$
      なお、これは\( \mu_1 \)の水準に無関係に決まる。

      ネイマン-ピアソンによる最強力検定の例(ベルヌーイ分布)

      確率pが未知で確率変数mが既知のケースにおいて、確率変数が\( X_1, \dots , X_n \sim B_N(m, p) \)として、以下の検定問題を考え、有意水準\( \alpha \)の最強力検定を求める。
      $$
      H_0 : p = p_0 \\
      H_1 : p < p_0
      $$
      これは単純仮説と単純対立仮説からなる検定問題なので、以下のように書き換えることができる。
      $$
      H_0 : p = p_0 \\
      H_1 : p = p_1 \\
      ( p_1 < p_0)
      $$

      仮定より、\( X = ( X_1, \dots , X_n ) \)の確率関数は
      $$\prod_{i=1}^n f(x_i ; p) = \prod_{i=1}^n {}_m C _{x_i}p^{x_i}(1-p)^{m-x_i} \\
      = p^{\sum_{i=1}^n x_i}(1-p)^{mn – \sum_{i=1}^n x_i} \prod_{i=1}^n {}_m C _{x_i} $$
      となる。
      対立仮説が正しいのであれば、尤度の比において以下の関係となる。

      $$\frac{\prod_{i=1}^n f(x_i ; p_1)}{ \prod_{i=1}^n f(x_i ; p_0) } > k \\
      \Leftrightarrow \frac{ p_1^{\sum_{i=1}^n x_i}(1-p_1)^{mn – \sum_{i=1}^n x_i} \prod_{i=1}^n {}_m C _{x_i}}{p_0^{\sum_{i=1}^n x_i}(1-p_0)^{mn – \sum_{i=1}^n x_i} \prod_{i=1}^n {}_m C _{x_i}} > k \\
      \Leftrightarrow \left ( \frac{p_1(1-p_0)}{p_0(1-p_1)} \right )^{\sum_{i=1}^nx_i} \times \left ( \frac{1-p_1}{1-p_0} \right )^{mn} > k \\
      \Leftrightarrow \left ( \frac{p_1(1-p_0)}{p_0(1-p_1)} \right ) ^{ \sum_{i=1}^n x_i } > k’ = \left ( \frac{1-p_0}{1-p_1} \right )^{mn}k $$

      \( p_1 < p_0 \)であることから、
      $$\left ( \frac{p_1(1-p_0)}{p_0(1-p_1)} \right ) < 1$$
      となり、1よりも小さな項の指数なので、これを満たす\( \sum_{i=1}^n x_i \)であるために、
      $$\Leftrightarrow \sum_{i=1}^n x_i < k”$$
      が満たされる必要がある。

      よってこの検定問題における有意水準\( \alpha \)の最強力検定はネイマン-ピアソンの基本定理に従うと、
      $$ \varphi_0 (x) = \begin{cases}
      1 \ ( \sum_{i=1}^n < k” ) \\
      \gamma \ ( \sum_{i=1}^n = k” ) \\
      0 \ ( \sum_{i=1}^n > k” )
      \end{cases}
      $$
      によってなされる。
      なお、
      $$ \gamma \in [0, 1] \\
      k” \in [0, \infty ] $$

      $$E_{p_0}[ \varphi_0 (X) ] = \alpha $$
      を満たす定数によって与えられる。

      帰無仮説\( H_0 \)が真のとき\( \sum_{i=1}^n X_i \)は\( B_N(mn, p_0) \)に従うので、\( \gamma, k” \)は以下の式を解くことで求められる。

      $$ E_{p_0} [ \varphi_0 (x)] = 1 \times P_{p_0} \left \{ \sum_{i=1}^n x_i < k” \right \} \\
      + \gamma \times P_{p_0} \left \{ \sum_{i=1}^n x_i = k” \right \} \\
      + 0 \times P_{p_0} \left \{ \sum_{i=1}^n x_i > k” \right \} \\
      = \sum_{j < k”} {}_{mn} C _j p_0^j (1-p_0)^{mn-j} + \gamma {}_{mn} C _{k”} p_0^{k”} (1- p_0)^{mn – k”} \\
      = \alpha $$

      \( B_N(mn, p_0) \)の分布関数を\( F(x ; p_0) \)としたとき、
      $$F(k-1; p_0) \leq \alpha $$
      を満たす最大の整数kを\( k” \)として、
      $$ \gamma = \frac{\alpha – F(k” – 1 ; p_0 )}{F(k”;p_0) – F(k” – 1 ; p_0)} $$
      ただし、\( F(0 ; p_0) \leq \alpha \)とする。\( F(0 ; p_0) \geq \alpha \)の場合は、
      $$k” = 0 \\
      \gamma = \frac{\alpha}{F(0; p_0)}$$
      となる。これらは対立仮説の\(p_1\)に関係なく決まるため、一様最強力検定となる。

      不偏検定

      単純仮説ではなく複合仮説の場合、先ほどのような検定を不偏検定と呼ぶ。
      検定問題は以下のようになる。

      $$ H_0 : \theta \in \Theta_0 \\
      H_1 : \theta \in \Theta_1 \ ( = \Theta – \Theta_0 ) $$

      この検定問題で水準\( \alpha \)の検定\( \varphi (X) \)でその検出力が\( \alpha \)以上、
      $$ \beta (\theta ; \varphi ) := E_{\theta} [ \varphi (X)] \geq \alpha, \forall \theta \in \Theta_1 $$
      を満たす検定\( \varphi (X) \)を水準\( \alpha \)の不偏検定と呼び、水準\( \alpha \)の不偏検定の中で検出力を一様に最大にするものを水準\( \alpha \)の一様最強力不偏検定と呼ぶ。

      不偏検定の際の検出力関数は以下の図のようになる。検出力関数がわかれば、\( \alpha \)(第一種の誤りをする確率)を与えたもとで棄却域が決まる。

      不偏検定に関する定理

      まず前提として、\( \theta \)を母数として、\( X = (X_1, \dots , X_n) \)の確率関数\( f(x, \theta) \)が以下の表現で表されるとする。
      $$ f(x, \theta) = C(\theta)e^{\theta T(x)}h(x) $$

      \( \theta_0, \theta_1, \theta_2 (\theta_1 < \theta_2) \)を与えられたものとして、以下の複合仮説の仮説検定問題を考える。

      (1)
      $$ H_0 : \theta_1 \leq \theta \leq \theta_2 \\
      H_1 : \theta < \theta_1 \ or \ \theta > \theta_2
      $$
      (2)
      $$ H_0 : \theta = \theta_0 \\
      H_1 : \theta \neq \theta_0
      $$
      このような仮説検定問題を両側検定問題と呼ぶ。

      検定問題(1),(2)に対する一様最強力不偏検定\( \varphi_0(X) \)は以下の式で与えられる。

      $$
      \varphi_0 (X) =
      \begin{cases}
      1 \ ( T(x) < k_1 \ or \ T(x) > k_2 ) \\
      \gamma \ ( T(x) = k_i , i=1, 2 ) \\
      0 \ (k_1 < T(x) < k_2)
      \end{cases} \\
      0 \leq \varphi (x) \leq 1
      $$

      ここで、定数\( \gamma_i \in [0,1], k_i \geq 0, (i = 1, 2) \)は検定問題(1)については以下の式により定まる。

      $$E_{\theta_1}[ \varphi_0 (X) ] = E_{\theta_2} [\varphi_0 (X)] = \alpha $$
      一方、検定問題(2)については以下の式により定まる。
      $$ E_{\theta_0}[\varphi_0 (X)] = \alpha \\
      E_{\theta_0}[T(X) \varphi_0 (X)] =E_{\theta_0}[T(X)] \alpha
      $$

      以前のネイマン-ピアソンの基本定理は確率の積和(尤度)があったが、この定理ではより一般的な記述がされており、帰無仮説や対立仮説のパラメータも記されていない。しかしながら、非常に似ている。

      不偏検定の例(正規分布(平均が未知、分散が既知のケース))

      確率変数が以下に従うと仮定する。
      $$X_1, X_2, \dots , X_n \sim N(\mu, \sigma_0^2)$$
      ここで以下の検定問題を考える。
      $$ H_0 : \mu = \mu_0 \\
      H_1 : \mu \neq \mu_0 $$
      この検定問題に対する、一様最強力不偏検定を行う。
      仮定より、\( X = (X_1,\dots , X_n) \)の確率密度関数の積は
      $$ \prod_{i=1}^n f(x_i ; \mu) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma_0}e^{-\frac{(x_i – \mu)^2}{2\sigma_0^2}} \\
      = \frac{1}{(\sqrt{2\pi}\sigma_0)^n}e^{-\frac{(x_i – \mu)^2}{2\sigma_0^2}} \\
      = \frac{1}{(\sqrt{2\pi}\sigma_0)^n} e^{-\frac{n\mu^2}{2\sigma_0^2}}e^{\mu \frac{\sum_{i=1}^n x_i}{\sigma_0^2}}e^{- \frac{\sum_{i=1}^n x_i^2}{2\sigma_0^2}} $$
      ここで、統計量として以下を定義する。
      $$ T(x) = \frac{1}{\sigma_0^2} \sum_{i=1}^n x_i $$
      これを用いると先ほどの尤度は以下のような指数型分布になる。
      $$ \prod_{i=1}^n f(x_i ; \mu) = C(\mu) e^{\mu T(x)}h(x) $$
      不偏検定に関する定理を用いることで、水準\( \alpha \)の一様最強力不偏検定の棄却域は以下で与えられる。

      $$ C = \left \{ x = (x_1, \dots, x_n) \in \mathbb R^n | \frac{1}{\sigma_0^2} \sum_{i=1}^n x_i < nk_1 \ or \ \frac{1}{\sigma_0^2} \sum_{i=1}^n x_i > nk_2 \right \} $$

      ここで、\(k_1, k_2\)は、
      $$ P_{\mu_0} \{ C \} = \alpha $$
      および、

      $$ E_{\mu_0} \left \{ \frac{1}{\sigma_0^2} \sum_{i=1}^n x_i \cdot 1_c(x) \right \} = E_{\mu_0} \left \{ \frac{1}{\sigma_0^2} \sum_{i=1}^n x_i \right \} \cdot \alpha $$

      を満たす定数となる。(\( 1_c(x) \)はインデックス関数。)

      各\( X_i \)が独立に\( N(\mu, \sigma_0^2) \)に従い、\( \bar X = \frac{1}{n} \sum_{i=1}^n \)が\( N(\mu, \frac{\sigma_0^2}{n}) \)に従う。

      ここで水準\( \alpha \)を以下のように書くことができる。

      $$ \alpha = P_{\mu_0} \left \{ \bar X < k_1 \sigma_0^2 \ or \ \bar X > k_2 \sigma_0^2 \right \} $$

      ここで、
      $$ k’_i = \frac{\sqrt {n} (k_i \sigma_0^2) – \mu_0}{\sigma_0}, i =1, 2 $$
      とおくと、

      $$ \alpha = P_{\mu_0} \left \{ \frac{\bar X – \mu_0}{ \sqrt{\frac{\sigma_0^2}{n}}} < k’_1 \right \} + P_{\mu_0} \left \{ \frac{\bar X – \mu_0}{ \sqrt{\frac{\sigma_0^2}{n}}} > k’_2 \right \} \\
      = \Phi (k’_1) + 1 – \Phi (k’_2) $$

      となる。
      他方、不偏検定に関する定理より、
      $$ E_{\mu_0} \left \{ \bar X \cdot 1_C(X) \right \} = E_{\mu_0} \left \{ \bar X \right \}\alpha $$
      ここで、
      $$ z = \frac{\sqrt{n}(\bar X – \mu_0)}{\sigma_0} $$
      とおくと、\( \bar X = \mu_0 + \frac{\sigma_0}{\sqrt{n}}Z \)なので、

      $$ E_{\mu_0} \left \{ \bar X \cdot 1_C(X) \right \} = E_{\mu_0} \left \{ ( \mu_0 + \frac{\sigma_0}{\sqrt{n}}Z ) \cdot 1_C(X) \right \} \\
      = \mu_0 E_{\mu_0} \left \{ 1_C(X) \right \} + \frac{\sigma_0}{\sqrt{n}} E_{\mu_0} \left \{z \cdot 1_C(X) \right \} \\
      = \alpha \mu_0 + \frac{\sigma_0}{\sqrt{n}} E_{\mu_0} \left \{z \cdot 1_C(X) \right \}$$

      よって、
      $$ E_{\mu_0} \left \{ z \cdot 1_C(X) \right \} = 0 $$
      となる。
      \( C = \left \{ z < k’_1 \ or \ z > k’_2 \right \} \)と書けることから、これを書き換えると、
      $$ \int_{-\infty}^{k’_1}z \phi (z) dz + \int_{k’_2}^{+\infty}z \phi (z) dz = 0 $$
      となる。
      すなわち、
      $$ \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{(k’_1)^2}{2}} = \frac{1}{\sqrt{2 \pi}} e^{ – \frac{(k’_2)^2}{2} } $$
      となる。これより、
      $$ | k’_1 | = | k’_2 | $$
      から、
      $$ k’_1 = -k’_2 < 0 $$
      となる。
      以上から、水準\( \alpha \)の一様最強力不偏検定の棄却域は
      $$ \bar X < k_1\sigma_0^2 \\
      \bar X > k_2 \sigma_0^2 $$
      で与えられる。
      ここで、
      $$ k_1 = \left ( \mu_0 + \frac{\sigma_0}{\sqrt{n}}k’_1 \right ) / \sigma_0^2 \\
      k_2 = \left ( \mu_0 – \frac{\sigma_0}{\sqrt{n}}k’_1 \right ) / \sigma_0^2 $$
      となり、\( k’_1 \)は
      $$ \Phi ( k’_1 ) = \frac{\alpha}{2} $$
      を満たす、標準正規分布の値、\( k’_1 = -z \left ( \frac{\alpha}{2} \right ) \)となる。

      不偏検定の例(正規分布(平均が既知、分散が未知のケース))

      確率変数が以下に従うと仮定する。
      $$X_1, X_2, \dots , X_n \sim N(\mu_0, \sigma^2)$$
      ここで以下の検定問題を考える。
      $$ H_0 : \sigma^2 = \sigma_0^2 \\
      H_1 : \sigma^2 \neq \sigma_0^2 $$
      この検定問題に対する、一様最強力不偏検定を行う。
      仮定より、\( X = (X_1,\dots , X_n) \)の確率密度関数の積は
      $$ \prod_{i=1}^n f(x_i ; \sigma^2) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x_i – \mu_0)^2}{2\sigma^2}} \\
      = \frac{1}{(\sqrt{2\pi}\sigma)^n}e^{-\frac{(x_i – \mu_0)^2}{2\sigma^2}} $$
      ここで、統計量として以下を定義する。
      $$ T(x) = \frac{ \sum_{i=1}^n (x_i – \mu_0)^2 }{\sigma_0^2}$$
      さらに、
      $$ \theta = – \frac{\sigma_0^2}{2\sigma_2} $$
      と考えると、先ほどの尤度は以下のような指数型分布になる。
      $$ \prod_{i=1}^n f(x_i ; \sigma^2) = C(\theta) e^{\theta T(x)}h(x) $$
      不偏検定に関する定理を用いることで、水準\( \alpha \)の一様最強力不偏検定の棄却域は以下で与えられる。

      $$ C = \left \{ x = (x_1, \dots, x_n) \in \mathbb R^n | T(x) < k_1 \ or \ T(x) > k_2 \right \} $$

      ここで、\(k_1, k_2\)は、
      $$ P_{\sigma_0^2} \{ C \} = \alpha $$
      および、
      $$ E_{\sigma_0^2} \left \{ T(x) \cdot 1_c(x) \right \} = E_{\sigma_0^2} \left \{ T(x) \right \} \cdot \alpha $$
      を満たす定数となる。(\( 1_c(x) \)はインデックス関数。)
      \(H_0\)が真のとき、\( T(x)= \frac{\sum_{i=1}^n (x_i – \mu_0)^2 }{\sigma_0^2} \)は自由度nの\( \chi ^2 \)分布に従うので、
      $$ P_{\sigma_0^2} \{ C \} = P_{\sigma_0^2} \{ T(x) < k_1 \ or \ T(x) > k_2 \} \\
      = \int_0^{k_1}g_n (x) dx + \int_{k_2}^{\infty}g_n (x) dx = \alpha $$
      となる。

      $$ E_{\sigma_0^2} \{ T(x) \cdot 1_c(x) \} = \int_{0}^{k_1}x g_n(x)dx + \int_{k_2}^{\infty}x g_n(x)dx $$

      ここで、

      $$ \int_{a}^{b}xg_n(x)dx = \int_{a}^{b} \frac{1}{\Gamma \left (
      \frac{n}{2} \right )2^{\frac{n}{2}} }x^{\frac{n}{2}}e^{-\frac{x}{2}} \\
      = \frac{\Gamma \left ( \frac{n+2}{2} \right ) 2^{\frac{n+2}{2}} }{\Gamma \left ( \frac{n}{2} \right ) 2^{\frac{n}{2}}} \int_{a}^{b}\frac{1}{\Gamma \left ( \frac{n+2}{2} \right ) 2^{\frac{n+2}{2}}} x^{\frac{n+2}{2}-1}e^{-\frac{x}{2}}dx \\
      = n \int_{a}^{b} g_{n+2}(x)dx $$

      とすると、先ほどの式は、

      $$ E_{\sigma_0^2} \{ T(x) \cdot 1_c(x) \} = n \int_{0}^{k_1}g_{n+2}(x)dx + n \int_{k_2}^{\infty}g_{n+2}(x)dx \\
      = n \alpha$$

      となる。以上より、一様最強力不偏検定の棄却域\( k_1, k_2 \)は、
      $$\int_{0}^{k_1}g_n(x)dx + \int_{k_2}^{\infty}g_n(x)dx = \alpha \\
      \int_{0}^{k_1}g_{n+2}(x)dx + \int_{k_2}^{\infty}g_{n+2}(x)dx = \alpha$$
      を同時に満たすような\( k_1, k_2 \)を選べば良い。
      実際には簡単にするために、
      $$ \int_{0}^{k_1}g_{n}(x) dx = \int_{k_2}^{\infty}g_{n}(x) dx = \frac{\alpha}{2} $$
      となるような\( k_1, k_2 \)を\( \chi^2 \)分布表から読み取って、
      $$ k_1 = \chi^2 \left ( 1 – \frac{\alpha}{2}; n \right ) \\
      k_2 = \chi^2 \left ( \frac{\alpha}{2}; n \right ) $$
      によって近似的に水準\( \alpha \)の一様最強力不偏検定の棄却域として、
      $$ \frac{\sum_{i=1}^n (x_i – \mu_0)^2 }{\sigma_0^2} < k_1 \\
      \frac{\sum_{i=1}^n (x_i – \mu_0)^2 }{\sigma_0^2} > k_2 $$
      が採用される。

      尤度比検定

      仮説検定論における普遍的な検定方式のこと。帰無仮説に従う際の尤度と、真のパラメータに従う際の尤度の比をもとに検定を行う。

      確率ベクトル(\( X = (X_1, \dots , X_n) \))が確率密度関数(\( f(X_1, \dots , X_n; \theta) \))に従うとして、
      真のパラメータ(\( \theta \in \Theta \))、
      帰無仮説に従うパラメータ(\( \theta_0 \ ( \neq \phi ) \ \in \Theta \))、
      対立仮説に従うパラメータ(\( \theta_1 = \Theta – \Theta_0 \ ( \neq \phi ) \))
      に対して、以下の仮説検定問題を考える。
      $$ H_0 : \theta \in \Theta_0 \\
      H_1 : \theta \in \Theta_1 $$

      ここで、固定された各標本(\( x = (x_1, \dots , x_n) \))に対して、尤度の比を求める。
      $$ \lambda (x) := \frac{ \sup_{\theta \in \Theta_0} f(x_1, \dots ,x_n ; \theta) }{ \sup_{\theta \in \Theta} f(x_1, \dots ,x_n ; \theta) } $$
      適当に定められた定数cに対して、
      $$ \lambda (x) < c $$
      となるとき、仮説\( H_0 \)を棄却し、そうでないとき\( H_0 \)を受容するという検定方式を考える。ここでcは
      $$ \sup_{\theta \in \Theta_0} P_{\theta} \left \{ \lambda (X) < c \right \} = \alpha $$
      を満たすものとする。このような検定方式を水準\( \alpha \)の尤度比検定と呼び、\( \lambda (x) \)を尤度比統計量と呼ぶ。

      尤度比検定は複合仮説からなる検定問題において、\( \Theta_0, \Theta_1 \)をそれぞれ1つの母数で代表させることで、複合仮説検定問題を単純仮説と単純対立仮説からなる検定問題に置き換えることができる。

      他にも尤度比検定は、データが十分に大きい際に、検出力が1に収束するという一致性や、漸近的に最大の検出力を持つという漸近有効性などの性質も持ち合わせている。

      尤度比検定に関する定理

      \( n \to + \infty \)とすると、統計量\( -2 log \lambda (x) \)は仮説\( H_0 \)のもとで自由度qの\( \chi^2 \)分布に収束する。ここでは\( g(t) \)がその確率密度関数となる。
      $$ \theta \in \Theta_0 \\
      \lim_{n \to \infty} P_{\theta} \left \{ -2 \log \lambda (x) \leq z \right \} \\
      = \int_{0}^{z}g(t)dt \\
      \forall z \in (0, \infty) $$

      尤度比検定(1標本のケース)

      ここではデータが1つの標本から独立にそれぞれ\( N(\mu, \sigma^2) \)に従うケースを想定する。
      その上で、
      「\(\mu\)が未知で\(\sigma^2 = \sigma_0^2 \)が既知の場合」
      「\(\mu, \sigma^2\)ともに未知の場合」についての尤度比検定を行う。

      \(\mu\)が未知で\(\sigma^2 = \sigma_0^2 \)が既知の場合

      ここでは以下の検定問題を扱う。
      $$ H_0 : \mu = \mu_0 \\
      H_1 : \mu \neq \mu_0 $$
      仮定から、確率ベクトル(標本確率変数)である\( X = (X_1, \dots, X_n) \)の確率密度関数は
      $$ f(x_1, \dots, x_n ; \theta) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma_0} e^{- \frac{(x_i – \mu)^2}{2\sigma_0^2}} \\
      = \left ( \frac{1}{\sqrt{2\pi} \sigma_0} \right )^n e^{ – \frac{\sum_{i=1}^n(x_i – \bar x)^2 + n(\bar x – \mu)^2 }{2\sigma_0^2} } $$

      \(\mu\)の\( \Theta = \mathbb R^1 \)における最尤解は\( \hat \mu = \bar X \)となるので、尤度比は帰無仮説に従う確率密度関数と真のパラメータの従う確率密度関数であるから、
      $$ \lambda (x) = \frac{f(x_1, \dots, x_n ; \mu_0)}{f(x_1, \dots, x_n ; \hat \mu)} \\
      = \frac{ \left ( \frac{1}{\sqrt{2\pi} \sigma_0} \right )^n e^{ – \frac{ \sum_{i=1}^n(x_i – \bar x)^2 + n(\bar x – \mu_0)^2 }{2\sigma_0^2} }}{ \left ( \frac{1}{\sqrt{2\pi} \sigma_0} \right )^n e^{ – \frac{\sum_{i=1}^n(x_i – \bar x)^2 + n(\bar x – \bar x)^2 }{2\sigma_0^2} }} \\
      = e^{- \frac{n(\bar x – \mu_0)^2}{2\sigma_0^2}} $$

      ここで、水準\(\alpha\)の尤度比検定の棄却域は最後の式の指数部分の
      $$ | \bar x – \mu_0 | > c$$
      で与えられる。
      定数cは
      $$ P_{\mu_0} \{ | \bar X – \mu_0 | > c \} = \alpha $$
      より定まる。

      帰無仮説\( H_0 \)のもとで\(\bar X\)は\( N \left ( \mu_0 , \frac{\sigma_0^2}{n} \right ) \)に従うので、
      $$ \alpha = P_{\mu_0} { \frac{ \sqrt{n} | \bar X – \mu_0 |}{\sigma_0} > \frac{\sqrt{n}c}{\sigma_0} } \\
      = \int_{k}^{\infty} \phi (t)dt + \int_{-\infty}^{-k} \phi (t) dt \\
      = 2 \left ( 1 – \Phi (k) \right ) $$

      したがって、標準正規分布から水準\(\alpha\)を満たすようなkの値を見つければよい。

      \(\mu, \sigma^2\)ともに未知の場合

      ここでは以下の二つの検定問題を扱う。
      「平均に関する検定」
      「分散に関する検定」

      平均に関する検定

      $$ H_0 : \mu = \mu_0 \\
      H_1 : \mu \neq \mu_0 $$

      真のパラメータは以下で表される。
      $$ \theta = \left \{ (\mu, \sigma^2) | \mu \in \mathbb R^1 , \sigma^2 > 0 \right \} $$
      同様に帰無仮説に従う場合のパラメータは、
      $$ \theta_0 = \left \{ (\mu_0, \sigma^2) | \sigma^2 > 0 \right \} $$
      対立仮説に従う場合のパラメータは、
      $$ \theta_1 = \left \{ (\mu, \sigma^2) | \mu (\neq \mu_0) \in \mathbb R^1 , \sigma^2 > 0 \right \} $$
      で表される。

      仮定から、確率ベクトル(標本確率変数)である\( X = (X_1, \dots, X_n) \)の確率密度関数は
      $$ f(x_1, \dots, x_n ; \mu, \sigma^2) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma} e^{- \frac{(x_i – \mu)^2}{2\sigma^2}} \\
      = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{\sigma^2} \right )^{\frac{n}{2}} e^{ – \frac{\sum_{i=1}^n(x_i – \bar x)^2 + n(\bar x – \mu)^2 }{2\sigma^2} } $$

      ここで、\( \theta = (\mu, \sigma^2) \)が\( \Theta_0\)上を動くときに、\( \mu = \mu_0\)なので、尤度を最大にする問題において以下のように書くことができる。

      $$ \sup_{\theta \in \Theta_0} f(x_1, \dots, x_n ; \mu, \sigma^2) = \sup_{\sigma^2 > 0} f(x_1, \dots, x_n ; \mu_0, \sigma^2) \\
      = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat \sigma^2} \right )^{\frac{n}{2}} e^{ – \frac{\sum_{i=1}^n(x_i – \mu_0)^2}{2 \hat \sigma^2} } $$

      ここで\( \hat \sigma^2 = \frac{1}{n} \sum_{i=1}^n (x_i – \mu_0)^2 \)とすると、
      $$ = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat \sigma^2} \right )^{\frac{n}{2}} e^{-\frac{n}{2}} $$
      一方、\( \theta = (\mu, \sigma^2) \)が真のパラメータ\( \Theta\)上を動くときに、最尤解\( (\hat \mu, \hat \sigma^2)\)が、以下で与えられる。
      $$ \hat \mu = \bar x \ \left ( \frac{1}{n} \sum_{i=1}^n \bar x \right ) \\
      \hat \sigma = \frac{1}{n}\sum_{i=1}^n (x_i – \bar x)^2 $$
      よって尤度を最大にする問題において以下のように書くことができる。

      $$ \sup_{\theta \in \Theta} f(x_1, \dots, x_n ; \mu, \sigma^2) = f(x_1, \dots, x_n ; \hat \mu,\hat{\hat{\sigma^2}}) \\
      = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat{\hat{\sigma^2}}} \right )^{\frac{n}{2}} e^{ – \frac{\sum_{i=1}^n(x_i – \hat \mu)^2}{2 \hat{\hat{\sigma^2}}} } $$

      ここで、先ほどと同じように\( \hat{\hat{\sigma^2}} = \frac{1}{n} \sum_{i=1}^n (x_i – \hat \mu)^2 \)とすると、
      $$ = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat{\hat{\sigma^2}}} \right )^{\frac{n}{2}} e^{-\frac{n}{2}} $$
      よって尤度比は帰無仮説の確率密度関数と真のパラメータの確率密度関数の比であるから、
      $$ \lambda (x) = \frac{f(x_1, \dots , x_n ; \mu_0,\hat{\hat{\sigma^2}} )}{f(x_1, \dots , x_n ; \mu_0,\hat \sigma^2 )} \\
      = \frac{\left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat{\hat{\sigma^2}}} \right )^{\frac{n}{2}} e^{-\frac{n}{2}}}{\left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat \sigma^2} \right )^{\frac{n}{2}} e^{-\frac{n}{2}}} \\
      = \left ( \frac{ \hat{\hat{\sigma^2}} }{\hat \sigma^2} \right )^{-\frac{n}{2}}$$
      となる。

      $$\frac{ \hat{\hat{\sigma^2}} }{\hat \sigma^2}= \frac{\frac{1}{n}\sum_{i=1}^n (x_i – \mu_0)^2}{\frac{1}{n}\sum_{i=1}^n(x_i – \bar x)^2}\\
      = \frac{ \sum_{i=1}^n \left ( (x_i – \bar x) – (\mu_0 – \bar x ) \right )^2 }{\sum_{i=1}^n(x_i – \bar x)^2} \\
      = \frac{ \sum_{i=1}^n (x_i – \bar x)^2 – n(\bar x – \mu_0)^2 }{\sum_{i=1}^n(x_i – \bar x)^2} \\
      = 1 + \frac{n(\bar x – \mu_0)^2}{\sum_{i=1}^n(x_i – \bar x)^2} $$
      ここで、2項目に注目すると、尤度比は

      $$ \lambda (x) < c \\
      \Leftrightarrow \frac{n(\bar x – \mu_0)^2}{ \sum_{i=1}^n (x_i – \bar x)^2 } > c’ \\
      \Leftrightarrow \frac{\sqrt{n} | \bar x – \mu_0 |}{u} > c” \\
      u = \sqrt{ \frac{1}{n-1} \sum_{i=1}^n (x_i – \bar x)^2 }$$
      となる。
      これは、\( H_0 \)が真のとき、
      $$ \frac{\sqrt{n} | \bar x – \mu_0 |}{u} \sim t(n-1) $$
      となるため、t分布表から棄却域を得ることができる。
      実際、t分布の上側\( \frac{\alpha}{2} \times 100% \)点である、\( t \left (
      \frac{\alpha}{2} ; n-1 \right ) \)を求めて、
      $$ c” = t \left ( \frac{\alpha}{2} ; n-1 \right ) $$
      とおくと、水準\( \alpha \)の尤度比検定の棄却域は
      $$ \frac{\sqrt{n} | \bar X – \mu_0 | }{U} > t \left ( \frac{\alpha}{2} ; n-1 \right ) $$
      で与えられる。これは両側のt検定に他ならない。

      このように尤度比検定は何を検定したいかで最終的に用いる分布が異なる。正規分布のときもあれば、t分布のときもある。

      分散に関する検定

      $$ H_0 : \sigma = \sigma_0^2 \\
      H_1 : \sigma \neq \sigma_0^2 $$

      真のパラメータが以下に従うとする。
      $$ \hat \mu = \bar x \\
      \hat \sigma^2 = \frac{1}{n}\sum_{i=1}^n(x_i – \bar x)^2$$

      その際、真のパラメータに従う確率密度関数の最尤解は

      $$ \sup_{ (\mu, \sigma^2) \in \Theta } f(x_1, \dots, x_n ; \mu, \sigma^2) = f(x_1, \dots, x_n ; \hat \mu, \hat \sigma^2)\\
      = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat \sigma^2 } \right )^{\frac{n}{2}} e^{-\frac{n}{2}} $$

      他方、帰無仮説に従う確率密度関数の最尤解は

      $$ \sup_{ (\mu, \sigma^2) \in \Theta_0 } f(x_1, \dots, x_n ; \mu, \sigma^2) = f(x_1, \dots, x_n ; \hat \mu, \sigma_0^2)\\
      = \left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \sigma_0^2 } \right )^{\frac{n}{2}} e^{-\frac{n \hat \sigma^2}{2\sigma_0^2}} $$

      となる。

      よって、尤度比は以下のようになる。
      $$ \lambda (x) = \frac{f(x_1, \dots, x_n ; \hat \mu, \sigma_0^2) }{f(x_1, \dots, x_n ; \hat \mu, \hat \sigma^2)} \\
      = \frac{\left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \sigma_0^2 } \right )^{\frac{n}{2}} e^{-\frac{n \hat \sigma^2}{2\sigma_0^2}}}{\left ( \frac{1}{\sqrt{2\pi}} \right )^n \left ( \frac{1}{ \hat \sigma^2 } \right )^{\frac{n}{2}} e^{-\frac{n}{2}}} \\
      = \left ( \frac{\hat \sigma^2}{\sigma_0^2} \right )^{\frac{n}{2}}e^{-\left ( \frac{\hat \sigma^2}{\sigma_0^2} – 1 \right ) \frac{n}{2} }
      $$
      これより、定数cについての不等号を記すと、

      $$ \lambda (x) < c \\
      \Leftrightarrow \log \lambda (x) = – \frac{n}{2} \left \{ \frac{\hat \sigma^2}{\sigma_0^2} – 1 – \log \left ( \frac{\hat \sigma^2}{\sigma_0^2} \right ) \right \} < c’ $$

      となる。これは結局のところ、統計検定の棄却域である\( k_1, k_2 \)を定めるように決まるので、以下のようになる。

      $$ \frac{\hat \sigma^2}{\sigma_0^2} < k_1 \ or \ \frac{\hat \sigma^2}{\sigma_0^2} > k_2 \\
      \frac{1}{\sigma_0^2}\sum_{i=1}^{n}(x_i-\bar x)^2 < k’_1 \ or \ \frac{1}{\sigma_0^2}\sum_{i=1}^{n}(x_i-\bar x)^2 > k’_2
      $$

      これは\( \chi^2 \)分布であるから、
      帰無仮説\( H_0 \)が真のとき、\( \frac{1}{\sigma_0^2}\sum_{i=1}^{n}(x_i-\bar x)^2 \)は\( \chi^2(n-1) \)に従う。
      定数\( k’_1,k’_2 \)を
      $$ \alpha = \int_{0}^{k’_1}g(x)dx + \int_{k’_2}^{\infty}g(x)dx $$
      より求めると、水準\( \alpha \)の尤度比検定の棄却域は

      $$ \frac{1}{\sigma_0^2}\sum_{i=1}^{n}(x_i-\bar x)^2 < k’_1 \ or \ \frac{1}{\sigma_0^2}\sum_{i=1}^{n}(x_i-\bar x)^2 > k’_2 $$

      で与えられる。棄却限界である\( k’_1,k’_2 \)は近似的に数値解を求めることで得られる。

      また、標本数nがある程度大きいとき、尤度比検定に関する定理より、統計量\( -2\log \lambda (x)\)が\( \chi^2 (1) \)に近似的に従うことより、

      $$ -2 \log \lambda (x) = n \left \{ \frac{\hat \sigma^2}{\sigma_0^2} – 1 – \log \left ( \frac{\hat \sigma^2}{\sigma_0^2} \right ) \right \} > \chi^2(\alpha ; 1) $$

      を棄却域とする検定方法を考えることができる。

      尤度比検定(2標本のケース)

      (一部数式の導出の解読中)
      三つのケースについて扱う。
      「A:\( \mu_1, \mu_2 \)は未知、\( \sigma_1^2, \sigma_2^2 \)は既知のケース」
      「B:\( \mu_1, \mu_2 \)は未知、\( \sigma_1^2 = \sigma_2^2 = \sigma^2 \)であるが\( \sigma^2 \)の値は未知のケース」
      「C::\( \mu_1, \mu_2, \sigma_1^2, \sigma_2^2 \)の全てが未知のケース」
      当然、ケースによって尤度比検定で従う分布が違ってくる。

      \( \chi^2 \)適合度検定

      ここではサイコロを複数回繰り返し投げたときの出目の出現回数について扱う。各出目の出現総数は多項分布に従うので、以下の確率密度関数を考える。
      $$ f(n_1, \dots , n_k ; \theta) = \frac{n!}{n_1! \dots n_k!}\theta_1^{n_1} \dots \theta_k^{n_k} $$
      ここでパラメータは以下に従うとする。

      $$ \Theta = \left \{ \theta = (\theta_1, \dots, \theta_k) | \sum_{i=1}^k \theta_i = 1 , 0 < \theta_i < 1, i=1, \dots , k
      \right \} $$

      ここでは以下の検定問題を考える。
      $$ H_0 : \theta_i = \theta_{i0} \ , \ i = 1, \dots ,k \\
      H_1 : \theta_i \neq \theta_{i0} \ (1 \leq \forall i \leq k) $$

      尤度比検定を考える際に、真のパラメータにおける最尤推定量を
      $$ \hat \theta_n = ( \hat \theta_{1n}, \dots, \hat \theta_{kn}) \\
      = \left ( \frac{n_1}{n},\dots , \frac{n_k}{n} \right ) $$
      とする。よって尤度比は
      $$ \lambda (n_1, \dots, n_k) = \frac{f(n_1, \dots , n_k ; \theta_0)}{f(n_1, \dots , n_k ; \hat \theta_n)} \\
      = \frac{\frac{n!}{n_1! \dots n_k!}\theta_{10}^{n_1} \dots \theta_{k0}^{n_k}}{\frac{n!}{n_1! \dots n_k!}\theta_{1n}^{n_1} \dots \theta_{kn}^{n_k}} \\
      = \left ( \frac{\theta_{10}}{\hat \theta_{1n}} \right )^{n_1} \dots \left ( \frac{\theta_{k0}}{\hat \theta_{kn}} \right )^{n_k} $$
      となる。
      両辺で対数をとり-2を掛けると、

      $$ -2 \log \lambda (n_1, \dots , n_k) = 2 \sum_{i=1}^k n_i \{ \log \hat \theta_{in} – \log \theta_{i0} \} $$

      この式より、水準\( \alpha \)の尤度比検定の棄却域は、適当な定数cに対して以下の式により与えられる。

      $$ -2 \log \lambda (n_1, \dots , n_k) = 2 \sum_{i=1}^k n_i \{ \log \hat \theta_{in} – \log \theta_{i0} \} > c $$

      ここで、cは
      $$ P_{\theta_{0}} \{ -2 \log \lambda (n_1, \dots, n_k) > c \} = \alpha $$
      を満たすように決まる。

      $$ -2 \log \lambda (N_1, \dots , N_k) = 2 \sum_{i=1}^k N_i \{ \log \hat \theta_{in} – \log \theta_{i0} \} \\
      = 2 \sum_{i=1}^k N_i \log \left ( \frac{\hat \theta_{in}}{\theta_{i0}} \right ) \\
      = 2 \sum_{i=1}^k N_i \log \left \{ 1+ \left ( \frac{\hat \theta_{in} – \theta_{i0} }{\theta_{i0}} \right ) \right \} $$

      ここでテイラー近似を行うと、

      $$ = 2 \sum_{i=1}^k N_i \left \{ \frac{\hat \theta_{in} – \theta_{i0} }{\theta_{i0}} – \frac{ (\hat \theta_{in} – \theta_{i0})^2 }{2\theta_{i0}^2} + O_{p} \left ( | \hat \theta_{in} – \theta_{i0} |^2 \right ) \right \} \\
      = 2 \sum_{i=1}^k N_i \left \{ \frac{N_{i} – n \theta_{i0} }{n \theta_{i0}} – \frac{ (N_{i} – n \theta_{i0})^2 }{2n^2\theta_{i0}^2} \right \} + O_{p}(1) \\
      = 2 \sum_{i=1}^k N_i \left \{ \frac{ (N_i – n \theta_{i0})^2 + n\theta_{i0}(N_i – n\theta_{i0}) }{n \theta_{i0}} – \frac{(N_i – n\theta_{i0})^3 + n\theta_{i0} (N_i – n \theta_{i0})^2}{2n^2\theta_{i0}^2} \right \} + O_{p}(1) \\
      = \sum_{i=1}^k N_i \frac{ (N_i – n \theta_{i0})^2}{n \theta_{i0}} – \sum_{i=1}^k \frac{(N_i – n\theta_{i0})^3}{n^2\theta_{i0}^2} + O_{p}(1) $$

      となる。ここで2項目に着目すると、

      $$ \sum_{i=1}^{k} \frac{ |N_i – n \theta_{i0}|^3 }{n^2 \theta_{i0}^2} = \sum_{i=1}^{k} \frac{(N_i – n\theta_{i0})^2}{n\theta_{i0}} \cdot \frac{|N_i – n \theta_{i0}|}{n\theta_{i0}} \\
      = \sum_{i=1}^{k} \frac{(N_i – n\theta_{i0})^2}{n\theta_{i0}} \cdot \frac{|\hat \theta_{i0} – \theta_{i0}|}{\theta_{i0}} \\
      \leq \epsilon_n \cdot \sum_{i=1}^{k} \frac{(N_i – n\theta_{i0})^2}{n\theta_{i0}} $$

      となる。ここで、
      $$ \epsilon_n = \max_{1 \leq i \leq k} \left \{ \frac{| \hat \theta_{i0} – \theta_{i0} |}{\theta_{i0}} \right \} $$
      かつ、
      $$ \epsilon_n = O_{p} (1) \ (n \to \infty) $$
      であることから、

      $$ -2 \log \lambda (N_1, \dots , N_k) = \sum_{i=1}^k N_i \frac{ (N_i – n \theta_{i0})^2}{n \theta_{i0}} + O_{p}(1) $$

      が成り立つ。
      これは\( \chi^2 \)分布に従うことから、帰無仮説\( H_0 \)のもとで、
      $$-2 \log \lambda (N_1, \dots , N_k) $$
      は\( n \to \infty \)のとき自由度k-1の\( \chi^2 \)分布に近似的に従う。

      よって、
      $$ \chi^2 = \sum_{i=1}^k N_i \frac{ (N_i – n \theta_{i0})^2}{n \theta_{i0}} > \chi^2 (\alpha ; k-1) $$
      のとき帰無仮説\( H_0 \)を棄却する。このような検定を、水準\( \alpha \times 100% \)の\( \chi^2 \)適合度検定と呼ぶ。これを実際の観測値に置き換えたものを修正\( \chi^2 \)検定と呼び以下のように表す。
      $$ \chi_0^2 = \sum_{i=1}^k \frac{(N_i – n \theta_{i0})^2}{N_i} $$

      分割表の検定

      (更新予定)

    参考情報

    [1]鈴木・山田(1996),『数理統計学―基礎から学ぶデータ解析』,内田老鶴圃

    [数理統計学]統計的推定のまとめ


    通勤電車のなかで私が勉強する用のシリーズ第4弾です。今回は統計的推定についてまとめておこうと思います。統計検定のための勉強のログのようなものです。概要と、簡単な例を可能な限り集めて載せていこうと思います。

    【これまでのシリーズへのリンク】


    目次

    不偏性

    • バイアスがゼロになるときに、その推定量は不偏性を持つという。不偏性のある推定量を不偏推定量(unbiased estimator)と呼ぶ。
    • バイアスは推定量の期待値と真のパラメータの値との差 \(E [\hat \theta] – \theta \)で表される。\( \hat \theta \)はサンプリングした分布の平均を意味しています。
    • 一意性はなく、多くの不偏推定量の中から良いものを選ぶ必要がある。任意の期待値0の確率変数を加えたとしても、不偏推定量になってしまうため、非常に弱い要請とされる。
    • 最尤推定量やモーメント推定量は一般には不偏性を持たない。
    • 一様最小分散不偏推定量(Uniformly Minimum-Variance Unbiased Estimator:UMVUE)は母数\(\theta\)の値によらず、一様に分散を最小化する不偏推定量のこと。
    • 標本平均、\( \bar X = \frac{1}{n} \sum_i X_i \)の期待値は\(E[\bar X] = \frac{1}{n}\sum_{i=1}^n E[X_i] = \frac{1}{n}n\mu = \mu \)となり、不偏推定量となる。
      標本分散、\( E \left[ \frac{1}{n} \sum_{i=1}^n (X_i – \bar X)^2 \right] = \frac{n-1}{n} \sigma^2 \)は\( \frac{\sigma^2}{n} \)だけ過少になるので、不偏推定量ではない。
      他方、\( E \left[ \frac{1}{n-1}\sum_{i=1}^n (X_i – \bar X)^2 \right] = \sigma^2 \)は不偏推定量となる。

    一致性

    • 標本サイズnの推定量\(\hat \theta_n\)がnの増大とともに\( \hat \theta_n \overset{p}{\to} \theta \)のように真のパラメータ\( \theta \)に確率収束することような推定量を一致推定量(consistent estimator)と呼ぶ。
    • 標本平均の場合、\( \hat X_n \overset{p}{\to} \mu \)となり、任意の正の\( \epsilon \)に対して\( \lim_{n \to \infty} P(| \hat X_n – \mu | > \epsilon ) = 0 \)が成り立つ。
    • 標本平均・標本分散などのモーメントの推定量は一致性を持つ。連続変形したモーメント推定量も一致性を持つ。
    • 最尤推定量は適当な正則条件のもと一致性を持ち、標本サイズが大きいとさらにより良い性質を持つ。
    • \( X_1, X_2, \dots , X_n \sim N(\mu, \sigma^2) \)として(ただし独立)、\( T_n (X_1, X_2, \dots, X_n) = \frac{1}{n} \sum_{i=1}^{n}(X_i – \bar X )^2 \)が\(\sigma^2\)の一致推定量であることを示す。
      大数の法則により、\( \bar X = \frac{1}{n} \sum_{i=1}^n X_i \)は\( \mu \)になることから、
      \( \frac{1}{n} \sum_{i=1}^n (X_i – \mu)^2 \\
      = \frac{1}{n}\sum_{i=1}^n X_i^2 – 2 \mu \times \frac{1}{n}\sum_{i=1}^n X_i + \mu^2 \\
      = \frac{1}{n}\sum_{i=1}^n ( (X_i – \mu ) + \mu )^2 – 2\mu \times \mu + \mu^2 \\
      = (\mu^2 + \sigma^2) – 2\mu \times \mu + \mu^2 = \sigma^2 \)
      となり、これは一致推定量となる。

    有効性

      ある不偏推定量\( \hat \theta \)が\( V_{\theta} [\hat \theta] = J_n (\theta)^{-1} \)のようにクラメール-ラオの下限を達成しているときに、有効推定量(efficient estimator)と呼ぶ。有効推定量は一様最小分散不偏推定量となる。

      ・正規分布において、標本平均は平均パラメータの有効推定量となる。
      ・バイアス補正を行った標本分散\( \frac{1}{n-1} \sum_{i=1}^n (X_i – \bar X)^2 \)は分散パラメータの一様最小分散不偏推定量だが、有効推定量ではない。

    推定量の相対効率

    • 2つの不偏推定量の分散の逆数の比のこと。
    • \( e( \hat \theta_1, \hat \theta_2 )= \frac{V_{\theta} [(\hat \theta_2)]}{V_{\theta} [(\hat \theta_1)]} \)
      などとおいて、\( e( \hat \theta_1, \hat \theta_2 ) > 1 \)のとき、\( \hat \theta_1 \)が優れている見なし、\( e( \hat \theta_1, \hat \theta_2 ) < 1 \)のとき、\( \hat \theta_2 \)が優れていると見なす。
    • \( e( \hat \theta_1, \hat \theta_2 )= \frac{E_{\theta}[(\hat \theta_2 – \theta )^2]}{E_{\theta}[(\hat \theta_1 – \theta)^2]} \)
      などとおいて、一般の推定量に拡張することもできる。

    最小分散不偏推定量

    クラメール・ラオの不等式

    • ある推定量が不偏推定量であるかどうかを確かめることができる条件式のこと。
    • 不偏推定量\( \hat \theta \)は適当な正則条件のもとで、
      \( V_{\theta}[\hat \theta] \geq J_n (\theta)^{-1} \)を満たすというもの。\( J_n(\theta) \)は対数尤度のパラメータに関する導関数の分散となっており、フィッシャー情報量と呼ばれる。

    フィッシャー情報量、フィッシャー情報行列

    • フィッシャー情報量は以下の式のことを指す。

      \( J_n(\theta) = E_{\theta} \left [ \left ( \frac{\partial}{\partial \theta} \log f(x;\theta) \right)^2 \right ] = V_{\theta} \left [ \left ( \frac{\partial}{\partial \theta} \log f(x;\theta) \right) \right ] \)

      データを観測したとき、確率密度関数の対数がパラメータが少し変わった際に、どれだけ変化するのかを測るもので、大きければ大きいほど推定量の分散の下限が小さくなる。パラメータがどれだけ効率的にもっともらしい値を捉えにいっているかを指し示していると思えばよい。行列の場合、クラメール-ラオの不等式は分散共分散行列からフィッシャー情報行列を差し引いたものが、非負値定符号であるかどうかを確かめるものとなる。
      \( Cov_{\theta} (T) \geq I^{-1}_{X_1, X_2, \dots, X_n}(\theta) \)
      これは同時に、以下の行列が非負値定符号であることを意味する。
      \( Cov_{\theta} (T) – I^{-1}_{X_1, X_2, \dots, X_n}(\theta) \)

      ここで、j行k列の成分は

      \( E_{\theta} \left [ \frac{ \partial }{ \partial \theta_j} \log \prod_{i=1}^n f(X_i , \theta) \frac{\partial}{\partial \theta_k} \log \prod_{i=1}^n f(X_i , \theta) \right ] \)

      となる。

    十分統計量

    • ある分布のパラメータ\( \theta \)を推定したいときに、分布から得られた標本Xのうち推定に十分な情報を含んだ統計量\( T = T(X) \)を十分統計量と呼ぶ。
    • \( P(X=x | T(x)=t, \theta) = P(X=x | T(x) = t) \)

      のように\( T(x) \)で条件づけたXの分布がパラメータによらないとき\( T(x) \)が十分統計量になる。

    • $$
      \begin{cases}
      X_i = 1 \ if \ ith \ coin \ is \ heads \\
      X_i = 0 \ if \ ith \ coin \ is \ tails
      \end{cases}
      $$
      ここでは\( T(X) = \sum_{i=1}^nX_i \)、つまり、コインの表の数の合計が十分統計量であることを示す。
      コイン投げの表の合計値は二項分布に従うので、表の出る割合\( \mu \)が与えられたもとでのXの条件付き確率は
      $$ P(X=x | \mu) = \prod_{i=1}^n \mu^{x_i}(1-\mu)^{1-x_i} \\
      = \mu^{T(x)}(1-\mu)^{n-T(x)} \\
      = (1 – \mu)^n \left ( \frac{\mu}{1-\mu} \right )^{T(x)}
      $$
      ここで\( T(X) \)に関する確率
      $$ P(T(X)=t | \mu) = \sum P(X=x|\mu) \\
      {}_n C _t (1 – \mu)^n \left ( \frac{\mu}{1-\mu} \right )^{T(x)}
      $$
      \( T(X) \)に関する条件付き確率は
      $$ P(X=x | T(X)=t, \mu) \\
      = \frac{P(X=x, T(X)=t | \mu)}{P(T(X)=t | \mu)} \\
      = \frac{(1 – \mu)^n \left ( \frac{\mu}{1-\mu} \right )^{t} 1(T(X)=t) }{{}_n C _t (1 – \mu)^n \left ( \frac{\mu}{1-\mu} \right )^{t} } \\
      = \frac{1}{{}_n C _t} 1(T(X)=t)
      $$
      となり、分布がパラメータによらず、十分統計量となっていることがわかる。( \( 1(T(X)=t) \) はインデックス関数)

    Newmanの因子分解定理

    • 十分統計量かどうかを簡単に判定するための方法で、\( T(X) \)が\( \theta \)の十分統計量であるとき、またそのときに限り、
      $$f(x; \theta) = h(x)g(T(X),\theta)$$
      となる関数hとgが存在するというもの。密度関数を\( \theta \)によらない関数と、\( \theta \)による関数の積に分解したとき、後者が\( T(X) \)のみを含むような分解が存在する。

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

    • 計算が容易ではない条件付き期待値\( E(S|T) \)の計算を行わずに、不偏推定量を求めることができる定理。
    • \( \delta (X) \)をパラメータ\( \theta \)のある推定量として、Tを\( \theta \)の十分統計量のうちの一つとする。十分統計量で条件付けた推定量の期待値を以下のように定義する。
      $$\delta_1 (T) = E_{\theta} [ \delta (X) | T ] $$
      その場合、ラオ・ブラックウェルの定理は以下を保証する。
      $$E_{\theta} [(\delta_1(T) – \theta )^2] \leq E_{\theta} [(\delta(X) – \theta )^2] $$
      つまり、\( \delta_1 \)の平均二乗誤差は\( \delta \)の平均二乗誤差以下になる。
    • ただし、十分統計量の関数になっている不偏推定量がいつも一様最小分散不偏推定量をもたらすとは限らず、それがラオ・ブラックウェルの定理の弱点とされている。
    • ケンブリッジ大学のラオ・ブラックウェルの定理に関する資料を見つけました。ポアソン分布と一様分布に関する例が載っています。

    順序統計量

    • 標本を昇順に並べ直したもの、もしくはそのうちの特定の順位の標本\( X_{(i)} \)に注目した統計量。
    • 中央値( \( X_{(\frac{n}{2})} \) )、最大値( \( X_{(n)} \) )、最小値( \( X_{(1)} \) )、四分位数などが順序統計量にあたる。
    • n個の標本を昇順に並べ替えた統計量\( (X_{(1)}, \dots, X_{(n)}) \)はフィッシャー・ネイマンの分解定理を用いることで統計量だけからなる関数に分解できることから、母集団の平均や分散の十分統計量になっている。
      $$f(x_1, \dots, x_n |\mu,\sigma) = \prod_{i=1}^n f(x_i | \mu, \sigma) = \prod_{j=1}^n f(x_{(j)}|\mu, \sigma) $$

      最後の項が並べ替えた統計量からなっている。

    最尤推定量

    • 尤度関数\( L(x, \theta) \)を最大化するパラメータ\( \hat \theta = \hat \theta (x) \)のことを最尤推定量と呼ぶ。
    • 最尤推定法は観測についての直感に基づいて定義された推定法とされる。推定の誤差についての記述は定義の中にない。

    尤度関数

    • 確率(密度)関数をxが固定された\(\theta\)の関数と考えたもの。
    • パラメータ\( \theta \)が持っている、観測値\( x_1, x_2, \dots , x_n \)を実現させるもっともらしさを表している。手元にある観測値は起こりにくい状態(パラメータ)から得られたと考えるより、起こるべくして起こったと考えるのが自然。

    対数尤度関数

    • 尤度関数の対数。\( l(\theta) = \log L(\theta) \)
    • 対数にすることで、独立な確率変数の同時確率を加算に直すことができたり、期待値の計算が簡単になったり、大数の法則や中心極限定理などを適用しやすくなる。

    有効スコア関数

    • 対数尤度がパラメータ\( \theta \)で偏微分できるとき、その微分値を\( \theta \)の関数としたもの。
      $$V(\theta) = V(x, \theta) = \frac{\partial}{\partial \theta} l(x, \theta) = \frac{\partial}{\partial \theta} \log f(x;\theta) $$
    • スコア関数の分散はフィッシャー情報量と呼ばれる。
    • 最尤推定値の定義から、任意のxに関してスコア関数は\( V(x, \hat \theta (x)) = 0 \)を満たす。(微分してゼロなので。)
    • \( E_{\theta}\left [ g(X) \right ] \)が\( \theta \)に依存しないような任意の関数\( g(x) \)について、
      $$ E_{\theta}\left [ g(X) V(X, \theta) \right ] \\
      = \int g(x)f(x;\theta)\frac{\partial}{\partial \theta} \log f(x;\theta)dx \\
      = \int g(x)f(x;\theta) \frac{ \frac{\partial}{\partial \theta} f(x;\theta) }{f(x;\theta)}dx \\
      = \frac{\partial}{\partial \theta} \int g(x)f(x ; \theta)dx \\
      = \frac{\partial}{\partial \theta} E_{\theta}\left [ g(X) \right ]=0 $$
      が成り立つ。これらの性質は情報量規準や漸近理論において役に立つ。

    モーメント法

    • モーメントを用いて、パラメータをモーメント関数で表し、それを推定値として用いる推定手法のこと。密度関数の形が複雑なときに、分布のパラメータ推定の手法として有力とされる。
    • m個のパラメータ\( \theta = (\theta_1, \dots, \theta_m) \)をもつ確率密度関数に関して、中心モーメントを以下のように定義する。
      $$ \mu_1 = \int x f(x;\theta)dx \\
      \mu_k = \int (x-\mu_1)^k f(x;\theta)dx \\
      k=2,3,\dots$$
      この場合、\( \mu_1 \)は平均、\( \mu_2 \)は分散となる。
      先ほどの中心モーメントを以下のように、パラメータの関数として表現することもできる。
      $$ \mu_1 = m_1(\theta) \\
      \mu_k = m_k(\theta)
      $$
      これらの式は\( \theta \)の連立方程式であるから、中心モーメントの関数としてパラメータを表すこともできる。
      $$ \theta_j = g_j (\mu_1, \dots, \mu_k) \\
      j = 1, \dots, m $$
      標本平均や標本分散を用いると、
      $$ \hat \theta_j = g_j (\hat \mu_1, \dots, \hat \mu_k) $$
      となり、実際はこれを推定値として使う。
    • 標本モーメントは真のモーメントに確率収束することから、\(g_j\)が連続関数のときに推定量の一致性が保証される。
    • ガンマ分布のパラメータ推定
      $$ f(x;\alpha, \lambda) = \frac{\lambda^\alpha}{\Gamma(\alpha)}x^{\alpha -1}e^{-\lambda x} \\
      (x \geq 0 ) $$
      期待値

      $$ \mu_1 = \frac{\alpha}{\lambda} \Rightarrow \alpha = \frac{\mu_1^2}{\mu_2} \Rightarrow \hat \alpha = \frac{\bar X^2}{\sum (X_i – \bar X)^2/n} $$

      分散

      $$ \mu_2 = \frac{\alpha}{\lambda^2} \Rightarrow \lambda = \frac{\mu_1}{\mu_2} \Rightarrow \hat \lambda = \frac{\bar X}{\sum (X_i – \bar X)^2/n} $$

    最小二乗法

    • 残差の二乗和\( \sum^n ( Y_i – f(X_i) )^2 \)を最小にするような関数fを選ぶ手法。二乗であると計算が簡単になる、あるいは、正規分布の撹乱項が加わっている\( Y_i = f(X_i) + \epsilon_i \)というモデルの最尤法と一致するなどの利便性がある。

    最良線形不偏推定量(BLUE)

    • 線形不偏推定量の中で平均二乗誤差\( E[(\hat \beta – \beta )^2] \)を最小化するような\( \hat \beta \)となる推定量のこと。
    • 各iについて、\( E(\epsilon_i) = 0, V(\epsilon_i)=\sigma^2 < \infty \)が共通で、さらに、\( i \neq j \)のとき、\( E(\epsilon_i \epsilon_j) = 0 \)とするときに、最小二乗推定量\( \hat \beta _{OLS} \)とBLUEは一致するという定理を、ガウス・マルコフの定理と呼ぶ。誤差項の期待値に関する仮定しかないという点でシンプルな定理となっている。
      $$ \hat \beta _{OLS} = \frac{\sum_i x_i Y_i }{ \sum_i x_i^2 } $$
      $$ E[\hat \beta _{OLS}] = E \left [ \frac{ \sum_i x_i(\beta x_i + \epsilon_i)}{\sum_i x_i^2} \right ] \\
      = \beta + \frac{\sum x_i E( \epsilon_i )}{\sum x_i^2} \\
      = \beta $$
    • 平均二乗誤差\( E_{\beta} [ (\hat \beta – \beta)^2 ] \)は未知である、真のパラメータ\( \beta \)を用いているので、実際に計算することはできないので、標本X,Yから計算できる二乗誤差を最小化することとなる。

    区間推定

    • 平均などのパラメータを推定する際に、ある一点ではなく区間でもって推定すること。
    • ある確率分布に従う標本\( X \sim P_{\theta} \)に対して、ある関数L, Uに対して
      $$P_{\theta}( L(X) < \theta < U(X)) \geq 1 – \alpha $$
      が全ての\( \theta \)で成立しているときに、\( (L(X),U(X)) \)を信頼係数\( 1 – \alpha \)の信頼区間と呼びます。これは真のパラメータ\(\theta\)が\( (L(X),U(X)) \)に入る確率ではなく、\( (L(X),U(X)) \)が真のパラメータ\(\theta\)を含む確率と考えます。被覆確率とも呼ばれています。
    • 先ほど、信頼係数というものが出てきましたが、真の値が入る確率ではありません。信頼係数が95%あるということは、10000回標本を集めて、パラメータを推定し区間を求めた際に、500回は母集団平均を含む区間を得られない(9500回は母集団平均を含む範囲を得られる)ということを意味しています。ただ、実際は研究する際に、「標本を集めることはたかだか1回しかできず、その1回の推定が20回に1回は母集団平均を含んでくれないという前提」を知った上で、結果を甘んじて受け入れるだけです。
    • 正規分布に従う標本平均\( \bar X_n \sim N( \mu, \frac{1}{n} ) \)について考える。片側\( \frac{\alpha}{2} \)の点をおくと、以下のようになる。

      $$ – \frac{Z_{\frac{\alpha}{2}}}{\sqrt n } < \bar X – \mu < \frac{Z_{\frac{\alpha}{2}}}{\sqrt n } $$

      これは確率にすると\( \alpha \)を差っ引いた、\( 1 – \alpha \)がこの範囲の取りうる確率となる。これを書き換えると、以下のようになり、信頼区間の上限と下限を計算できる。
      $$ P_{\mu} \left ( \bar X – \frac{Z_{\frac{\alpha}{2}}}{\sqrt n } < \mu < \bar X \frac{Z_{\frac{\alpha}{2}}}{\sqrt n } \right ) = 1 – \alpha $$

    相関係数の区間推定

    • 相関係数の分布は歪んだ形をしているため、正規分布で近似することに問題がある。そこで、Fisherのz変換と呼ばれる手法を用いることで相関係数の分布が正規分布に従うようにする。
    • Fisherのz変換を行った場合、nが十分大きいときに近似的に相関係数は以下の正規分布に従う。
      \( \xi (r) = \frac{1}{2} \log \frac{1+r}{1-r} \)
      平均:\( \xi (\rho_{XY} ) = \frac{1}{2} \log \frac{1+\rho_{XY}}{1-\rho_{XY}} \)
      分散:\( \frac{1}{n-3} \)

      片側\( \frac{\alpha}{2} \)の点をおくと、以下のようになる。
      $$ – \frac{Z_{\frac{\alpha}{2}}}{\sqrt{n-3} } < \xi (r) – \xi (\rho_{XY}) < \frac{Z_{\frac{\alpha}{2}}}{ \sqrt{n-3} } $$
      これは確率にすると𝛼を差っ引いた、1–𝛼がこの範囲の取りうる確率となる。

      $$ P \left ( \xi(r) – \frac{Z_{\frac{\alpha}{2}}}{\sqrt{n-3} } < \xi (\rho_{XY}) < \xi (r) + \frac{Z_{\frac{\alpha}{2}}}{\sqrt{n-3} } \right ) \simeq 1 – \alpha $$

      これを書き換えると、以下のようになり、信頼区間の上限と下限を計算できる。

      $$ \xi^{-1} \left ( \xi(r) – \frac{Z_{\frac{\alpha}{2}}}{\sqrt{n-3} } \right ) < \rho_{XY} < \xi^{-1} \left ( \xi (r) + \frac{Z_{\frac{\alpha}{2}}}{\sqrt{n-3} } \right ) $$

      ここで逆関数は\( \xi^{-1}(x) = \frac{e^x – e^{-x}}{e^x + e^{-x}} \)となっているものとする。

    デルタ法

    • 推定量の漸近正規性を示すことができる方法の一つ。ある変数が正規分布に従う場合、その連続関数も漸近的に正規分布に従うという理論。GreenのAppendixに載っている記述によると、
      \( \sqrt{n}(z_n – \mu ) \overset{d}{\to} N(0, \sigma^2) \)かつ\( g(z_n) \)がnと関係しない連続関数であれば、
      $$ \sqrt{n}[ g(z_n) – g(\mu) ] \overset{d}{\to} N(0, \{g'(\mu)\}^2\sigma^2) $$
      が成り立つと記されています。これはよく見ると、線形のテイラー近似における平均と分散に他なりません。
      $$ g(z_n) \simeq g(\mu) + g'(\mu)(z_n – \mu) $$

    AIC

    • 赤池情報量規準。多項式の中で何次の多項式を用い、さらにどのような係数パラメータを用いれば良い推定を行えるのか(モデル選択の問題)を追究することが目的の指標。対数尤度などからなり、それに負号をかけているため小さければ小さいほど良い指標となる。
      AICは以下の式で表される。

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

      第1項は対数尤度、第2項は罰則項。2倍されているのは、モデルが正しいときに、第1項が漸近的にカイ二乗分布に従うため。より小さい値を求めているため、対数尤度は大きい方がいいし、罰則項は小さい方がいい。

    • AICで比較できるモデルは、それぞれが包含関係にあるようなもののみとされているので扱う際は注意が必要。
    • 通っていた大学の教授が、赤池先生は生きていればノーベル経済学賞を受賞できたのではないかと仰っていたのを思い出した。
    • 多項式回帰の場合のAIC
      まず、対数尤度は以下の形を考える。
      $$ \log f(y ; \hat \alpha, \hat \sigma) \\
      = \frac{n}{2} \log (2\pi \hat \sigma^2) – \sum_{i=1}^n (Y_i – f(X_i ; \hat \alpha))^2 \\
      = -\frac{n}{2} \log (2\pi) – \frac{n}{2} \log \hat \sigma^2 – \frac{n}{2} $$
      この場合、AICは以下のようになる。

      $$ AIC = n \log 2\pi + n + n \log \hat \sigma^2 + 2(d+2) $$

      最後の項は次数dの多項式のd、切片、標本分散の数からd+2となる。

    カルバック-ライブラー情報量とAICの導出

    • カルバック=ライブラー情報量(KLダイバージェンス)は二つの分布(例えばF,Gとする)の間の近さを測るための指標で、それぞれの密度関数をf,gとすると、
      $$KL( f \| g ) = \int f(x) \log \frac{f(x)}{g(x)} dx $$
      で定義される。ただし、\( – 0 \log 0 = 0 \)となる。カルバック=ライブラー情報量には対称性はなく数学的な距離を測るものではない。
      カルバック=ライブラー情報量の下限が0であることを示すために、-logの凸性よりイェンセンの不等式を用いると、
      $$KL( f \| g ) = \int f(x) \left ( – \log \frac{g(x)}{f(x)} \right ) dx \\
      \geq – \log \left ( \int f(x) \frac{g(x)}{f(x)}dx \right ) \\
      = – \log 1 = 0 $$
      となる。また、カルバック=ライブラー情報量の最大化は平均対数尤度\( E_f [\log g(X)] \)の最大化と対応している。
    • カルバック=ライブラー情報量はAICの導出において重要な役目を果たす。
      カルバック=ライブラー情報量で真の分布と真の分布に近づけようとする推定量\( \hat \theta \)をもった分布との距離を計算すると、

      $$ KL(f(y;\theta) \| f(y; \hat \theta (X)) )\\
      = \int f(y;\theta)\log f(y ; \theta)dy – \int f(y;\theta) \log f(y;\hat \theta(x))dy $$

      となる。第1項目は推定量\( \hat \theta \)とは無関係に定まるが、第2項には含まれている。カルバック=ライブラー情報量の下限が0であることから、
      $$\int f(y;\theta) \log f(y;\hat \theta(x))dy$$
      を最大化するような\( \hat \theta \)が良い推定量となる。
      実際に計算できる最大対数尤度は、
      $$ \sum_{i=1}^n \log f(X_i ; \hat \theta_n^{ML}(X) ) $$
      であり、ここで、真のパラメータ\( \theta \)を持つ分布から標本\( Y_1, \dots, Y_n \)が無作為抽出で得られたとして、XとYが同じ分布に従うことから、

      $$ \sum_{i=1}^n \log f(X_i ; \hat \theta_n^{ML}(X) ) \sim \sum_{i=1}^n \log f(Y_i ; \hat \theta_n^{ML}(Y) ) $$

      となる。しかしながら、先ほどカルバック=ライブラー情報量にあった推定値は、
      $$\sum_{i=1}^n \log f(Y_i ; \hat \theta_n^{ML}(X) )$$
      であるから、これは最適ではないため以下の不等号となる。

      $$ \sum_{i=1}^n \log f(X_i ; \hat \theta_n^{ML}(X) ) \sim \sum_{i=1}^n \log f(Y_i ; \hat \theta_n^{ML}(Y) ) \\
      \geq \sum_{i=1}^n \log f(Y_i ; \hat \theta_n^{ML}(X) ) $$

      となる。AICの第1項目に現れる最大対数尤度は、真のパラメータに従う分布から得られたYをもとに計算した平均対数尤度よりも大きくなってしまうため、つまり、負号をつけると、不当に小さくなってしまうため、補正するための項が必要になる。それがAICの第2項目となる。
      実際のところ、AICの第2項目は漸近理論において、標本サイズnが十分に大きいときに、補正のための項がパラメータの次元\( dim(\theta) \)で近似されるところに由来している。

    参考文献

    [1]
    [2]日本統計学会(2013), 『日本統計学会公式認定 統計検定1級対応 統計学』, 東京図書
    [3]鈴木・山田(1996),『数理統計学―基礎から学ぶデータ解析』,内田老鶴圃
    [4]William H. Greene (2019),”Econometric Analysis, Global Edition”, Pearson Education

    BLSTMを用いた文書分類でデザイナーズマンション予測に再挑戦

    はじめに

    仕事で深層学習を使うことはないのですが、扱える技術の幅を広げていきたいなと思い、BLSTMを用いた文書分類についてkerasでの簡単なサンプルコードをやってみようと思います。データは以前集めた、某不動産紹介サイトの賃貸マンションの設備に関する文書とそれがデザイナーズマンションかどうかのラベルを使います。そして文書内の単語からその文書がデザイナーズマンションかどうかを予測します。前回はAUCで83%だったので、それを超えれると良いですね。

    目次

    単純なRNNとは

    • モチベーション
      • フィードフォワード型のニューラルネットワークではうまく扱うことができない時系列データをうまく扱えるようにすること。
    •  特徴
      • 入力が互いに関係している(多層パーセプトロンの際は置かれていない仮定)
        • 直訳すると循環するニューラルネットワークとなる。
      • 最初の文の単語が2つ目、3つ目の単語に影響を与える可能性を考慮。
      具体的な関数の形としては、 $$ h_t = \tanh (W h_{t-1} + U x_t) \\\ y_t = softmax(V h_t) $$ で与えられる。
      \( h_t \)は隠れ層の状態を、\( x_t \)は入力変数を、\( y_t \)は出力ベクトルを表している。(他のアルファベットは重み行列)

    LSTMとは

    • モチベーション
      • 単純なRNNの勾配消失問題を解決するために提案された手法。
    • 特徴
      • 単純なRNNに置き換えることで性能が大幅に向上することも珍しくない。
        • 時系列データ、長い文章、録音データからなる長期的なパターンを取り出すことを得意としている手法。
      • 勾配消失問題に強い
      • Long Short-Term Memory:長短期記憶
        • 長期依存性を学習できるRNNの亜種。
        • 入力ゲート、忘却ゲート、出力ゲート、内部隠れ層という4つの層が相互に関わり合う。重要な入力を認識し、それを長期状態に格納すること、必要な限りで記憶を保持すること、必要なときに記憶を取り出すことを学習する。
          • 入力ゲート:内部隠れ層のどの部分を長期状態に加えるかを決める。
          • 忘却ゲート:長期状態のどの部分を消去するか決める。
          • 出力ゲート:各タイムステップで、長期状態のどの部分を読み出し、出力するかを決める。
    $$
    i = \sigma (W_i h_{t-1} + U_i x_t) \\
    f = \sigma (W_f h_{t-1} + U_f x_t) \\
    o = \sigma (W_o h_{t-1} + U_ox_t ) \\
    g = \tanh (W_g h_{t-1} + U_g x_t) \\
    c_t = (c_{t-1} \otimes f ) \otimes ( g \otimes i) \\
    h_t = \tanh (c_t) \otimes o
    $$ i:入力ゲート(input)
    f:忘却ゲート(forget)
    o:出力ゲート(output)
    \( \sigma \):シグモイド関数
    g:内部隠れ層状態
    \(c_t \):時刻tにおけるセル状態
    \(h_t \):時刻tにおける隠れ状態

    Bidirectional LSTMとは

    • モチベーション
      •  従来のRNNの制約を緩和するために導入。
    • 特徴
      • ある特定の時点で過去と将来の利用可能な入力情報を用いて訓練するネットワーク
        • 従来のRNNのニューロンをフォワード(未来)なものとバックワード(過去)なものとの2つに分ける。
          • 2つのLSTMを訓練している。
        • 全体のコンテキストを利用して推定することができるのが良いらしい。
          • 文章で言うと、文章の前(過去)と後(未来)を考慮して分類などを行うイメージ。
        • 人間もコンテキストを先読みしながら解釈することもある。
      • BLSTMは全てのシークエンスの予測問題において有効ではないが、適切に扱えるドメインで良い結果を残す。
      • 1997年とかなり歴史があるものらしい
    出典:Deep Dive into Bidirectional LSTM

    Bidirectional LSTMで文書分類

    今回は、BLSTMを使って、マンションの設備に関するテキスト情報から、そのマンションがデザイナーズマンションかどうかを予測します。全体のソースコードはGoogle Colabを御覧ください。

    データ

    データを確認してみると、
     
    今回のデータがテキストとラベルからなっていることがわかる。なお、データ数は1864件あり、そのうち16%がデザイナーズマンションです。

    前処理

    Google Colab上で形態素解析を行うために、MeCabをインストールする。
     
    これでMaCab NeologdをGoogle ColabのPythonで実行できます。
    テキストデータの前処理を行うために名詞だけを抽出してリスト形式で返す関数を定義します。
     
    ここで、語彙に関するデータを作成します。
     
    続いて、単語とそのインデックスからなるdict形式のデータを作成します。
     
    最後に、これまでに作成した語彙とそのインデックスのデータを用いて、実際に学習に使う入力データと出力データを作成します。
     

    実行

    先程作成したデータを訓練データとテストデータに分けます。
    ここで、AUCを計算するための関数を定義します。(まるまる拝借しました)
    続いて、kerasを用いて、ネットワークを構築し、コンパイルを行います。
      テストデータでの予測精度を確認します。
    実際のラベルとテキストデータを見比べてみます。
    LSTMの場合は、BLSTMに関する記述(1行)をネットワークから除外すれば良いので簡単に試すことができます。 せっかくなので比較してみたところ、AUCは0.841でした。BLSTMは0.844だったので、若干上回る程度のようです。 以前扱った記事ではデザイナーズマンション分類はAUCが0.83程度だったので、LSTMを使ったほうが精度が少しだけ高いようです。

    追記

    TensorBoardの検証のloglossを見てみます。 エポックに対してloglossが非常に荒れています。学習曲線としてはセオリー的にアウトなようです。一方で、検証のAUCはエポックに従って高まるようです。 AUCは良くなっているけどloglossが増え続けている。loglossが下がらないと過学習している可能性が高いので、これは過学習しているだけなのだろう。

    参考文献

     

    [数理統計学]正規分布から導かれる分布(カイ二乗分布/t分布/F分布)の期待値と分散の導出まとめ

    はじめに

    前回は、連続型確率分布に関するよくある分布とその平均と分散の導出についてひたすら記しました。今回は同じく連続型ながら、正規分布から導かれる確率分布に関して同様に記していこうと思います。大学のときに数理統計学の授業でがっつりやったはずの領域だったのですが、見事に忘れていたので、電車の中などでしっかりと復習していこうと思います。

    ※PC版でないと数式が切れてしまうので、SP版での閲覧はおすすめしません。私、SKUEのスマホだけは表示崩れがないようにチェックしております。通勤電車で定期的に読むので。

    今回登場する連続分布

    『数理統計学―基礎から学ぶデータ解析』という本に登場するものを式を補いながら紹介します。

    1.カイ二乗分布

    ガンマ分布\( \Gamma \left ( \frac{n}{2}, \frac{1}{2} \right ) \)、すなわち、密度関数が以下の式で与えられるとき、この確率分布を自由度nの\( \chi^2 \)分布という。

    $$ f(x ; n) = \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}} \\
    x > 0 \\
    n = 1, 2, \dots
    $$

    \( \chi^2 \)分布の再生性

    XとYがカイ二乗分布に従うとして、それらの和がカイ二乗分布になる性質を再生性と呼ぶ。
    $$ X \sim \chi^2 (m) \\
    Y \sim \chi^2 (n) \\
    X + Y \sim \chi^2 (m+n)
    $$

    以下ではこの再生性に関して証明を行う。
    ここで\( z = X + Y \)とする。zの密度関数は、以下で表される。

    $$ g(z) = \int_0^z f(x ; m) f(y ; n) dx \\
    = \int_0^z f(x ; m) f(z-x ; n) dx \\
    = \int_0^z \frac{ x^{\frac{m}{2}-1} e^{-\frac{x}{2}} }{2^{\frac{m}{2}} \Gamma \left ( \frac{m}{2} \right ) } \frac{ (z-x)^{\frac{n}{2}-1} e^{-\frac{z-x}{2}} }{2^{\frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) } dx
    $$

    ここで、\( u = \frac{x}{z} \)として、変数変換を行う。( \( \frac{dx}{du} = z \quad , x = zu \) )

    $$ = \int_0^1 \frac{ (zu)^{\frac{m}{2}-1} e^{-\frac{zu}{2}} }{2^{\frac{m}{2}} \Gamma \left ( \frac{m}{2} \right ) } \frac{ (z-zu)^{\frac{n}{2}-1} e^{-\frac{z(1-u)}{2}} }{2^{\frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) } z du $$

    $$ = \frac{z^{\frac{m}{2} + \frac{n}{2} -1} e^{-\frac{z}{2}} }{ 2^{\frac{m}{2} + \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \int_0^1 u^{\frac{u}{2}-1}(1-u)^{\frac{u}{2}-1}du $$

    積分記号以降の部分はベータ分布であるから、

    $$ = \frac{ e^{-\frac{z}{2}} z^{\frac{m+n}{2} -1} }{ 2^{\frac{m+n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } B \left ( \frac{m}{2}, \frac{n}{2} \right ) $$

    以前記したベータ分布とガンマ関数との関係( \( B(x, y) = \frac{\Gamma(x) \Gamma(y)}{\Gamma(x+y)} \) )より、

    $$ = \frac{ e^{-\frac{z}{2}} z^{\frac{m+n}{2} -1} }{ \Gamma \left ( \frac{m+n}{2} \right ) 2^{\frac{m+n}{2} } } $$

    これは\( \chi^2 (m+n) \)の確率密度関数である。

    定理:\( X_1, X_2, \dots , X_n \)は独立に正規分布\( N(0, 1) \)に従うとき、\( Y = X_1^2 + X_2^2 + \dots + X_n^2 \)は\( \chi^2(n) \)に従う。

    以下はUCLAのレクチャーノートに従っています。

    まず、\( X_i^2 \)がカイ二乗分布に従うことを示します。
    標準正規分布に従う変数zの二乗としてxを定義します。
    $$ x = z^2 \\
    f(z) = \frac{1}{\sqrt {2 \pi } } e^{-\frac{1}{2}z^2}
    $$
    xの確率は以下のように表される。
    $$ F_X (x) = P( X \le x ) \\
    = P(Z^2 \le x ) \\
    = P( – \sqrt {x} \le Z \le \sqrt {x} )
    $$

    この場合、求めたい確率は1から両端を差し引いたものとなる。

    $$ {left} = F_z(-\sqrt {x}) \\
    {right} = 1 – F_z(\sqrt {x}) \\
    F_X (x) = 1 – {left} – {right} \\
    = F_z(\sqrt {x}) – F_z(-\sqrt {x})
    $$

    確率密度を求めるために両辺をxで微分し整理する。

    $$ \frac{d F_X (x) }{dx} = \frac{d F_z(\sqrt {x}) }{ d \sqrt {x} } \frac{d \sqrt {x}}{dx} – \frac{d F_z( – \sqrt {x}) }{ d – \sqrt {x} } \frac{d – \sqrt {x}}{dx} \\\\
    = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x}\frac{1}{2}x^{-\frac{1}{2}} – (-1) \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x}\frac{1}{2}x^{-\frac{1}{2}} \\\\
    = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x}x^{-\frac{1}{2}} \\\\
    = \frac{x^{-\frac{1}{2}} e^{-\frac{x}{2}}}{2^{\frac{1}{2}} \Gamma \left ( \frac{1}{2} \right ) } \sim \chi_1^2
    $$

    これは自由度1のカイ二乗分布であり、正規分布に従う変数の二乗は自由度1のカイ二乗分布に従うことが示された。

    ここで、ガンマ分布の積率母関数について導出を行う。

    $$ M_{X}(t) = E \left ( e^{tX} \right ) \\
    = \int _{0}^{\infty} e^{tx} \frac{x^{\alpha – 1} e^{-\frac{x}{\beta}}}{\beta^{\alpha} \Gamma \left ( \alpha \right )} dx \\
    = \frac{1}{\beta^{\alpha} \Gamma \left ( \alpha \right )}\int _{0}^{\infty} x^{\alpha – 1} e^{-x \left ( \frac{1-\beta t}{\beta} \right ) }dx
    $$

    ここで以下の変数変換を行う。
    $$ y = x \left ( \frac{1-\beta t}{\beta} \right ) $$
    つまり、\( x = \frac{\beta}{1-\beta t} y \)、\( \frac{dx}{dy} = \frac{\beta}{1-\beta t} \)であるから、これらを代入すると、

    $$ M_X(t) = \frac{1}{\beta^{\alpha} \Gamma \left ( \alpha \right )} \int _{0}^{\infty} \left ( \frac{\beta}{1-\beta t} \right )^{\alpha – 1}y^{\alpha – 1}e^{-y}\frac{\beta}{1-\beta t}dy \\\\
    = \frac{1}{\beta^{\alpha} \Gamma \left ( \alpha \right )} \left ( \frac{\beta}{1-\beta t} \right )^{\alpha – 1} \frac{\beta}{1-\beta t} \int _{0}^{\infty} y^{\alpha – 1}e^{-y} dy
    \\\\
    = (1-\beta t)^{-\alpha}
    $$

    よって、先程求めた自由度1のカイ二乗分布は積率母関数で表すと以下のようになる。
    $$ M_X(t) = M_{Z^2}(t) = (1-2t)^{-\frac{1}{2}} $$

    今回の定理において、\( X_i^2 \)は各々が独立なので、先程求めた自由度1のカイ二乗分布の積率母関数を掛け合わせることで\( Y = \sum_{i}^{n} X_i^2 \)の積率母関数は求まる。

    $$ M_Y(t) = M_{Z_1^2}(t) \times M_{Z_2^2}(t) \times \dots \times M_{Z_n^2}(t) \\\\
    = (1-2t)^{-\frac{1}{2}} \times (1-2t)^{-\frac{1}{2}} \times \dots \times (1-2t)^{-\frac{1}{2}} \\\\
    = (1-2t)^{-\frac{n}{2}}
    $$

    これは自由度nのカイ二乗分布の積率母関数と同じであるから、\( Y = \sum_{i}^{n} X_i^2 \)が自由度nのカイ二乗分布に従うことが示された。

    カイ二乗分布の平均

    $$ E \left ( X \right ) = \int _{0}^{\infty} x \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\
    = \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}} e^{-\frac{x}{2}}dx
    $$

    ここで、部分積分の公式より

    $$ \int _{0}^{\infty} x^{\frac{n}{2}} e^{-\frac{x}{2}}dx \\\\
    = \left [ \frac{2}{n+2}x^{\frac{n}{2}+1} e^{-\frac{x}{2}} \right ]_0^{\infty}- \frac{n}{2}(-2) \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx
    \\\\
    = n \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx
    $$

    となるから、カイ二乗分布の積分が1になることから、平均は

    $$ E \left ( X \right ) = n \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\
    = n
    $$
    となる。

    カイ二乗分布の分散

    分散は
    $$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 $$
    より求まるので、

    $$ E(X^2) = \int _{0}^{\infty} x^2 \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\
    = \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}+1} e^{-\frac{x}{2}}dx
    $$

    ここで、部分積分の公式より

    $$ \int _{0}^{\infty} x^{\frac{n}{2}+1} e^{-\frac{x}{2}}dx \\\\
    = \left [ \frac{2}{n+4}x^{\frac{n}{2}+2} e^{-\frac{x}{2}} \right ]_0^{\infty}- \frac{n+2}{2}(-2) \int _{0}^{\infty} x^{\frac{n}{2}} e^{-\frac{x}{2}}dx \\\\
    = (n+2)\int _{0}^{\infty} x^{\frac{n}{2}} e^{-\frac{x}{2}}dx \\\\
    = (n+2) \left [ \left [ \frac{2}{n+2}x^{\frac{n}{2}+1} e^{-\frac{x}{2}} \right ]_0^{\infty}- \frac{n}{2}(-2) \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \right ] \\\\
    = n(n+2) \int _{0}^{\infty} x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx
    $$

    となるから、カイ二乗分布の積分が1になることから、

    $$ E(X^2) = n(n+2) \int _{0}^{\infty} \frac{1}{ 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right ) }x^{\frac{n}{2}-1} e^{-\frac{x}{2}}dx \\\\
    = n(n+2)
    $$

    よってカイ二乗分布の分散は
    $$ V(X) = n(n+2) – n^2 \\
    = 2n $$
    となる。

    2.t分布

    以下の確率密度関数で与えられる分布をt分布と呼ぶ。

    $$ f(x ; n) = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } \\
    -\infty < x < \infty \\
    n = 1, 2, \dots
    $$

    確率密度関数から明らかなように、パラメータは自由度nだけからなる。
    なお、自由度が1のとき、t分布はコーシー分布と同じになる。
    $$ f(x ; 1) = \frac{1}{\pi} \frac{1}{1 + x^2} $$

    t分布の特徴としては、
    ・y軸に関して対称
    ・標準正規分布よりも左右に裾が重い
    ・自由度が1よりも大きいとき、xf(x)は可積分となり、平均は0となる。
    ・自由度が1の場合はコーシー分布と同じであるため平均も分散も存在しない。
    などがあげられる。

    定理:Xは\( N(0, 1) \)に、Yは\( \chi^2(n) \)にそれぞれ従い、かつ互いに独立とする。このとき、\( t = \frac{x}{\sqrt {y/n}} \)は\( t(n) \)に従う。

    XとYの結合密度関数は
    $$ f(x, y) = \frac{1}{\sqrt {2\pi}} e^{ -\frac{x^2}{2}} \frac{ y^{ \frac{n}{2} -1} }{ 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) } e^{ – \frac{y}{2} } $$
    となる。
    $$ t = \frac{x}{\sqrt {y / n} } \\
    x = t \sqrt {y / n} = t \sqrt {s / n} \\
    y = s
    $$
    とおいて、確率の変数変換を行うと、tとsの結合密度関数は、

    $$ g(t, s) = f \left ( t \sqrt {s / n}, s \right ) | J | \\
    = f \left ( t \sqrt {s / n}, s \right ) \begin{bmatrix} \frac{\partial x}{\partial t} & \frac{\partial x}{\partial s} \\ \frac{\partial y}{\partial t} & \frac{\partial y}{\partial s} \end{bmatrix} \\
    = f \left ( t \sqrt {s / n}, s \right ) \sqrt {s / n} \\
    = \frac{1}{\sqrt {2\pi}} e^{ -\frac{t^2 \frac{s}{n} }{2}} \frac{ s^{ \frac{n}{2} -1} }{ 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) } e^{ – \frac{s}{2} } \sqrt {s / n} \\
    = \frac{ s^{ \frac{1}{2}(n-1)} e^{-\frac{1}{2} \left ( 1 + \frac{t^2}{n}s \right ) } }{ \sqrt {2\pi} 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) \sqrt {n} }
    $$

    ここで\( u = \frac{1}{2} \left ( 1 + \frac{t^2}{n} \right ) s \)とおく。なお、\( s = \frac{2u}{ \left ( 1 + \frac{t^2}{n} \right ) } \)、\( \frac{ds}{du} = \frac{2}{1 + \frac{t^2}{n}} \)であるからこれらを代入し、分子について定積分をとりtを固定し、sを変換する。

    $$ \int _{0}^{\infty} s^{ \frac{1}{2} (n-1)} e^{ – \frac{1}{2} \left ( 1 + \frac{t^2}{n} \right )s }ds \\
    = \int _{0}^{\infty} \left ( \frac{2u}{ 1 + \frac{t^2}{n} } \right )^{ \frac{1}{2}(n-1) }e^{-u} \frac{ds}{du}du \\
    = \int _{0}^{\infty} \left ( \frac{2u}{ 1 + \frac{t^2}{n} } \right )^{ \frac{1}{2}(n-1) }e^{-u} \frac{2}{1 + \frac{t^2}{n}} du \\
    = \frac{ 2^{ \frac{1}{2} (n-1) + 1 } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n-1) + 1 } } \int _{0}^{\infty} u^{ \frac{1}{2}(n-1)} e^{-u}du \\
    = \frac{ 2^{ \frac{1}{2} (n+1) } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n+1)} } \int _{0}^{\infty} u^{ \frac{1}{2}(n-1)} e^{-u}du
    $$

    ここで、ガンマ関数の定義より、
    $$ \Gamma ( \alpha ) = \int_{0}^{\infty} u^{ \alpha – 1} e^{-u}du \\
    \alpha – 1 = \frac{1}{2}(n-1) \\
    \alpha = \frac{1}{2}(n+1)
    $$
    であるから、これを考慮すると、
    $$ = \frac{ 2^{ \frac{1}{2} (n+1) } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n+1)} } \Gamma \left ( \frac{1}{2} (n+1) \right ) $$
    となる。

    以上より、tの密度関数は

    $$ g(t) = \int_{0}^{\infty} g(t, s) ds \\\\
    = \frac{ 2^{ \frac{1}{2} (n+1) } }{ \left ( 1 + \frac{t^2}{n} \right )^{ \frac{1}{2}(n+1)} } \Gamma \left ( \frac{1}{2} (n+1) \right ) \frac{1}{\sqrt {2\pi} 2^{ \frac{n}{2} } \Gamma \left ( \frac{n}{2} \right ) \sqrt {n}} \\\\
    = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt{n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{t^2}{n} \right )^{ \frac{n+1}{2}} }
    $$

    これはt分布の密度関数と同じである。

    t分布の平均

    $$ E \left ( X \right ) = \int _{-\infty}^{\infty} x \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx $$
    奇関数であることから
    $$ = 0 $$
    となる。

    t分布の分散

    分散は
    $$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 \\
    = E \left ( X^2 \right ) $$
    より求まる。

    $$ E \left ( X^2 \right ) = \int _{-\infty}^{\infty} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx \\\\
    = \int _{-\infty}^{0} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx + \int _{0}^{\infty} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx
    $$

    偶関数であることから以下のように表される。
    $$ = 2 \int _{0}^{\infty} x^2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right ) \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} } dx \\
    = 2 \frac{ \Gamma \left ( \frac{n+1}{2} \right ) }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \int _{0}^{\infty} \frac{x^2}{ \left ( 1 + \frac{x^2}{n} \right )^{\frac{n+1}{2}} }dx
    $$

    ここで、\( y = \frac{1}{1+\frac{x^2}{n}} \)により変数変換を行う。
    なお、\( \left ( 1 + \frac{x^2}{n} \right ) = \frac{1}{y} \)、\( x = \left [ \frac{(1-y)n}{y} \right ] ^{\frac{1}{2}} \)、
    \( \frac{dx}{dy} = \frac{ (1-y)^{-\frac{1}{2}} n^{\frac{1}{2}}(-1)y^{-\frac{3}{2}}}{2} \)となる。

    $$ = (-1)2\frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \int _{0}^{\infty} \frac{(1-y)^{\frac{1}{2}} y^{\frac{n+1}{2}-1-\frac{3}{2}} }{2}dy \\\\
    = 2\frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \int _{0}^{\infty} \frac{(1-y)^{\frac{1}{2}} y^{\frac{n}{2}-2} }{2}dy
    $$

    ベータ関数( \( B( \alpha, \beta ) = \int_0^{1} t^{\alpha – 1}(1-t)^{ \beta – 1} dt \) )より、

    $$ = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} B\left ( \frac{n-2}{2}, \frac{3}{2} \right ) $$

    ベータ関数とガンマ関数の関係より、
    $$ B\left ( \frac{n-2}{2}, \frac{3}{2} \right ) = \frac{ \Gamma \left ( \frac{n-2}{2} \right ) \Gamma \left ( \frac{3}{2} \right ) }{ \Gamma \left ( \frac{n-2}{2} + \frac{3}{2} \right ) } \frac{1}{2} \\
    = \frac{ \Gamma \left ( \frac{n-2}{2} \right ) \Gamma \left ( \frac{3}{2} \right ) }{ \Gamma \left ( \frac{n}{2} + \frac{1}{2} \right ) }
    $$

    よって、
    $$ = \frac{ \Gamma \left ( \frac{n+1}{2} \right ) n^{\frac{3}{2}} }{ \sqrt {n \pi} \Gamma \left ( \frac{n}{2} \right )} \frac{ \Gamma \left ( \frac{n-2}{2} \right ) \Gamma \left ( \frac{3}{2} \right ) }{ \Gamma \left ( \frac{n}{2} + \frac{1}{2} \right ) } $$

    ガンマ関数は\( \Gamma (n) = (n-1) \Gamma (n-1) \)となることから、\( \Gamma (n+1) = n \Gamma (n) \)あるいは\( \Gamma \left( \frac{n}{2} – 1 \right) = \frac{ \Gamma \left ( \frac{n}{2}
    \right ) }{\frac{n}{2} – 1} \)であることを用いると、
    $$ = \frac{n^{ \frac{3}{2}} \Gamma \left ( \frac{n}{2} -1 \right ) \Gamma \left ( \frac{n}{2} + 1 \right ) }{ n^{\frac{1}{2}} \sqrt{\pi} \Gamma \left ( \frac{n}{2} \right ) } \\
    = \frac{n \Gamma \left ( \frac{n}{2} \right ) }{ \pi^{\frac{1}{2}} \Gamma \left ( \frac{n}{2} \right ) \left ( \frac{n}{2} – 1 \right )} \frac{1}{2} \Gamma \left ( \frac{n}{2} \right ) \\
    = \frac{n}{n-2}
    $$
    となる。

    よって、
    $$ V(X) = \frac{n}{n-2} $$

    3.F分布

    以下の確率密度関数で与えられる分布をF分布と呼ぶ。
    $$ f(x;m;n) = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } \\
    x > 0 \\
    m, n = 1, 2, \dots
    $$

    定理:Xは\( \chi^2(m) \)に、Yは\( \chi^2(n) \)に独立に従うとする。このとき、\( F = \frac{ \frac{1}{m}X }{ \frac{1}{n}Y } \)は\( F(m, n) \)に従う。

    ここで、\( (X, Y) \rightarrow (F, U) \)と確率の変数変換を行う。

    \(
    \begin{cases}
    F = \frac{ \frac{1}{m}X }{ \frac{1}{n}Y } \\
    U = Y \\
    \end{cases}
    \)
    また、\( X = \frac{1}{n}YFM \)、\( Y = U \)となる。

    よって確率の変数変換の際のヤコビアンは、

    $$ \begin{bmatrix} \frac{\partial X}{\partial F} & \frac{\partial X}{\partial U} \\ \frac{\partial Y}{\partial F} & \frac{\partial Y}{\partial U} \end{bmatrix} = \begin{bmatrix} \frac{1}{n}Um & 0 \\ 0 & 1 \end{bmatrix} = \frac{m}{n}U $$
    となる。

    FとUの結合分布の密度関数は
    $$ h(x, y) = \frac{ x^{ \frac{m}{2} -1} e^{ – \frac{x}{2} }} { 2^{ \frac{m}{2}} \Gamma \left ( \frac{m}{2} \right )} \frac{ y^{ \frac{n}{2} -1} e^{ – \frac{y}{2} }} { 2^{ \frac{n}{2}} \Gamma \left ( \frac{n}{2} \right )} $$
    となる。
    確率の変数変換を行い整理すると、
    $$ g(f,u) = h(x, y) |J | \\
    = h \left ( \frac{1}{n}ufm , u \right ) \left | \frac{m}{n}u \right | \\
    = \frac{ u^{\frac{m+n}{2} -1} n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} e^{-\frac{u}{2} \left ( \frac{mf}{n} + 1 \right ) } }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) }
    $$

    ここでfを固定して、\( t = \frac{1}{2} \left ( \frac{mf}{n} + 1 \right ) u \)として変数変換を行う。なお、\( u = \frac{2t}{ \frac{mf}{n}+1 } \)、\( \frac{du}{dt} = \frac{2}{ \frac{mf}{n} + 1 } \)となる。

    $$ g(f) = \int_0^{\infty} g(f, u)du \\\\
    = \frac{ n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \int_0^{\infty} e^{-\frac{u}{2} \left ( \frac{mf}{n} + 1 \right ) } u^{\frac{m+n}{2} -1} du
    $$

    ここで積分の中についてのみ注目すると、

    $$ \int_0^{\infty} e^{-\frac{u}{2} \left ( \frac{mf}{n} + 1 \right ) } u^{\frac{m+n}{2} -1} du \\\\
    = \int_0^{\infty} e^{-t} \frac{2}{ \frac{mf}{n} + 1 } \left ( \frac{2t}{ \frac{mf}{n} + 1} \right )^{ \frac{m+n}{2}-1 }dt \\\\
    = \left ( \frac{2n}{mf+n} \right )^{\frac{m+n}{2}} \int_0^{\infty} t^{ \frac{m+n}{2}-1}e^{-t}dt
    $$

    よって、

    $$ g(f) = \frac{ n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \left ( \frac{2n}{mf+n} \right )^{\frac{m+n}{2}} \int_0^{\infty} t^{ \frac{m+n}{2}-1}e^{-t}dt $$

    ガンマ関数の定義より、
    $$ \Gamma ( \alpha ) = \int_{0}^{\infty} u^{ \alpha – 1} e^{-u}du \\
    \alpha – 1 = \frac{m+n}{2} – 1 \\
    \alpha = \frac{m+n}{2}
    $$
    であるから、これを考慮し整理する。

    $$ g(f) = \frac{ n^{-\frac{m}{2}} m^{\frac{m}{2}} f^{\frac{m}{2} -1} }{ 2^{ \frac{m}{2} } 2^{ \frac{n}{2} } \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \left ( \frac{2n}{mf+n} \right )^{\frac{m+n}{2}} \Gamma \left ( \frac{m+n}{2} \right ) \\\\
    = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ f^{ \frac{m}{2} -1}}{ (mf+n)^{ \frac{m+n}{2} } }
    $$

    これはF分布の確率密度関数である。

    F分布の平均

    $$ E(X) = \int_0^{\infty} x \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx \\\\
    = \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2}}}{ (mx+n)^{ \frac{m+n}{2} } } dx
    $$

    ここで、\( y = \frac{x}{mx+n} \)と変数変換すると、\( x = \frac{ny}{1-my} \)、\( \frac{dx}{dy} = \frac{n}{(1-my)^2} \)、\( mx + n = \frac{nym}{1-my} + n = \frac{n}{1-my} \)であるから、これらを代入すると、

    $$ E(X) = \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ \left ( \frac{ny}{1-my} \right )^{\frac{m}{2}} }{ \left ( \frac{n}{1-my} \right )^{\frac{m+n}{2}} } \frac{n^2}{n (1-my)^2 }dy \\\\
    = \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } y^{\frac{m+n}{2}} \left ( \frac{ny}{1-my} \right )^{-\frac{n}{2}} \frac{1}{n} \left ( \frac{n}{1-my} \right )^2 dy
    $$

    ここで、\( z=my \)と変数変換すると、\( y = \frac{z}{m} \)、\( \frac{dy}{dz} = \frac{1}{m} \)であるから、これらを代入すると、

    $$ E(X) = \int_0^{\infty} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \left ( \frac{z}{m} \right )^{\frac{m+n}{2}} \left ( \frac{ny}{1-z} \right )^{-\frac{n}{2}} \frac{1}{n} \left ( \frac{n}{1-z} \right )^{2}\frac{1}{m}dz \\\\
    = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) }m^{-1}n \int_0^{\infty} z^{ \frac{m}{2}}(1-z)^{\frac{n}{2}-2}dz
    $$

    ベータ関数、\( B(\alpha , \beta) = \int_0^{1} t^{\alpha – 1}(1-t)^{\beta -1}dt \)の定義より、\( \alpha = \frac{m}{2} + 1 \)、\( \beta = \frac{n}{2} -1 \)であるから、
    $$ B \left ( \frac{m}{2}+1, \frac{n}{2} -1 \right ) = \int_0^{\infty} z^{ \frac{m}{2}}(1-z)^{\frac{n}{2}-2}dz $$
    ベータ関数とガンマ関数の関係式より、
    $$ B \left ( \frac{m}{2}+1, \frac{n}{2} -1 \right ) = \frac{ \Gamma \left ( \frac{m}{2} + 1 \right ) \Gamma \left ( \frac{n}{2} – 1 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } $$
    となる。
    よって、
    $$ E(X) \\
    = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) }m^{-1}n \frac{ \Gamma \left ( \frac{m}{2} + 1 \right ) \Gamma \left ( \frac{n}{2} – 1 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } \\
    = \frac{n}{n-2}
    $$
    となる。

    F分布の分散

    分散は
    $$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 $$
    より求まる。

    $$ E(X^2) = \int_0^{\infty} x^2 \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx \\\\
    = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \int_0^{\infty} \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx
    $$

    ここで、先ほどと同様に変数変換を行うが、ここでは積分以降にのみ注目する。

    $$ \int_0^{\infty} \frac{ x^{ \frac{m}{2} -1}}{ (mx+n)^{ \frac{m+n}{2} } } dx \\\\
    = \int_0^{\infty} \left ( \frac{ny}{1-my} \right )^{ \frac{m}{2}+1 } \left ( \frac{n}{1-my} \right )^{ -\frac{m+n}{2}}\frac{1}{n} \left ( \frac{n}{1-my} \right )^{2} dy \\\\
    = \int_0^{\infty} \left ( \frac{ny}{1-my} \right )^{ -\frac{n}{2}+1 } \frac{1}{n} \left ( \frac{n}{1-my} \right )^{2} dy
    $$

    ここで、\( z = my \)と変数変換を行う。なお、\( y = \frac{z}{m} \)、\( \frac{dy}{dz} = \frac{1}{m} \)であるから、これらを代入すると、

    $$ = \int_0^{\infty} \left ( \frac{z}{m} \right )^{ \frac{m+n}{2}} \left ( \frac{nz/m}{1-z} \right )^{ – \frac{n}{2} +1} \frac{1}{n} \left ( \frac{n}{1-z} \right )^{2} \frac{1}{m} dz \\\\
    = m^{ -\frac{m}{2} -2 }n^{- \frac{n}{2} + 2} \int_0^{\infty} z^{ \frac{m}{2} +1 }(1-z)^{ \frac{n}{2} -3 }dz
    $$

    ベータ関数、\( B(\alpha , \beta) = \int_0^{1} t^{\alpha – 1}(1-t)^{\beta -1}dt \)の定義より、\( \alpha = \frac{m}{2} + 2 \)、\( \beta = \frac{n}{2} -2 \)であるから、

    $$ B \left ( \frac{m}{2}+2, \frac{n}{2} -2 \right ) = \int_0^{\infty} z^{ \frac{m}{2} +1 }(1-z)^{ \frac{n}{2} -3 }dz $$

    よって、

    $$ E(X^2) = \frac{ \Gamma \left ( \frac{m+n}{2} \right ) m^{ \frac{m}{2} } n^{ \frac{n}{2} } }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } m^{ -\frac{m}{2} -2 }n^{- \frac{n}{2} + 2} B \left ( \frac{m}{2}+2, \frac{n}{2} -2 \right ) $$

    ベータ関数とガンマ関数の関係式より、
    $$ B \left ( \frac{m}{2}+2, \frac{n}{2} -2 \right ) = \frac{ \Gamma \left ( \frac{m}{2} + 2 \right ) \Gamma \left ( \frac{n}{2} – 2 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } $$
    となる。
    よって、

    $$ E(X^2) = \frac{n^2}{m^2} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ \Gamma \left ( \frac{m}{2} + 2 \right ) \Gamma \left ( \frac{n}{2} – 2 \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } \\\\
    = \frac{n^2}{m^2} \frac{ \Gamma \left ( \frac{m+n}{2} \right ) }{ \Gamma \left ( \frac{m}{2} \right ) \Gamma \left ( \frac{n}{2} \right ) } \frac{ \left ( \frac{m}{2} + 2 \right ) \frac{m}{2} \Gamma \left ( \frac{m}{2} \right ) }{ \Gamma \left ( \frac{m+n}{2} \right ) } \frac{ \Gamma \left ( \frac{n}{2} \right ) }{ \left ( \frac{n}{2} – 2 \right ) \left ( \frac{n}{2} -1 \right ) } \\\\
    = \frac{n^2 \left ( \frac{m+2}{2} \right ) \frac{m}{2} }{m^2 \left ( \frac{n-4}{2} \right ) \left ( \frac{n-2}{2} \right ) } \\\\
    = \frac{n^2(m+2)}{m(n-4)(n-2)}
    $$

    分散の式にあてはめると、
    $$ V(X) = E \left ( X^2 \right ) – \left [ E(X) \right ]^2 \\
    = \frac{n^2(m+2)}{m(n-4)(n-2)} – \frac{n^2}{(n-2)^2} \\
    = \frac{2n^2(n+m-2)}{m(n-4)(n-2)^2}
    $$
    となる。

    参考文献

    [1] Distributions related to the normal distribution
    [2] 鈴木・山田 (1996), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃

    [R]ボージョレ・ヌーボーのコメントに対してLDATSパッケージを使って時系列トピックモデルを扱う

    はじめに

    先日、某勉強会でLTをしました。その際に10秒だけ紹介したRのパッケージについて記事を書いてみようと思います。

    LDATSパッケージについて

    時系列でのトピックモデルを推定することができるパッケージです。
    やっていることとしてはLDAでトピックを推定して次元を減らし、そのトピックの多変量時系列に関してベイズ手法による変化点検知のためのパラメータ推定を行っているようです。GitHubの該当しそうなソースコードに多変量のデータに対するsoftmax関数での回帰をやっているとの記述がある。(multinomial Bayesian Time Series analysis

    元となっている論文を見る限り、BoW(Bag of Words)を想定して作っておらず、20~30程度のグループからなるデータに対して適用するのがちょうど良いです。アクセスログのページカテゴリや、マーケティングの顧客セグメントであればそんなに数は多くないので扱いやすいと思います。

    データ

    Webサイトから集めてきたボージョレ・ヌーボーのキャッチコピー14年分を今回は扱います。実は販売店側のキャッチコピーとワイン委員会が決めた評価が存在します。私の知っている世界は販売店側のキャッチコピーだけでした。

    試してみた

    今回はとにかく動くことだけを考えて、汚いコードとなっております。やっていることとしては、キャッチコピーを販売側とワイン委員会側のものを一つにつないで、数字を正規表現で「数字」に変換し、RMeCabで形態素解析をし、LDATS向けの形式のデータを作成していきます。
    途中で、日本語の文字化け問題を回避するためにGoogle翻訳を使って単語名を置き換えています。
    1時系列につき1文書となるようにデータを作っていく必要があるのですが、今回はボージョレ・ヌーボーのキャッチコピーなので最初から1時系列につき1文書となっているため都合が良いです。
    データとソースコードはこちら

    こちらは論文の図と同じものだとドキュメントの説明にあったので、論文の説明を見る限り、表すものとしては以下のようです。

    • 一番上の積み上げグラフはトピックごとの単語の割合を表しています。
    • 二番目の折れ線グラフはLDAによって推定されたトピックの時系列推移です。
    • 三番目のヒストグラムは二番目の時系列における変化点を集計したものです。
    • 四番目の折れ線グラフはモデルが推定したトピック割合の変化点の前後での推移です。

    今回の図では文字が潰れていて見にくいですが、

    • トピック1はボキャブラリーが比較的リッチなコメント(「フルーティー」「フレグランス」「複雑」)
    • トピック2は数字を用いたコメント(「何年に一度の!」みたいな)
    • トピック3はボキャブラリーが貧相なコメント(「すごい!」みたいな)

    のようです。
    二番目の折れ線グラフを見る限り、周期的に数字を用いたコメントが現れているように思われます。四番目の折れ線グラフの変化点を見る限り、近年は数字を用いたコメントが相対的に減ってきて、リッチなボキャブラリーになってきているようです。

    おわりに

    時系列トピックモデルをカジュアルに試せる面白そうなパッケージだなと思い、LDATSパッケージを触ってみましたが、そもそもBoWなどを想定して作られているパッケージではないので、単語数が多いような分析ではそもそも可視化ができず使いにくいだろうなと思いました。マーケティングなどでユーザーのセグメントの推移を分析したい場合などにちょうど良いのだろうと思われます。

    参考情報

    [1] Long‐term community change through multiple rapid transitions in a desert rodent community
    [2] Latent Dirichlet Allocation coupled with Bayesian Time Series analyses
    [3] Package ‘LDATS’

    RStudioをdockerで使える、Rockerを触ってみた

    はじめに

    dockerは会社のエンジニア向けの勉強会で紹介されていて、存在は知ってはいたが使っていませんでした。先日参加したTokyo.RでRockerというものが紹介されており、よし使ってみようと思うに至りました。
    今回は遅ればせながら、Rockerを使ってみて、何かを回して、保存するという一連のチュートリアルをやってみます。

    用語のざっくり理解

    専門家の方に怒られたらあれですが、私の理解はこんな感じです。

    • Dockerfile
      • OS、言語、ライブラリ、バージョンなどをレシピの食材っぽい感じでテキスト形式で書かれたもの
    • イメージ
      • コンテナの元。Dockerfileをもとにdockerで生成できる。あるいはdocker hubとか言うのサイトから入手できる。ここからコンテナ(RStudioとかを回したい環境のこと)を起動できる。
    • コンテナ
      • RStudioとかを起動させておく実行環境。
    • 永続化
      • ファイルの内容が消えないように残すこと。立ち上げたコンテナを永続化しないと都度データやインストールしたパッケージが消える。

    Rockerはdocker hubで提供されているイメージです。RockerのGitHubにはDockerfileも当然あげられていました。自分であれやこれや設定を追加するにはこれを参考にいじることになるのだと思います。

    とにかくやってみる

    まずは、dockerをインストールします。

    Install Docker Desktop on Mac
    docker自体はだいぶ前にインストールしたので、入れ方を忘れていますが、ここを見ればわかるはず。

    ターミナルで以下を実行

    ブラウザなどで以下を表示

    DALEXのサンプルコードを回してみる。DALEXはTokyo.Rでも紹介されている機械学習の解釈可能性に関する手法をあらかた実践できるパッケージです。
    moDel Agnostic Language for Exploration and eXplanation

    以下の図はDALEXなどを使ってGLMで推定した結果の各特徴量の影響度をプロットしたものです。

    ターミナルでコンテナのIDを確認する。

    分析結果やソースコードやデータを保存したいので、コンテナの永続化を行います。

    ここではユーザー名をperio、イメージ名をrocker_hadleyにしています。

    永続化したコンテナのイメージがきちんと存在するか確認してみます。

    ありますね。

    さて、先程のコンテナを止めてみましょう。コンテナIDを指定するだけです。

    永続化したコンテナを読み込んでRstudioを起動してみます。

    実は、誤って最新の状態で保存せずに再度立ち上げてしまいましたwDALEXのサンプルコードを回すために色々パッケージインストールしたので、それが消えてしまった。残念。
    元のimageにはないRDataがあるので、一応、永続化はできているので安心です。うむ、次からはきちんと保存をしよう。

    これでどのPCでも同じ条件で分析していけるので、良さそうです。
    ただ、一部で入らないパッケージ(sfパッケージとか)があったので、より良いDockerfileを探し求めたいとも思いますね。

    u_riboさんがsfパッケージが動くdocker imageを教えてくれました。これまで、一つの環境にギュッとパッケージを詰め込んでやってきたんですが、目的に応じてdocker imageを選ぶのが良いのでしょう。

    せっかくなので先日作成した事故物件分析用のコードを回してみましょう。

    mapviewもインストールできて、問題なく使えました。u_riboさんありがとう。

    参考情報

    [1] The Rocker Project Docker Containers for the R Environment
    [2] The Rocker Images: choosing a container
    [3] rocker/tidyverse
    [4] How to save data
    [5] DALEXverse and fraud detection

    2019年に読んだデータ分析系の本の振り返り(21+1冊)

    はじめに

    2020年、あけましておめでとうございます。年末に自分自身を振り返ろうと思ったのですが、結局データ分析と勉強しかしていないわけで、書籍を振り返ろうと思うに至りました。私の知り合いのデータサイエンティストはだいたい全冊持っているであろうと思われますが、良い本だと思うので思い出していただければ幸いです。

    1.『ベイズモデリングの世界』(岩波書店)

    基本的に階層ベイズモデルを使って、個体ごとの異質性を考慮した分析手法が提案されています。前半はオムニバス形式で様々な先生がモデルの適用について執筆されており、後半では伊庭先生による階層ベイズモデルの講義になっています。途中でスタイン統計量による縮小推定の話があげられ、柔軟なモデリングのためには「階層化した方が少なくとも望ましい推定量が得られる」という数学的証明を捨てることもやむを得ないと書かれています。

    2.『トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ) 』(コロナ社)

    この本はトピックモデルの教科書というよりも、ベイズ推定の教科書という側面が強い印象があります。途中で出てくる数式は流し読みするのは難しく、最低2冊以上のノートが別途必要になると思います。一度でもLDAのパラメータを導出してみたいという方には良い教科書だと思います。疑似コードが提供されているので、それをもとにRやPythonでコーディングしていけば、一番シンプルなLDAが非常に短い行で実行できてしまうことに驚かれるかもしれません。人間が手を動かして推定アルゴリズムを導出しているからこそ、短いコードで済むということを実感できるはずです。

    3.『構造的因果モデルの基礎』(共立出版)

    グラフィカルなアプローチで因果推論を扱っている書籍です。Judea Pearl流の因果推論アプローチについて記すことを目的に書かれています。基礎と書かれていますが決して簡単ではありません。ただ、扱われる数学のレベルとしては確率と線形代数がわかれば大丈夫だと思われます。余談ではありますが、1章の相関関係と因果関係の事例紹介で「おむつとビールの話」が都市伝説ではなくきちんと記事としてWall Street Journalという雑誌に掲載されていたことが明らかにされています。

    4.『現場で使える!PyTorch開発入門 深層学習モデルの作成とアプリケーションへの実装 (AI & TECHNOLOGY)』(翔泳社)

    PyTorchを触ったことがないが、深層学習の手法について知っている層を対象とした本です。6章まではGoogleのColabで動かせるのでGoogleに課金することなく深層学習による回帰、CNN、GAN、RNN、Encoder-Decoderモデル、ニューラル行列因子分解をPyTorchで試すことができます。写経したものはこちら。転移学習や高解像度化や画像生成、文章のクラス分類、文書生成、機械翻訳などもできるので、PyTorchでこれくらいの量をコーディングしたらこれくらいのことができるのかという学びや、他の人の書いたPyTorchコードを読みやすくなるなどの便益は十分にあると思いました。

    5.『作ってわかる! アンサンブル学習アルゴリズム入門』(シーアンドアール研究所)

    会社で行っているPythonもくもく会用に買った本で、scikit-learnを使わずに機械学習のアルゴリズム(アンサンブル系)をコーディングするための本です。pythonのコードについて逐次、細かい解説が行われているわけではないので、1行1行自分でコメントを加えながら写経をしていけば力が付くという本かなと思われます。sklearnはそれはそれで素晴らしいですが、こういう本でフルスクラッチで修行できるのはいいですね。

    6.『数理統計学―基礎から学ぶデータ解析』(内田老鶴圃)

    統計検定1級を合格された方のブログで紹介されていた教科書です。理系の大学生レベルの数学知識があれば、数理統計学の基礎を学べると思います。中心極限定理の証明や、様々な分布の期待値や分散、様々な分布の性質について数式を用いてしっかり理解することができます。数式もほどよく端折られているので、無論ですがノートが数冊必要になります。各章毎にある練習問題も解くことで力が付くと思います。日本の大学の授業の教科書がこれだったらジェノサイド(再履修者の大量発生)が起きるんだろうなと思ってしまった。

    7.『44の例題で学ぶ統計的検定と推定の解き方』(オーム社)

    統計の検定に関してだけ扱った珍しい本です。第3部までは統計学の普通の教科書ですが、それ以降であらゆる検定の例題が44件も載せられています。パラメトリックな検定から、ノンパラメトリックな検定まで幅広く扱われています。一番気にいっているのは仮説検定法の分類の表です。これさえあれば、どのデータに対してどの検定を行えばいいかが一目瞭然です。

    8.『わけがわかる機械学習 ── 現実の問題を解くために、しくみを理解する』(技術評論社)

    機械学習の原理を手早く数式を交えて学べる本です。かゆいところに手が届いていると言うか、既出の教科書では捨象されがちな、条件付き確率における2変数以上の条件づけでの表現に紙面を割いていたりしてくれるのが嬉しいです。ある程度数学の話はわかるが、だいぶ忘れているビジネスパーソンには大変にありがたいコンテンツと言えると思います。ベイズ線形回帰に関しても行列を用いた、わかりやすい導出方法が紹介されています。またコラムで紹介されている、測度論にどう向き合えばいいかの著者の見解は参考になります。

    9.『Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

    R言語とstanを用いてベイズ統計学を入門レベルから学べる本です。各トピックごとにそれなりの紙面が割かれています。例題も豊富にあるので、線形回帰・MCMC・情報量基準・階層ベイズモデルまで、ベイズ統計学を基礎から応用までしっかりと学べると思います。youtubeで著者の講義も配信されているので、留学気分を味わえます。

    10.『scikit-learnとTensorFlowによる実践機械学習』(オライリージャパン)

    2019年に日本で開かれたML SummitでTFの開発者がおすすめしていた教科書です。前半部分で機械学習の入門から応用までをわかりやすい説明で学ぶことができます。数式は少ないですが、図とソースコード(Python)がちりばめられており、手を動かして理解を進めることができます。後半部分はTensorFlowを用いた深層学習の基礎を同様に手を動かして学ぶことができます。ただ、TFのバージョンも変わってきているので前半の説明をアテにして読むのも良いと思います。

    11.『AIアルゴリズムマーケティング 自動化のための機械学習/経済モデル、ベストプラクティス、アーキテクチャ (impress top gear)

    マーケティングへのデータサイエンスの適用に関する珍しい書籍です。ソースコードはついていないですが、業務で使う際のアイデアが手に入ることもあります。一般的な回帰、生存時間分析、オークション、アトリビューション分析、アップリフトモデリング以外にも、情報検索やレコメンデーションやトピックモデルなどマーケティングながら学際的なトピックも扱われています。レコメンドなどで使われる、ランク学習に関して詳しく書かれた書籍をあまり知らないので、この本はその点においてもありがたい本でもあります。

    12.『入門 統計的因果推論』(朝倉書店)

    ほぼ全ての章でグラフィカルなアプローチで因果推論を扱っています。例題も豊富なので、一つ一つ丁寧にやれば理解が捗ります。おそらく、例題の多さを含め一番丁寧にd分離性、do演算子、バックドア基準、フロントドア基準に関する説明をしてくれている本なのかなと思いました。グラフでの因果推論に関して初めての人でも、確率さえ知っていれば読み進めることができるはずです。また、途中で操作変数法の紹介もされ、経済学出身者としては読みやすい。ただ、傾向スコアのくだりや、DIDなどのくだりはあまり出てきません。あと、やってないですが章末の練習問題に対するSolution Manualが提供されているようです。

    13.『実践 ベイズモデリング -解析技法と認知モデル-』(朝倉書店)

    ベイズモデリングを様々な事例に適用する方法がオムニバス形式で記された本です。ワイブル分布や異質性を考慮した二項分布、無制限複数選択形式のアンケートデータに対する手法、トピックモデル、項目反応理論などが扱われています。マーケティングの実務で使える事例が多いように感じました。こちらはサポートサイトでRコードとstanコードが提供されています。あと、appendixにあるプレート表現の見方も参考になります。

    14.『機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

    機械学習などで用いるベイズ推論を扱った教科書です。入門とありますが、入門者は書かれた数式をそのまま見ていても頭に入らないのではないでしょうか。手を動かしてなんぼの本だと思います。ノート2冊は絶対に必要です。たぶん、数式の展開を丁寧に記すと倍以上の厚みの本になると思います。各々のモデルに関してグラフィカルモデルが記されているのや、サンプルコードとしてGitHubにJuliaで書かれたソースコードが提供されているのも良いです。

    15.『その問題、数理モデルが解決します』(ベレ出版)

    物語形式で、様々な問題に対して数理モデリングのアプローチが紹介されています。途中でマッチング理論やゲーム理論やオークションなども登場することから、経済学出身者も喜ぶ内容かもしれません。社会人になってからナッシュ均衡という言葉が書かれた本は中々出会って来なかった。

    16.『ヤバい予測学 ― 「何を買うか」から「いつ死ぬか」まであなたの行動はすべて読まれている』(CCCメディアハウス)

    2013年と結構古い本ですが、データ分析を様々な事象に対して適用した事例紹介本です。アップリフトモデリングへの言及もあり、こういったものに関して日本は何年も遅れてブームが来るんだなという実感を与えてくれた本でもありました。appendixに分析事例が147個ほどあげられているのも参考になります。

    17.『たのしいベイズモデリング2: 事例で拓く研究のフロンティア』(北大路書房)

    主にstanを用いたベイズモデリングによる分析事例が1と2で38本もオムニバス形式で載っています。ほとんどの事例で階層ベイズモデルが扱われています。2では若干マーケティングに近い内容の題材も扱われ、データサイエンティストの人にも嬉しい内容かもしれません。もちろんデータとstanとRのコードがサポートサイトで提供されています。

    18.『カルマンフィルタ ―Rを使った時系列予測と状態空間モデル― (統計学One Point 2)』(共立出版)

    状態空間モデルで時系列予測を行うための手法が記されている本です。RのKFASパッケージが全面に渡って扱われています。トレンドを考慮したり、カレンダー効果を追加したり、共変量を追加したりなど様々なアプローチが紹介されコードも伴っているわけですから、業務でも抜群に役に立ちました。

    19.『機械学習のエッセンス -実装しながら学ぶPython,数学,アルゴリズム- (Machine Learning)』(SBクリエイティブ)

    自分のいる会社で最低限の数学がわかると思われる若いメンバーに買ってもらうように言っている本です。微積分・線形代数だけでなく、カルシュ・キューン・タッカー条件(最適化数学)に関しても扱ってくれているので、ここで出てくる数学がわかれば大体の論文に立ち向かえると思います。さらに、Pythonの基礎もこれで学ぶことができるので一石二鳥な素敵な本ですね。また、最後の方でスクラッチでアルゴリズムを書くパートがあり、こちらも勉強になります。

    20.『機械学習のための特徴量エンジニアリング ―その原理とPythonによる実践 (オライリー・ジャパン)』(オライリー・ジャパン)

    機械学習における前処理の指針を与えてくれる本です。Pythonのコードが提供されています。例えばですが、「テストデータにだけある、新しい単語は取り除いてしまえばいい」などの細かいアドバイスが何気に嬉しいです。「Effectコーディング」「特徴量ハッシング」「ビンカウンティング」「バックオフ」「leakage-proof統計量」などは読むまで知らないところだったので勉強になりました。

    21.『データサイエンスのための統計学入門 ―予測、分類、統計モデリング、統計的機械学習とRプログラミング』(オライリージャパン)

    データ分析の仕事をする上で最低限必要な知識を幅広く抑えることができる本です。数式は少なく、ところどころ出てくるコードはR言語です。参考文献などがブログだったりするため厳密さがめちゃあるわけではないですが、業務で使う分には問題ないと思います。分類問題において、AUCなどの評価指標だけでなく、予測値自体の探索的分析のすすめなどが書かれており、参考になりました。また、特徴量エンジンとしてのk-NN法の話も面白いと思いました。

    [+α]『プログラマのためのGoogle Cloud Platform入門 サービスの全体像からクラウドネイティブアプリケーション構築まで』(翔泳社)

    Google Cloud Platformを初めて触るデータ分析者にはちょうど良い本です。説明もわかりやすいので、いきなりアカウントを作ってドキュメントを解読するよりかは戸惑いは減るはずです。この本を土台に、GCS・GCEを駆使してML系のAPIを呼び出して使うなどの最低限の操作は私でもできるようになりました。GCPの画面や機能もどんどん変わっていくので書籍を買ってもアレなんですが、歴史的な背景も若干記述されているので、それはそれで勉強になります。ただ、エンジニアにこの本を買うべきか聞いた際にネガティブな意見があったのですが、たぶん現役プログラマからすると簡単過ぎるからなんだろうなと思います。

    終わりに

    2019年もぼちぼち勉強できましたが、2020年もこれまで同様にノートとペンを大事にする勉強を続けていき、コーディングも分析ももっともっと数をこなして会社や社会に求められるようなデータ分析官を目指していこうと思います。あぁ、英会話などの勉強をする時間を作るのが難しい。

    R advent calendar 2019 RSelenium、jpmesh、sfパッケージで東京23区の事故物件を分析してみよう!

    はじめに

    今回で3回目となるR advent Calendarですが、前回は「一発屋芸人の検索トレンドのデータ」を扱い、前々回は「ポケモンのデータ」を扱いました。今回は人の命に関わるようなデータを扱ってみたいと思い、某サイトから東京都の23区内における事故物件の住所と詳細を集めてきました。どのようなエリアが事故が起きやすいのかの分析を行います。(以下では、事故物件をAP(Accident Property)と呼びます。)

    分析工程

    ・データの収集
    ・データの整形
    ・可視化
    ・分析

    データの収集

    APに関する情報を某サイトより集める必要があります。そこで必要なライブラリとしては、RSeleniumやtwitteRがあげられます。
    twitteRが必要な理由は、APに関するサイトにAPの一覧ページがなく、公式アカウントがAPに関するページのリンクを投稿しているところにあります。ただ、私が以前使っていた頃とはTwitterAPIの仕様が変わり、3ヶ月よりも前の情報にアクセスできなくなっていました。そのため、今後のデータに関してはTwitterAPIでいいのですが、過去のものに関しては別アプローチが必要となります。
    また、APに関するサイトはJavaScriptで地図が表示されているだけなので、RSeleniumを使って地図をクリックさせ、表示された情報をスクレイピングするなどの処理が必要となります。
    当初の想定ではTwitterのデータ収集、リンク先の情報をRSeleniumでスクレイピングするだけの簡単な仕事だと思っていたのですが、過去のデータにアクセスできないので、地図上で一つ一つ見つけていくためにRSeleniumだけで頑張ることにしました。(私の過去のアドカレ史上、一番面倒なデータとなりました。)

    誰もやらないと思うのですが、一応手順を記しておきます。

    RSeleniumだけでMapからAPの情報を抽出するための手順
    1.都内の住所一覧を収集
    2.検索窓に住所を入力
    3.検索結果一覧の上位5件をクリック
    4.一度地図を引くことでAPを広い範囲で捉えれるようにする
    5.APのマークの要素を取得し、1件ずつクリックし、表示されたAPの情報をデータフレームに格納する

    こちらが取得できたデータです。

    RSeleniumでAPの
    ・住所
    ・発生時期(フリーテキスト)
    ・AP詳細
    ・投稿日
    を集めることができるので、その住所データに対して、Yahoo!のジオコードに関するAPIを利用します。(利用申請が必要なはずです。4年前くらいに申請していたのでそのまま使えました。)
    Yahoo!のAPIを使えば、住所から緯度経度の情報を取得することができます。

    APの緯度経度がわかれば、jpmeshパッケージを用いて1kmメッシュやら10kmメッシュやらのメッシュデータに変換することができます。
    jpmeshを用いてメッシュデータに変換し、メッシュ単位でAPの発生件数を集計します。

    データ収集用のソースコードは思いのほか長くなってしまったので、GitHubにあげておきました。
    https://github.com/KamonohashiPerry/r_advent_calendar_2019

    ここで再度、手順を整理しておきましょう。

    Twitterに出てきたものだけを取得(直近3ヶ月〜)する場合、
    run_tweet_collect.RでTweetを収集

    run_selenium.RでAPの情報をスクレイピング

    run_map_api.Rで住所から緯度経度の取得

    making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

    直近3ヶ月以前のものを取得する場合、
    run_selenium_from_map.Rで地図上から直接APの情報を取得する

    making_mesh_data_and_download_other_data.Rで緯度経度からメッシュデータへの変換、その他の人口データや地価データと接続をします。

    データの整形

    APの1kmメッシュデータを手に入れたら、kokudosuuchiパッケージを使って国土地理院の収集したデータをつなぎこみます。手順としては、以下のとおりです。

    まずは推計人口というそのエリアの人口の予測値です。今回は2010年のものを抽出しました。こちらは1kmメッシュのデータなので、変換することなく使えて都合が良いです。

    続いて、2015年の公示地価を抽出しました。

    こちらはメッシュデータではないので、緯度経度の情報から1kmメッシュのデータに変換する必要があります。後で行います。

    単純集計・可視化

    今回のデータセットのデータ数は3919件です。本当は7000件以上はあると思われますが、マップから取ってくるという勝手上、なかなか全てを取り切ることができませんでした。

    まずは、1kmメッシュごとのAP発生件数のヒストグラムです。

    1kmメッシュにおける人口のヒストグラムです。

    公示地価のヒストグラムです。

    1kmメッシュにおける人口あたりのAP件数のヒストグラムです。

    分析

    ここでは、色々な軸でAPのデータに向き合ってみようと思います。

    APの発生件数の集計

    人口が多いところがAPの発生件数が多いところだと思われますが、とりあえず確認します。

    世田谷区は最も人口が多いことから、AP発生件数では一番となっています。続いて、歌舞伎町などがある新宿が来ています。しかしながら、人口に占めるAP発生件数で言うと、港区がかなり高く出ているのがわかります。

    人口あたりのAP件数

    ここでは、メッシュデータをsfパッケージ用のオブジェクトに変換して、1kmにおける人口あたりのAP発生割合を可視化しています。

    こちらのmapviewパッケージで作ったマップはインタラクティブにいじることができます。ぜひ関心のあるエリアでいじってみてください。

    1kmメッシュ人口あたりのAP発生件数(×100)の可視化

    1kmメッシュでのAP発生件数の可視化

    比率ベースで、色の明るいメッシュのところを見ると、港区、中央区、新宿区、渋谷区などがAPが発生しやすいようです。件数ベースで言うと新宿が一番多いですね。
    一番色が明るい港区はてっきり六本木ではないかと思ったのですが、新橋から日比谷にかけたエリアでした。会社員による自○が多いようです。恐ろしいものです。

    APの名前の集計

    APの名前を集計してみます。これは別にこの名前だからAPになりやすいというわけではなく、単純に数が多いだけの可能性がありますし、実際にそうだろうと思われます。AP発生率を知るには、APではないものも含めた全物件名に占めるAP発生件名を手に入れないといけませんが、全物件名を収集するのが難しいことから単純に頻度の集計となります。今回は、wordcloud2パッケージを使って、ワードクラウドにしてみます。文字が大きいと頻度が高いものとなります。

    ハイツ、荘、コーポ、マンション、号棟、アパート、ハウスなどが多く出現しているようです。ただ、物件の名前としても頻度が高いとも考えられますね。

    地価と人口あたりのAP発生件数の関係

    ここでは地価のデータとAP発生の関係性について見てみます。

    地価の階級値(10個のパーセンタイルに分割)を横軸に、縦軸に人口あたりAP発生数をおくと、地価が上がるに従い人口あたりAP発生数が高まる傾向があります。これは、人口密度が高く地価の高いところではAPが発生しやすいということを示しているのではないでしょうか。人口密度が高いと地価があがる、人口密度が高いと治安が悪くなるという可能性が考えられます。

    APの詳細の集計

    ここではAPになってしまった詳細の内容について先ほどと同様に形態素解析を行いワードクラウドにしてみます。

    どうやら孤独死が多いようです。高齢者の人口構成比が関係しているのだろうと思われます。

    APの発生時期に関するテキストマイニング

    ここでは、発生時期に含まれる四桁の数字を集計して、何年くらいのAPが多いのかをざっくりと掴みます。

    どうやら昔のデータはあまり登録されていないようです。記憶が確かではないかもしれませんし、古すぎるものは消されているのかもしれませんね。あのサイトはユーザー生成コンテンツ(UGC)なので、投稿する人はそこまで昔のことをわざわざ投稿しないのかもしれないですね。

    APの詳細に関するテキストマイニング

    ここではトピック数10として、topicmodelsパッケージを使いLDAを行います。

    なかなか恐ろしいキーワードが多いですが、なんとなくですがうまく分類されているのではないかと思われます。

    トピック1は男性の不幸
    トピック2は不動産屋に言われた告知事項
    トピック3は孤独死
    トピック4は病死
    トピック5は火災・転落・事故
    トピック6は事故のあった建物に関する記載
    トピック7は腐乱した事例
    トピック8は建物に関して不明であることの記載
    トピック9は心理的瑕疵あり
    トピック10は自○

    となっているように思われます。まさかこのようなデータにトピックモデルを使うことになるとは。

    おわりに

    今回はR言語のみを用いて、APに関するデータを収集し、地図にプロットしたり他のメッシュデータとつなぎ合わせて分析をするなどしました。APが発生しやすいエリア、APと地価との関係、APのテキストマイニングなど興味深い結果が得られたと思います。
    一つ残念なのは、時系列情報がフリーテキストなので、APがどのエリアでどの頻度で発生していくのかの分析のコストが高く、今回は時系列情報を用いた分析にチャレンジできませんでした。
    今後はタクシーの需要推定の分析で行われているように、メッシュ単位でのAP発生確率の推定などを機械学習で行えると面白いなと思います。どなたか一緒にアノテーションしましょう!

    それでは、どうか良い年末をお過ごし下さい!
    メリークリスマス!

    参考情報

    数多くの方々の記事を見てどうにか仕上げることができました。感謝します。

    [1]【追記あり】sfパッケージでシェープファイルを読み込んでmapviewパッケージで可視化するまで
    [2]How to use mesh cord in R
    [3]Rを使ってワードクラウドを作ってみました
    [4]国土数値情報ダウンロードサービスWeb APIからデータを取得するためのRパッケージです
    [5]東京の地価公示データを眺める
    [6]Chapter 1 Introduction to spatial data in R
    [7][翻訳] RSelenium vignette: RSeleniumの基本
    [8]RからYahoo!のジオコーディングを利用する方法
    [9]EMBEDDING A LEAFLET MAP ON WORDPRESS
    [10]mapview advanced controls
    [11]RSeleniumでChromeからファイルをダウンロードするディレクトリを指定する方法
    [12]Selenium Serverが立ち上がらないときはportが被っているかも!?
    [13]brew install selenium-server-standalone
    [14]ナウでヤングなRの環境変数管理方法
    [15]タクシードライバー向け需要予測について
    [16]LDA with topicmodels package for R, how do I get the topic probability for each term?
    [17]dplyr — 高速data.frame処理