FacebookのMMMのOSS「Robyn」のチュートリアルさわってみた

目次
Robynとは
とりあえずチュートリアルやってみる
向き合い方
参考情報

Robynとは

Robyn(ロビン)はFacebook(META)が開発しているMarketing-Mix-Modeling(以降、MMM)のオープンソース(https://facebookexperimental.github.io/Robyn/)です。主にR言語で開発されています。(Python版は目下開発中らしいです。)

MMMは、マーケティングの広告投資の予算を、効果を最大化するためにどこに配分するかを数理的なモデルで決めようとするものです。
そもそも、MMMにそれほど明るくない方もおられると思いますが、その際は、『データ活用のための数理モデリング入門』の100ページ目に目的やら簡単な概要が載っていますので、参考にされると良いと思います。

また、RobynにはRobyn Open Source MMM Usersというユーザーグループがあるようです。そこそこ活発に運営されているようです。

私は以前、このコミュニティーのイベント(Facebook Project Robyn Open Source MMM 2021 Community Summit)があったので聴講しました。英語が苦手なので何度も聞きなおしましたが。。

この会では、世界中のマーケターがRobynを使ってプロモーションの分析をしているのを知れました。彼らはRobynのアーリーアダプターなんだろうなと思いました。

とりあえずチュートリアルやってみる

とりあえず、Robynがどういうツールなのか知るためにチュートリアル(https://facebookexperimental.github.io/Robyn/docs/quick-start/)をやってみることにしました。

まずは、Rを最新版にします。今回は「R version 4.2.0 (2022-04-22) — “Vigorous Calisthenics”」にしました。

結構な数の依存パッケージがインストールされているようです。スクロールバーがめちゃ小さくなりました。時間もそれなりにかかるようです。

Robynパッケージを読み込み、サンプルデータを呼び出します。

モデルを作る際の入力変数の設定をします。途中でProphetが使われているようです。確かに、Prophetは時系列の分析にちょうど適したライブラリではあります。

さて、指定する入力変数は結構多くて、以下の通りです。
・データセット
・従属変数
・従属変数のタイプ
・Prophetでの周期性などのオプション
・国
・競合の情報やイベントなどのコンテキスト変数
・ペイドメディアの支出
・ペイドメディアのインプレッションやクリックなど
・オーガニックな変数
・コンテキスト変数のなかでオーガニックなもの
・期初、期末
・広告の残存効果

実務でMMMを使っているものからすると馴染み深いですので、「ああ、このデータですね」とスッと頭に入ってきます。コードでいろいろ個別にやるととっ散らかるので、こういうカタマリで処理を実行できるのはいいですね。MMMが初めての方は、扱うデータセットの概要をよく調べてから入力変数にするようにしたほうがいいと思います。

次に、ハイパーパラメータの設定を行います。『StanとRでベイズ統計モデリング (Wonderful R)』の作法に従うならば、あまりハイパーパラメータを恣意的に決めて収束させるのはやりたくないですが、明らかに符号がおかしいとかの制約は付けてもいいのかなと思います。

ハイパーパラメータを設定したら、アルゴリズムを実行します。裏側でベイズ推定をしていることから、結構時間がかかります。Prophetを動かすということはStanを動かしていることと同義ですから。

イタレーションごとのモデルの目的関数の事後分布を可視化します。徐々に収束してそうに見えます。

モデルがいろいろと求まったので、パレート最適な組み合わせの計算をします。

パレート最適な組み合わせで返された複数のモデルから一つを選びます。

選んだモデルの係数などを確認します。

この選んだモデルをもとに、最適なアロケーションを計算します。

推定した、選んだモデルでの最適な広告のアロケーション結果を出力します。予算を削った方がいい広告経路、増やした方がいい広告経路などが示されます。

続いて、支出の上限を決めた上での、7日間でのアロケーションを行います。

こちらが、出力した結果です。

続いて、特定の広告経路の目的関数に対しての影響度が支出に応じてどう変わっていくか、つまりサチっているかどうかを見てみます。

支出に関して、目的関数がサチっているかどうかを見てみます。

新しいデータで、現在のモデルをアップデートします。

アップデートした場合、先ほどと同様に、推定結果や予算に応じたアロケーションを出力します。

向き合い方

以上、チュートリアルを行いましたが、過去にMMMを実務で使ったことがあるものとしては、Robynはかなりオートマチックなツールだなぁと思いました。時系列のベイズモデリングに対してProhpetに感じた感情と似ているかもしれません。
パレート最適なものを見つけたり、アロケーションをどうするかを決めたりする関数までもが用意されており、適切にモデルを作成することさえできれば、データサイエンティストの業務時間をかなり削減することができると思います。
ただ、残存効果をカスタマイズしたり、独自のモデルをやる自由度はある程度犠牲にしていると思うので、当てはまりにこだわる場合、これまで通りStanなどで独自にアルゴリズムを書くこともあってしかるべきかなと思います。

参考情報

https://facebookexperimental.github.io/Robyn/
データ活用のための数理モデリング入門
Robyn Open Source MMM Users
Facebook Project Robyn Open Source MMM 2021 Community Summit
https://facebookexperimental.github.io/Robyn/docs/quick-start/

ディリクレNBDモデルのマーケティング分野での適用に関して色々調べてみた


どの企業のマーケティング担当者も、自分が属する市場や自分の会社の商品の需要予測ないし購買行動の予測に関心があると思います。それらを予測できる手法として、ディリクレNBDモデルが候補にあがると思います。
以前、『確率思考の戦略論 USJでも実証された数学マーケティングの力』を読んだ際に、需要のシミュレーションにディリクレNBDモデルが使えるということを知りましたが、繰り返し利用するサービスや消費財向けのアプローチなのかなと思って、自分が扱うサービスへの適用可能性の低さから距離を置いてきました。
ところが、文献によると様々な分野でディリクレNBDモデルを使って生命保険やクレジットカードなどのサービス市場需要の予測がうまくいっているケースがあるらしく、真剣に向き合ってみようと思い調べてみることにしました。

目次

歴史
理論の概要
適用例
データについて
実践するコード
参考文献

歴史

マーケティングの分野では、消費者の繰り返し購買行動やダブルジョパディの法則など様々な現象を説明するための研究と実践がなされてきました。
そのために確率的なモデルが色々と適用されていきました。

“The Dirichlet model in marketing”という論文の図を引用しますが、1952年ごろからブランドロイヤルティーの研究が始まり、1958年には線形モデルが用いられ、1959年にはNBD(Negative Binomial Distribution:負の二項分布)を用いた購買予測などがなされるようになりました。
それらのアプローチを様々な対象に適用しつつ、1975年にはBBD(Beta Binomial Distribution:ベータ二項分布)が、ついに1984年にはNBD-DirichletモデルがGoodhardtらによって考案されるに至りました。

理論の概要

本家の論文(The Dirichlet: A Comprehensive Model of Buying Behaviour)や“Calculation of Theoretical Brand Performance Measures from the Parameters of the Dirichlet Model”を読まれることをお勧めしますが、この理論が想定しているものは以下の通りです。

・各々の消費者がポワソン過程に従って、負の二項分布に従い購買行動をする。
・消費者ごとの購買率はガンマ分布に従う。
・消費者の利用可能なブランドの選択は多項分布に従う。
・ブランドの選択確率自体は様々な消費者ごとに、多変量ベータ分布、あるいはディリクレ分布に従う。

なお、ディリクレモデルは負の二項分布とディリクレ多変量分布の二つの確率密度関数の結合からなり、二つの分布はそれぞれ独立であることを想定しています。
マーケットにおける購買自体の行動をつかさどる分布(負の二項分布)と、何を選ぶかをつかさどる分布(ディリクレ多変量分布)が独立しているということになります。

まず、負の二項分布についてですが、確率密度関数は以下の通りとなります。
$$ f_{\gamma, \beta}(k) = \frac{\Gamma (\gamma + k)}{\Gamma(\gamma)k!}\frac{\beta^k}{(1+\beta)^{\gamma + k}} \\ \textrm{ for } k = 0, 1, 2,…$$

負の二項分布は二つのパラメータを持ちます。
・形状(shape)パラメータ(\(\gamma\))
・尺度(scale)パラメータ(\(\beta\))

続いて、ディリクレ多変量分布ですが、確率密度関数は以下の通りとなります。

$$ f_{\alpha_1, \alpha_2, \dots, \alpha_h}(r_1, r_2, \dots, r_h | r_1 + r_2 + \dots + r_h = k ) \\
\frac{\Gamma \left ( \sum_{j=1}^{h} \alpha_j \right ) k! }{ \Gamma \left ( \sum_{j=1}^{h} \alpha_j + k \right ) } \prod_{j=1}^{h} \frac{\Gamma \left ( \alpha_j + r_j \right ) }{r_j! \Gamma \left ( \alpha_j \right )}$$

\(r_h\)などはそのマーケットにおける特定ブランドの購買率で、ブランド数がhだけあるとしています。
なお、
$$ r_1 + r_2 + \dots + r_h = k $$
で、足し合わせるとkがマーケット自体の購買率になるようになっています。ディリクレ多変量分布はブランドの数だけ、h個の正のパラメータ\( \alpha_1,\alpha_2, \dots, \alpha_h \)を持ちます。

ディリクレモデルは、以上の負の二項分布とディリクレ多変量分布を結合したもので、確率密度関数は以下の通りになります。

$$ f_{\gamma, \beta, \alpha_1, \alpha_2, \dots, \alpha_h}(r_1, r_2, \dots, r_h ) = \\ f_{\gamma, \beta}(k) f_{\alpha_1, \alpha_2, \dots, \alpha_h}(r_1, r_2, \dots, r_h | r_1 + r_2 + \dots + r_h = k )$$

この確率密度関数に色々かけて期待値を計算することで以下のようなものを算出することができます。
・あるブランドの購買者ごとの、あるブランドの購買の数の理論値
・プロダクトクラスでの購買の数の平均値
・あるブランドのみを買う母集団の割合
・購買者の平均購買頻度

諸々の式の導出は「[確率思考の戦略論] 1.確率理論の導入とプレファレンスの数学的説明」に詳しく書かれているので、見てみることをお勧めします。
数理統計学の知識とか、ベータ関数・ガンマ関数などの知識があれば式変形などに関してスムーズに理解できると思います。

ただ、この手法自体への批判などもないわけではなく、先ほどあげた論文、“The Dirichlet model in marketing”では、消費者行動を理解するのにつながるモデルではないとか、意思決定を扱っているモデルではないとする批判などがあるようです。ただ、幅広いマーケットで予測されてきた実績もあることから、モデルの限界を知りつつ適用するマーケットを増やしていくことが実務家に求められているのかなと思われます。

適用例

ディリクレモデルが様々な市場で適用されている事例が、『ディリクレモデルの境界条件 ― サービスへの適用可能性と限界 ―』という論文に記されていました。

当初は、インスタントコーヒー市場、衣料用洗剤市場、シャンプー市場などの最寄品(計画的に購入されることが少なく、単価は低く、何度も繰り返し購入される製品・サービス)のみで近似できるとされていましたが、
近年では銀行市場、クレジットカード市場、ガソリンスタンド市場、スーパーマーケット市場、テレビ番組市場、ファーストフード市場、オーケストラ市場、フットボールリーグ市場、野球市場などのサービス業や、スポーツウェアや自動車などの比較的高額な市場においても近似できるとされています。

データについて

ディリクレモデルの理論値を計算するには、4つのデータが必要とされています。
・1.そのカテゴリを購買した人々の全体に占める割合
・2.カテゴリに占める、いずれかの製品を購買した人々に対して記録された当該製品カテゴリの購買回数の平均
・3.各ブランドを一度でも購入した人々の割合
・4.各ブランドを購買した人々による各ブランドの購買回数の平均値

これらのデータはアンケートなどで収集する必要があるため、ディリクレモデルは総じて時間もお金もかかる手法であると言うこともできると思います。先ほどあげた、『ディリクレモデルの境界条件 ― サービスへの適用可能性と限界 ―』にはWEB調査でどんなアンケートを集めるのかの詳細まで書かれているので、自社の属するマーケットでこのモデルを適用する際は、それらを参考にすると良いと思います。

実践するコード

コードに関しては、こちらのColab(2_nb.ipynb)を書かれている方がおられたので、それを試すのも良いと思います。

ライブラリに関して、2016年と決して新しくないですが、R言語ではNBDdirichletというライブラリがあります。手法が1984年のものですから別に問題はないと思います。

こちらを実行した結果は、以下の通りです。

このdirichlet関数を使うことで、データを与えるだけで、ディリクレモデルのパラメータMやKやSが簡単に求まりました。

summaryで推定した諸々の値を確認できます。引数で色々と調整するので、実際に使う際はリファレンス(Package ‘NBDdirichlet’)を見てください。頻度のカットオフ値やヘビーユーザーの域値や、ブランド重複を見る際の基準ブランドなどを指定できます。

・buy:理論的なブランド浸透度、購買率、そのブランドの買い手によるマーケット内での購買率
・freq:そのブランドの購買の数の分布(以下の場合、6以上はまとめられている)
・heavy:特定の頻度の買い手におけるマーケットでの理論的なブランド浸透度、購買率
・dup:特定のブランドの購買者が他のブランドも購買している割合

この推定結果を使って、売上がどれくらいになるかが高精度にわかるのであれば、ブランドの浸透度が上がるような施策を打つことで、あるいは競合がそれらを打つことで売上がどう変わっていくかをシミュレーションすることもできるかもしれません。
その背景にはダブルジョパディの法則などがあると思うと結構面白そうですね。

参考文献

The Dirichlet: A Comprehensive Model of Buying Behaviour
Calculation of Theoretical Brand Performance Measures from the Parameters of the Dirichlet Model
The Dirichlet model in marketing
NBDdirichlet-package: NBD-Dirichlet model of consumer buying behavior
Package ‘NBDdirichlet’
ディリクレモデルの境界条件 ― サービスへの適用可能性と限界 ―
[確率思考の戦略論] 1.確率理論の導入とプレファレンスの数学的説明
Negative Binomial distribution
Dirichlet distribution

[Python]データ分析業務で使いそうなコードまとめ(随時更新)

はじめに

仕事で使いそうなPythonのコードを残しておくドキュメントが欲しいなと思ったので、よく使うものをこちらに貯めていこうと思います。まだ19個しかないですが、30個を目標に追記していきます。

フォーマットとしては、
1.やりたい処理
2.コード
3.参考情報のリンク
の3つを1セットにしています。
まずは、自分自身や周りで仕事をしている人が楽をできるドキュメントになればいいなと思って作っていきます。

目次

・重複削除
・階級のデータを作りたい
・再起的にリストをコピーしたい
・ピボットテーブルでアクセスログから特徴量用のデータを作りたい
・データフレームのある列に対して特定の文字列で分割し、その任意のN番目のものを抽出したい
・データフレームの行単位で割り算を行いたい
・列単位で行を足し合わせて文字列を作りたい
・グループのレコードごとにidを割り当てたい
・正規表現で任意の合致した文字列で挟まれた文字列を抽出したい
・複数列の集計値同士で除した値が欲しい
・等間隔の実数の数列を作成したい
・スペース区切りの文字列において、ある単語の前後N個の語を抽出したい
・グループごとに文字列を結合したい
・文字列の置換をdict形式のデータで行いたい
・テキストを形態素解析してBoW形式のデータにして特徴量にしたい
・訓練データとテストデータで特徴量の数が足りない時に不足した変数を追加したい
・アクセスログのユーザーごとのレコードごとの累積和を計算したい
・選択したカラムで一気にダミー変数にしたい
・グループ化して平均値を計算し、プロットし、2軸でそれぞれのデータ数も載せたい
・参考情報

データセット

全体を通じて使うデータセットを生成する関数を用意しておきます。一つ目の
データはkaggleのBoston Housingのデータセットです。

二つ目はアクセスログのデータが欲しかったので、Google Analytics APIを使ってデータを抽出したものとなります。ここではソースコードを載せますが、皆さんは当然見れません。一応、GitHubにデータを上げておきます。

三つ目はテキストデータです。

前処理

重複削除

  • やりたいこと
    重複したデータを除外したい。最初に現れたレコードだけを残したい。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.duplicated.html

階級のデータを作りたい

  • やりたいこと
    任意の範囲に応じたラベルをつけた階級値を作りたい。
  • コード
  • 参考情報
    https://pandas.pydata.org/docs/reference/api/pandas.cut.html

再起的にリストをコピーしたい

  • やりたいこと
    リストをコピーした際に、コピーしたものを変更した際に、もとのものも変更されてしまうので、それを防ぎたい。
  • コード
  • 参考情報
    https://docs.python.org/3/library/copy.html

ピボットテーブルでアクセスログから特徴量用のデータを作りたい

  • やりたいこと
    集計されていないアクセスログをもとに、セッションidごとの触れたページを集計し、Bag of Words形式のデータフレームにしたい。これは機械学習の特徴量に使えるので便利。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.pivot_table.html

データフレームのある列に対して特定の文字列で分割し、その任意のN番目のものを抽出したい

  • やりたいこと
    pandasのデータフレームである文字列の列に対して特定のキーワードで分割し、そのN番目のものをデータフレームの列として持たせる。今回はURLのパラメータの後ろについている文字列を抽出している。
  • コード
  • 参考情報

データフレームの行単位で割り算を行いたい

  • やりたいこと
    各行の合計で、各行を割ることでレコード単位で重み付けを行いたい時や確率を計算したいときに良く使います。データはアクセスログをBag of Word形式したものを使います。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.div.html

列単位で行を足し合わせて文字列を作りたい

  • やりたいこと
    データフレームの列の要素をつなぎ合わせて文字列を作成したい。複合キーを作らないとデータを繋げない時などによく使います。
  • コード
  • 参考情報
    https://www.kato-eng.info/entry/pandas-concat-pk

グループのレコードごとにidを割り当てたい

  • やりたいこと
    アクセスログなどで、任意のユーザーの初回に触れたページや、N番目に触れたページなどを抽出しやすいように、ユーザーごとのログにidを付与したい。
  • コード
  • 参考情報
    https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.factorize.html

正規表現で任意の合致した文字列で挟まれた文字列を抽出したい

  • やりたいこと
    正規表現を用いて、文字列に対して「指定した文字列と文字列の間にある表現」を抽出するための処理です。URLのデータや規則性のあるテキストデータに対して特徴量を作りたい時に使えます。
  • コード
  • 参考情報
    https://docs.python.org/ja/3/library/re.html

複数列の集計値同士で除した値が欲しい

  • やりたいこと
    複数の列で集計した値同士を用いて、さらにそこから比率のデータを作りたい。ターゲットエンコーディングを行いたい場合に使える。
  • コード
  • 参考情報
    https://stackoverflow.com/questions/45075626/summarize-using-multiple-columns-in-python-pandas-dataframe

等間隔の実数の数列を作成したい

  • やりたいこと
    パラメータチューニングをする際に、各パラメータについて等間隔の実数列を作りたいときがあるので、その時に使うコード。
  • コード
  • 参考情報
    https://www.pythoncentral.io/pythons-range-function-explained/

スペース区切りの文字列において、ある単語の前後N個の単語を抽出したい

  • やりたいこと
    テキストマイニングをする際に、関心のある単語の前後N文字を抽出したい時がある。もちろん、関心のある単語の周辺の単語というものを特徴量にすることも良いと思われる。ここではMeCabによる形態素解析により分かち書きにする関数と共に紹介する。
  • コード
  • 参考情報
    https://stackoverflow.com/questions/17645701/extract-words-surrounding-a-search-word

グループごとに文字列を結合したい

  • やりたいこと
    テキストマイニングなどで、任意のグループごとのテキストデータをスペース区切りで繋ぎたいときがあります。任意のグループごとのテキストを一気に連結できて便利です。これを形態素解析して特徴量にするのも良いと思います。
  • コード
  • 参考情報
    https://stackoverflow.com/questions/17841149/pandas-groupby-how-to-get-a-union-of-strings

文字列の置換をdict形式のデータで行いたい

  • やりたいこと
    事前に用意した置換のリストを使って、文字列の置換を一気に行いたい。
  • コード
  • 参考情報
    http://omoplatta.blogspot.com/2010/10/python_30.html

テキストを形態素解析してBoW形式のデータにして特徴量にしたい

  • やりたいこと
    任意の日本語のテキストとラベルの組み合わせがあるとして、テキストを形態素解析してTF-IDFを計算し、それを特徴量としてラベルを予測するということをやりたい。ここでは以前集めて、GitHubに載せてあるデザイナーズマンションのデータを使い、MeCabで形態素解析をしたあとに、TF-IDFを計算し、それを特徴量にしてsklearnで簡単に予測をし、リフト曲線を描いている。
  • コード
  • 参考情報
    https://scikit-plot.readthedocs.io/en/stable/metrics.html

訓練データとテストデータで特徴量の数が足りない時に不足した変数を追加したい

  • やりたいこと
    訓練データとテストデータの変数の数が違う時に、そのままではモデルを使えないので、補うために使う。
  • コード
  • 参考情報
    なし

アクセスログのユーザーごとのレコードごとの累積和を計算したい

  • やりたいこと
    アクセスログや購買データなどで、レコードごとの累積の値を計算したいときに使う。時系列に従いどんどん足されていく。
  • コード
  • 参考情報
    https://ohke.hateblo.jp/entry/2018/05/05/230000

選択したカラムで一気にダミー変数にしたい

  • やりたいこと
    ダミー変数を一気に生成したい。名前も一気に付けたい。
  • コード
  • 参考情報
    なし

グループ化して平均値を計算し、プロットし、2軸でそれぞれのデータ数も載せたい

  • やりたいこと
    EDAの時に、グループ化して平均値などを計算するだけでなく、該当するデータ数も2軸でプロットしたい。
  • コード
  • 参考情報
    なし

参考情報

前処理大全[データ分析のためのSQL/R/Python実践テクニック]
はじめてのアナリティクス API: サービス アカウント向け Python クイックスタート
カモノハシペリー
ペリー|フィニアスとファーブ|ディズニーキッズ公式

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


通勤電車のなかで私が勉強する用のシリーズ第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

[数理統計学]正規分布から導かれる分布(カイ二乗分布/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), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃

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

[数理統計学]連続型確率分布の期待値と分散の導出まとめ

はじめに

前回は、離散型確率分布に関するよくある分布とその平均と分散の導出についてひたすら記しました。今回は連続型の確率分布に関して同様に記していこうと思います。

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

今回登場する連続分布

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

1.一様分布

離散型の一様分布に対応しており、

$$ f(x; \alpha, \beta) = \frac{1}{\beta – \alpha}, \alpha \leq x \leq \beta $$
で定義されます。

ここで、ある区間(a,b)が区間(α,β)に含まれるとした場合、Xが一様分布に従う際の確率は、

$$ P(X \in (a,b)) = \int_a^b \frac{1}{\beta – \alpha} dx = \frac{b-a}{\beta – \alpha} $$
となり、確率は区間の長さ(b-a)に比例し、その区間が(α,β)のどこにあってもよい。つまり、一様にどこでも一定の確率となることを示している。

また、一様分布から生成される変数は、ある確率分布に従うような乱数を、確率分布関数の逆数から生成させることができる。そのような乱数生成方法を逆関数法と呼ぶ。

逆関数法
[0,1]上の一様分布に従う確率変数をUとする。このとき、目的の確率分布の関数の逆関数に一様乱数Uを代入したものは、目的の確率分布に従う。

Uが一様分布に従うことを式で示すと、
$$ P(U \leq u) = u, 0 \leq u \leq 1$$
ここで、
$$ U \leq u \Leftrightarrow F^{-1}(U) \leq F^{-1}(u) $$
であるから、
$$ P( F^{-1}(U) \leq F^{-1}(u) ) = u $$
となる。
さらに、\( F^{-1}(u) = x \)とおくと、
$$ P( F^{-1}(U) \leq x ) = F(x) $$
これは確率変数、\( F^{-1}(U) \)が累積分布関数が、\( F(x) \)であるような確率分布に従うことを表している。

一様分布の平均

$$ E(X) = \int_a^b xf(x) dx $$
$$ = \int_a^b x \frac{1}{b-a} dx$$
$$ = \left[ \frac{x^2}{2} \frac{1}{b-a} \right]^b_a $$
$$ = \frac{b^2 – a^2}{2(b-a)} $$
$$ = \frac{(b-a)(b+a)}{2(b-a)} $$
$$ = \frac{b+a}{2} $$

一様分布の分散

$$ E(X^2) = \int_a^b x^2f(x) dx $$
$$ = \int_a^b x^2 \frac{1}{b-a} dx $$
$$ = \left[ \frac{x^3}{3} \frac{1}{b-a} \right]^b_a $$
$$ = \frac{(b-a)(b^2 + ab + a^2)}{3(b-a)} $$
$$ = \frac{(b^2 + ab + a^2)}{3} $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ = \frac{(b^2 + ab + a^2)}{3} – \left( \frac{b+a}{2} \right)^2$$
$$ = \frac{(b^2 – 2ab + a^2)}{12} $$
$$ = \frac{(b – a)^2}{12} $$

2.正規分布

正規分布は以下のような式で定義される確率分布のことを指す。

$$ f(x; \mu, \sigma^2) = \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } $$
$$ -\infty < x < \infty , -\infty < \mu < \infty, \sigma > 0 $$

上の式で、μ=0、σ=1の場合は標準正規分布と呼ぶ。

実際に、正規分布に従う変数Xを標準化したもの\( Z = \frac{X-\mu}{\sigma} \)が標準正規分布に従うことを以下で示す。
$$ P(Z \leq z) = P(X \leq \mu + \sigma z) $$
$$ = \int_{-\infty}^{\mu + \sigma z} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx $$
ここで\( y = \frac{x-\mu}{\sigma} \)とすると、
$$ = \int_{-\infty}^{z} \frac{1}{ \sqrt{2 \pi} } e^{ – \frac{y^2}{2} }dy $$
となり、標準正規分布に従っていることがわかる。最後の積分の上限はもともとのxに対して、μを引き、σで割るという計算をすることで得られる。

また、正規分布の平均値や分散を計算する上で避けられない、ガウス積分について記しておく。

まず、
$$ I = \int_{-\infty}^{\infty} e^{-\frac{x^2}{2}} dx $$
とおき、xとyについての積分で

$$ I^2 = \left[ \int_{-\infty}^{\infty} e^{-\frac{x^2}{2}} dx \right] \left[ \int_{-\infty}^{\infty} e^{-\frac{y^2}{2}} dy \right] $$

$$ = \int_{-\infty}^{\infty} \!\!\! \int_{-\infty}^{\infty} e^{-\frac{x^2+y^2}{2}}dxdy $$
となるようにする。
ここで、重積分の変数変換を行う。\( x = r \cos \theta \)、\( y = r \sin \theta \)とし、dxdyをdθdrに変換するためにヤコビアンを用いて、
$$ = \int_{-\infty}^{\infty} \!\!\! \int_{-\infty}^{\infty} e^{-\frac{r^2}{2}} \left| \begin{array}{ccc} \frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\ \frac{\partial x}{\partial \theta} & \frac{\partial y}{\partial \theta} \end{array} \right| d\theta dr $$
と置き換える。

$$ = \int_{-\infty}^{\infty} \!\!\! \int_{-\infty}^{\infty} e^{-\frac{r^2}{2}} \left| \begin{array}{ccc} \cos \theta & \sin \theta \\ -r \sin \theta & r \cos \theta \end{array} \right| d\theta dr $$

$$ = \int_{-\infty}^{2\pi} \!\!\! \int_{0}^{\infty} e^{-\frac{r^2}{2}} r d\theta dr $$

$$ = \left[ \theta \int_{0}^{\infty} e^{-\frac{r^2}{2}} r dr \right]^{2\pi}_0$$
$$ = 2\pi \int_{0}^{\infty} e^{-\frac{r^2}{2}} r dr $$
$$ = 2\pi \left[ – e^{ -\frac{r^2}{2} } \right]^{\infty}_0$$
$$ = 2\pi $$
よって、
$$ I = \sqrt{2\pi} $$
となる。

さっそく、ガウス積分を使って、正規分布の積分が1になることを示す。

$$ \int_{-\infty}^{\infty} f(x) dx $$
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx $$
ここで、\( z = \frac{x-\mu}{\sigma}\)とおき、\( \frac{dx}{dz}=\sigma\)から、
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} } e^{ – \frac{z^2}{2} } dz $$
となる。
ガウス積分の公式(\( \int_{-\infty}^{\infty} e^{-\frac{z^2}{2}} dz = \sqrt{2\pi} \) )より、

$$ = \frac{1}{\sqrt{2 \pi}} \times \sqrt{2\pi} = 1 $$

正規分布の平均

$$ E(X) = \int_{-\infty}^{\infty} x \times \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx $$
ここで、\( z = \frac{x-\mu}{\sigma}\)とおき、\( \frac{dx}{dz}=\sigma\)から、
$$ = \int_{-\infty}^{\infty} (\mu + \sigma z) \times \frac{1}{ \sqrt{2 \pi} } e^{ – \frac{z^2}{2} } dz $$
第一項目は定積分が1になることから、第二項目は奇関数であることから定積分が0になるため、
$$ = \mu $$
となる。

正規分布の分散

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \int_{-\infty}^{\infty} x^2 \times \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } dx – \mu^2 $$
ここで、\( z = \frac{x-\mu}{\sigma}\)とおき、\( x = \sigma z + \mu \)、\( \frac{dx}{dz}=\sigma\)から、
$$ = \int_{-\infty}^{\infty} (\sigma z + \mu)^2 \times \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{z^2}{2} } dz – \mu^2 $$

$$ = \int_{-\infty}^{\infty} (\sigma^2 z^2 + 2\sigma \mu z + \mu^2) \times \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{z^2}{2} } dz – \mu^2 $$

ここで、第二項目は奇関数であることから定積分が0になるため、第三項目は定積分が1になることから、第三項目と第四項目が消し合うため、
$$ = \int_{-\infty}^{\infty} \sigma^2 z^2 \frac{1}{ \sqrt{2 \pi}} e^{ – \frac{z^2}{2} } dz $$
となる。
部分積分の公式より、

$$ = \frac{\sigma^2}{\sqrt{2 \pi}} \left( \left[ \frac{z^2}{-z} e^{ – \frac{z^2}{2} } \right]^\infty_-\infty – \int_{-\infty}^{\infty} 2z \times \frac{1}{-z} e^{ – \frac{z^2}{2}} \right) $$

第一項目はガウス積分に他ならないため、
$$ = \frac{\sigma^2}{\sqrt{2 \pi}} ( – \sqrt{2 \pi} + 2 \sqrt{2 \pi} ) $$
となり、
$$ = \sigma^2 $$
となる。

3.対数正規分布

正の値をとる確率変数Xを対数変換したものが正規分布に従うとした分布を対数正規分布と呼ぶ。

$$ f(x; \mu, \sigma^2) = \frac{1}{ \sqrt{2 \pi} \sigma x } e^{ – \frac{( \log x-\mu)^2}{2\sigma^2} } , x > 0$$

対数正規分布は以下のように正規分布から導出することができる。

\( f_x(x) = \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(x-\mu)^2}{2\sigma^2} } \)とし、\( x = g(y) = e^y \)とおき、以下の分布関数について以下のように変形する。

$$ F_x (x) = \int_{-\infty}^{x} f_x(t)dt $$
$$ = \int_{-\infty}^{x} f_y(g^{-1}(x)) \frac{dt}{dx} dx $$
$$ = \int_{-\infty}^{x} f_y(g^{-1}(x)) \frac{d g^{-1}(x) }{dx} dx $$

ここで両辺をxで微分すると、

$$ f_x(x) = f_y(g^{-1}x) \frac{d g^{-1}(x) }{dx} $$
$$ = f_y( \log x) \frac{1}{x} $$
$$ = \frac{1}{ \sqrt{2 \pi} \sigma x } e^{ – \frac{( \log x-\mu)^2}{2\sigma^2} } $$

対数正規分布の平均

$$ E(X) = E(e^y) $$
$$ = \int_{-\infty}^{\infty} e^y g(y) dy $$
$$ = \int_{-\infty}^{\infty} e^y \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(y-\mu)^2}{2\sigma^2} } dy $$
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{( y – \mu – \sigma^2 )^2}{2\sigma^2} + \mu + \frac{\sigma^2}{2} } dy $$
ここで、\( y – \mu – \sigma^2 = t \)とすると、
$$ = e^{\mu + \frac{\sigma^2}{2}} \frac{1}{ \sqrt{2 \pi} \sigma } \int_{-\infty}^{\infty} e^{-\frac{t^2}{2\sigma^2} } dt $$
となり、定積分の値はガウス積分の公式より\( \sqrt{2\pi} \sigma \)なので、
$$ = e^{\mu + \frac{\sigma^2}{2}} $$

対数正規分布の分散

$$ E(X^2) = E(e^{2y}) $$
$$ = \int_{-\infty}^{\infty} e^{2y} g(y) dy $$
$$ = \int_{-\infty}^{\infty} e^{2y} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{(y-\mu)^2}{2\sigma^2} } dy $$
$$ = \int_{-\infty}^{\infty} \frac{1}{ \sqrt{2 \pi} \sigma } e^{ – \frac{1}{2\sigma^2} \left( y – (\mu + 2\sigma^2 ) \right)^2 + 2\mu + 2\sigma^2 } dy $$
ここで、\( y – 2\sigma^2 = s \)とすると、
$$ = e^{2\mu + 2\sigma^2} \frac{1}{ \sqrt{2 \pi} \sigma } \int_{-\infty}^{\infty} e^{ – \frac{(s-\mu)^2}{2\sigma^2} } ds $$
ここで、\( s – \mu = t \)とおき、ガウス積分の公式を用いると、
$$ = e^{2\mu + 2\sigma^2} $$

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = e^{2\mu + 2\sigma^2} – e^{2\mu + \sigma^2} $$
$$ = e^{2\mu + \sigma^2} (e^{\sigma^2} – 1) $$

4.ガンマ分布

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

と表される確率密度関数をもつ分布をガンマ分布と呼ぶ。
\( \alpha, \beta \)の2つのパラメータを持ち、様々な形になり、正の値をとる変数の確率分布のモデルとして使われる。

ガンマ分布はポアソン過程から導出することができる。その際は、ガンマ関数が階乗の一般化であるという性質を使う。

まずガンマ関数が階乗の一般化であることの証明を行う。

まず、ガンマ関数は以下のように定義される。
$$ \Gamma (\alpha) = \int_{0}^{\infty} e^{-t} t^{\alpha – 1} dt $$
\( \alpha = 1 \)の場合、
$$ \Gamma (1) = \int_{0}^{\infty} e^{-t} dt $$
$$ = \left[ -e^{-t} \right]^{\infty}_{0} = 1 $$

任意の正の整数nに対して、
$$ \Gamma (n ) = \int_{0}^{\infty} e^{-t} t^{ n – 1} dt $$
部分積分の公式より、

$$ = \left[ -e^{-t} t^{n-1} \right]^{\infty}_{0} + (n-1)\int_{0}^{\infty} e^{-t}t^{n-2}dt$$

第1項目は0に、第2項目はガンマ関数の定義から以下のようになる。
$$ = (n-1) \Gamma (n-1) $$
以上より、逐次的に代入することで、
$$ \Gamma (n+1) = n! \Gamma (1) = n! $$
となる。

以下、ポアソン過程からガンマ分布を導出する。
まず、分布関数
$$ F_n(t) = P( Y_n \le t ) = P(X_t \ge n) $$
を考える。tや\( Y_n \)はn回の事象が生じるまでの時間を表し、\( X_t \)は事象の回数を表している。
\( F_n \)はポアソン分布に従うため、
$$ = \sum _{x=n}^{\infty} \frac{e^{-\lambda t} (\lambda t)^x }{x!} $$
と表され、全事象からn-1回分までの確率を差し引くことによっても表される。
$$ = 1 – \sum _{x=0}^{n-1} \frac{e^{-\lambda t} (\lambda t)^x }{x!} $$

\(Y_n \)の密度関数は分布関数をtで微分することによって得られる。

$$ \frac{d F_n (t)}{dt} = \sum _{x=0}^{n-1} \frac{\lambda e^{-\lambda t} (\lambda t)^x }{x!} – \sum _{x=1}^{n-1} \frac{\lambda e^{-\lambda t} \lambda^x t^{x-1} }{(x-1)!} $$
最後の項以外はキャンセルアウトされるため、
$$ = \frac{\lambda^n t^{n-1}e^{-\lambda t}}{(n-1)!} $$
となり、先程示したガンマ関数が階乗の一般系であることの証明から、
$$ = \frac{\lambda^n}{\Gamma (n) } t^{n-1}e^{-\lambda t} $$
となり、これはガンマ分布である。つまり、ポアソン過程に従う確率分布を時間で微分することでガンマ分布は得られる。

また、ガンマ分布には再生性があり、モーメント母関数を用いたり、確率分布の畳込みを用いての証明などがある。

ガンマ分布の平均

$$ E(X) = \int_{0}^{\infty} x \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha – 1}e^{-\beta x} dx $$
$$ = \int_{0}^{\infty} \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha}e^{-\beta x} dx $$
部分積分の公式より、

$$ \frac{ \beta^{\alpha} }{ \Gamma ( \alpha )} \left[ \frac{ x^{\alpha} e^{-\beta x} }{ -\beta } \right]^{\infty}_0 – \frac{ \beta^{\alpha} }{ \Gamma ( \alpha )} \int_{0}^{\infty} \frac{ \alpha x^{\alpha – 1}e^{-\beta x}}{-\beta} dx $$

第1項目は0になるため、
$$ = \frac{\alpha}{\beta} \int_{0}^{\infty} \frac{\beta^{\alpha} x^{\alpha-1}e^{-\beta x}}{\Gamma (\alpha) } $$
となり、ガンマ分布の定積分は1であることから、
$$ = \frac{\alpha}{\beta} $$
となる。

ガンマ分布の分散

$$ E(X^2) = \int_{0}^{\infty} x^{2} \frac{\beta^{\alpha} }{ \Gamma (\alpha) } x^{\alpha – 1}e^{-\beta x} dx $$
$$ = \frac{\beta^{\alpha}}{\Gamma (\alpha) } \int_{0}^{\infty} x^{\alpha + 1}e^{-\beta x} dx $$
部分積分の公式より、

$$ = \frac{\beta^{\alpha}}{\Gamma (\alpha) } \left[ \frac{x^{\alpha + 1}e^{-\beta x}}{- \beta } \right]^{\infty}_0 $$

$$ – \frac{\beta^{\alpha}}{\Gamma (\alpha) } \int_{0}^{\infty} \frac{(\alpha + 1)}{- \beta} x^{\alpha}e^{-\beta x} dx $$

第1項目は0になるため、
$$ = \frac{ (\alpha + 1) \beta^{\alpha-1}}{\Gamma (\alpha) } \int_{0}^{\infty} x^{\alpha}e^{-\beta x} dx $$
となり、ガンマ関数の性質より、\( \Gamma (\alpha + 1) = \alpha \Gamma (\alpha) \)から、
$$ = \alpha \beta^{-2} (\alpha + 1) \int_{0}^{\infty} \frac{\beta^{\alpha + 1}}{\Gamma (\alpha + 1)}x^{\alpha} e^{-\beta x} dx $$
ガンマ分布の定積分が1であることから、
$$ = \alpha \beta^{-2} (\alpha + 1) $$
となる。

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \alpha \beta^{-2} (\alpha + 1) – \frac{\alpha^{2}}{\beta^{2}} $$
$$ = \frac{\alpha}{\beta^2} $$

5.指数分布

ガンマ分布の特別なケースとして、\( \alpha = 1, \beta = \frac{1}{\theta} \) とした、
$$ f(x; \theta) = \frac{1}{\theta} e^{-\frac{1}{\theta}x} , x > 0, \theta > 0 $$
の形をした分布を指数分布と呼ぶ。これはポアソン過程における\( n=1 \)、つまり初めて事象が起こるまでの時間の分布を表している。加えて、無記憶性があるため、幾何分布の連続型版と考えることもできる。

無記憶性について

$$ P(X > x + h | X > x) = \frac{\int^{\infty}_{x+h} \frac{1}{\theta} e^{-\frac{y}{\theta}} dy }{ \int^{\infty}_{x} \frac{1}{\theta} e^{-\frac{y}{\theta}} dy } $$

$$ = \frac{ \left[ \left( – \frac{1}{\theta} \right) e^{-\frac{y}{\theta}} \right]^{\infty}_{x+h} }{ \left[ \left( – \frac{1}{\theta} \right) e^{-\frac{y}{\theta}} \right]^{\infty}_{x} } $$
$$ = \frac{-e^{-\frac{x+h}{\theta}}}{-e^{-\frac{x}{\theta}}} = e^{-\frac{h}{\theta}} $$
よって、確率密度がこれまでの時間には依存しないことがわかる。

ハザードレートについて

あるt時点での確率密度を、t時点までの確率分布を1から引いたもので割ったものをハザードレートとする。(事象が起きるまでの確率1単位あたりの確率密度)
$$ \frac{f(t; \theta)}{1 – F(t;\theta)} = \frac{ \frac{1}{\theta} e^{-\frac{1}{\theta} t} }{ \int^{\infty}_{t} \frac{1}{\theta} e^{-\frac{1}{\theta} x} dx } $$
$$ = \frac{ \frac{1}{\theta} e^{-\frac{1}{\theta} t} }{ \left[ -e^{\frac{1}{\theta}x} \right]^{\infty}_{t} } = \frac{1}{\theta} $$
よって、tによらず一定であることがわかる。これは製品の寿命の分布と考えるに際して、時間を通じて寿命をむかえる確率が一定という指数分布の性質はクレイジーなものと思われますね。

ここで、指数分布の平均や分散を計算する上で必要な極限について取り上げておきます。

\( xe^{-x} \)の極限について

\( x > 0 \)として、
$$ x < 2^{x} $$
とし、両辺に\( e^{-x} \)をかけると、
$$ x e^{-x} < 2^{x} e^{-x} $$

ここで右辺に関して極限を取ると、
$$ \lim_{x \to \infty} 2^{x} e^{-x} = \lim_{x \to \infty} \left( \frac{2}{e} \right)^x = 0 $$
となる。
よって、
$$ \lim_{x \to \infty} x e^{-x} = 0 $$
となる。

\( x^n e^{-x} \)の極限について

\( x > 0 \)として、
$$ f(x) = x^{n+1}e^{-x} $$
を考える。
これを微分すると、
$$ f'(x) = (n+1)x^{n}e^{-x} – x^{n+1}e^{-x} $$
$$ = x^ne^{-x}(n+1-x) $$
ここで、\( f'(x) = 0 \)とすると、
$$ x = n+1 $$
となり、最大値は
$$ f(n+1) = (n+1)^{n+1} e^{-(n+1)} $$
となる。そこで以下の不等式から

$$ 0 < x^ne^{-x} = x^{n+1}e^{-x}x^{-1} \le (n+1)^{n+1}e^{-(n+1)}x^{-x} $$

最後の式について極限を取ると、
$$ \lim_{x \to \infty} (n+1)^{n+1}e^{-(n+1)}x^{-x} = 0 $$
となり、はさみうちされるため、
$$ \lim_{x \to \infty} x^n e^{-x} = 0$$
となる。

指数分布の平均

$$ E(X) = \int_{0}^{\infty} x \frac{1}{\theta} e^{ – \frac{1}{ \theta } x } dx $$
部分積分の公式より、
$$ = \left[ x \frac{-\theta}{\theta} e^{- \frac{1}{\theta}x} \right]^{\infty}_{0} -\int_{0}^{\infty} \frac{1}{\theta} (-\theta) e^{-\frac{1}{\theta}x} dx $$
\( xe^{-x} \)の極限より、第1項目は0になるため、
$$ = \int_{0}^{\infty} e^{-\frac{1}{\theta}x}dx $$
また指数分布の定積分が1であることから、
$$ = \theta \int_{0}^{\infty} \frac{1}{\theta} e^{-\frac{1}{\theta}x}dx = \theta $$
となる。

指数分布の分散

$$ E(X^2) = \int_{0}^{\infty} x^2 \frac{1}{\theta} e^{ – \frac{1}{ \theta } x } dx $$
部分積分の公式より、
$$ = \left[ x^2 \frac{-\theta}{\theta} e^{- \frac{1}{\theta}x} \right]^{\infty}_{0} – \int_{0}^{\infty} \frac{1}{\theta} (-\theta) 2x e^{-\frac{1}{\theta}x} dx $$

\( x^n e^{-x} \)の極限より、第1項目は0になるため、
$$ = \int_{0}^{\infty} 2x e^{-\frac{1}{\theta}x} dx $$
再び、部分積分の公式より、
$$ = 2 \left[ x( – \theta )e^{ – \frac{1}{ \theta } x } \right]^{\infty}_{0} – 2 \int_{0}^{\infty} ( – \theta )e^{-\frac{1}{\theta}x}dx $$

\( x e^{-x} \)の極限より、第1項目は0になるため、
$$ = 2 \int_{0}^{\infty} (\theta)e^{-\frac{1}{\theta}x}dx $$
$$ = 2 \theta^2 \int_{0}^{\infty} \frac{1}{\theta} e^{-\frac{1}{\theta}x}dx $$
指数分布の定積分が1であることから、
$$ = 2\theta^2 $$

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = 2\theta^2 – \theta^2 = \theta^2 $$

6.ワイブル分布

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

の確率密度に従う分布をワイブル分布と呼び、\( \alpha = 1 \)のときは指数分布になる。

ハザードレートについて

あるt時点での確率密度を、t時点までの確率分布を1から引いたもので割ったものをハザードレートとする。(事象が起きるまでの確率1単位あたりの確率密度)

$$ \frac{f(t; \alpha, \beta)}{1 – F(t;\alpha, \beta)} = \frac{\beta \alpha t^{\alpha – 1} e^{-\beta t^{\alpha}}}{\int_{t}^{\infty}\beta \alpha x^{\alpha – 1}e^{-\beta x^{\alpha}}dx} $$

ここで、分母についてのみ注目して、\( z = x^{\alpha} \)と変数変換すると、\( \frac{dx}{dz} = \frac{1}{\alpha} z^{\frac{1}{\alpha}-1} \)から、

$$ \int_{t}^{\infty}\beta \alpha x^{\alpha – 1}e^{-\beta x^{\alpha}}dx = \beta \int_{t^{\alpha}}^{\infty} e^{-\beta z}dz $$

$$ = \beta \left[ – \frac{1}{\beta} e^{-\beta z} \right]^{\infty}_{t^{\alpha}} = e^{\beta t^{\alpha}} $$
よってハザードレートは、
$$ \frac{f(t; \alpha, \beta)}{1 – F(t;\alpha, \beta)} = \beta \alpha t^{\alpha -1} $$
となる。これは時間とともにハザードレートが高くなるという性質を示しており、全ての対象に対して、寿命に関して時間に対して単調増加を想定するのは誤りだろう。

ワイブル分布の平均

$$ E(X) = \int_{-\infty}^{\infty} x \beta \alpha x^{\alpha-1}e^{-\beta x^{\alpha}} dx $$
\( z = \beta x^{\alpha} \)として変数変換すると、\( \frac{dx}{dz} = \frac{1}{\alpha} z^{\frac{1}{\alpha}-1} \beta ^{-\frac{1}{\alpha}} \)から、
$$ = \int_{-\infty}^{\infty} \left( \frac{z}{\beta} \right)^{\frac{1}{\alpha}} \beta \alpha \left( \frac{z}{\beta} \right)^{\frac{\alpha – 1}{\alpha}} e^{-z} \frac{1}{\alpha}z^{\frac{1}{\alpha}-1} \beta ^{-\frac{1}{\alpha}} dz $$
$$ = \beta^{-\frac{1}{\alpha}} \int_{-\infty}^{\infty} z^{\frac{1}{\alpha}} e^{-z}dz $$
ガンマ関数の定義より、
$$ = \beta^{-\frac{1}{\alpha}} \Gamma \left( \frac{1}{\alpha} + 1 \right) $$

ワイブル分布の分散

$$ E(X^2) = \int_{-\infty}^{\infty} x^2 \beta \alpha x^{\alpha-1}e^{-\beta x^{\alpha}} dx $$
\( z = \beta x^{\alpha} \)として変数変換すると、
$$ = \int_{-\infty}^{\infty} \left( \frac{z}{\beta} \right)^{\frac{2}{\alpha}} \beta \alpha \left( \frac{z}{\beta} \right)^{\frac{\alpha – 1}{\alpha}} e^{-z} \frac{1}{\alpha}z^{\frac{1}{\alpha}-1} \beta ^{-\frac{1}{\alpha}} dz $$
$$ = \beta^{-\frac{2}{\alpha}} \int_{-\infty}^{\infty} z^{\frac{2}{\alpha}} e^{-z}dz $$
ガンマ関数の定義より、
$$ = \beta^{-\frac{2}{\alpha}} \Gamma \left( \frac{2}{\alpha} + 1 \right) $$
となる。

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \beta^{-\frac{2}{\alpha}} \Gamma \left( \frac{2}{\alpha} + 1 \right) – \left[ \beta^{-\frac{1}{\alpha}} \Gamma \left( \frac{1}{\alpha} + 1 \right) \right]^2 $$
$$ = \beta^{-\frac{2}{\alpha}} \left[ \Gamma \left( \frac{2}{\alpha} + 1 \right) – \Gamma^2 \left( \frac{1}{\alpha} + 1 \right) \right] $$

7.ベータ分布

$$ f(x;\alpha, \beta) = \frac{1}{B(\alpha, \beta)}x^{\alpha – 1}(1-x)^{\beta – 1} , 0 < x < 1, \alpha > 0, \beta > 0 $$

ただし、\( B(\alpha, \beta) = \int_{0}^{1} x^{\alpha – 1}(1-x)^{\beta – 1} dx \)とする。この確率密度に従う分布をベータ分布と呼び、(0, 1)上の確率現象などをモデリングする際に用いられる。

ベータ関数の復習

ベータ関数は任意の2つの正定数x,yに対して定義される2変数関数で、
$$ B(x, y) = \int_{0}^{1} t^{x – 1}(1-t)^{y – 1} dt $$
で表され、収束する。
証明は区間を\( (0,\frac{1}{2}] \)と\( [ \frac{1}{2}, 1) \)に分けた積分について、各項が収束することを示すことによりなされる。

・ベータ関数の非負性、対称性

$$ B(x,y) > 0 $$

$$ B(x,y) = B(y,x) $$

対称性の証明は\( s = 1 – t \)と変数変換することで容易に示される。

・ベータ関数を三角関数の積分に変数変換

$$ B(x, y) = 2 \int_{0}^{\frac{\pi}{2}} \sin ^{2x-1} \theta \cos ^{2y – 1} \theta d \theta $$

証明は、\( t = \sin^2 \theta \)として変数変換をすることで容易に示される。

・ベータ関数とガンマ関数の関係について

$$ B(x, y) = \frac{ \Gamma (x) \Gamma (y) }{ \Gamma (x+y) } $$

証明はガンマ関数の定義式について、\( t = s^2 \)と変数変換を行い、極座標の変換公式を用いてから、三角関数で表現できるベータ関数を用いることで示される。

ベータ分布の平均

$$ E(X) = \int_{0}^{1} x \frac{1}{B(\alpha, \beta)}x^{\alpha – 1}(1-x)^{\beta – 1} dx $$

ベータ関数とガンマ関数の関係、\( B(x, y) = \frac{ \Gamma (x) \Gamma (y) }{ \Gamma (x+y) } \)から、
$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + \beta ) }{ \Gamma (\alpha) \Gamma (\beta) } x^{\alpha}(1-x)^{\beta – 1} dx $$
ガンマ関数の性質\( \Gamma ( x + 1 ) = x \Gamma (x) \)より、

$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + 1 + \beta ) }{ \alpha + \beta } \frac{\alpha}{\Gamma (\alpha + 1) \Gamma (\beta)} x^{\alpha}(1-x)^{\beta – 1} dx $$

$$ = \frac{\alpha}{\alpha + \beta} \int_{0}^{1} \frac{ \Gamma (\alpha + 1 + \beta ) }{ \Gamma (\alpha + 1) \Gamma (\beta) } x^{\alpha}(1-x)^{\beta – 1} dx $$
ベータ分布の定義より、定積分が1になることから、
$$ = \frac{\alpha}{\alpha + \beta} $$

ベータ分布の分散

$$ E(X^2) = \int_{0}^{1} x^2 \frac{1}{B(\alpha, \beta)}x^{\alpha – 1}(1-x)^{\beta – 1} dx $$
$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + \beta ) }{ \Gamma (\alpha) \Gamma (\beta) } x^{\alpha+1}(1-x)^{\beta – 1} dx $$
ガンマ関数の性質より、

$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + 1 + \beta ) }{ \alpha + \beta } \frac{\alpha}{\Gamma (\alpha + 1) \Gamma (\beta)} x^{\alpha+1}(1-x)^{\beta – 1} dx $$

さらに、ガンマ関数の性質より、

$$ = \int_{0}^{1} \frac{ \Gamma (\alpha + 2 + \beta ) }{ (\alpha + \beta) (\alpha + \beta + 1) } \frac{\alpha(\alpha + 1)}{\Gamma (\alpha + 2) \Gamma (\beta)} x^{\alpha+1}(1-x)^{\beta – 1} dx $$

整理すると、

$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} \int_{0}^{1} \frac{ \Gamma (\alpha + 2 + \beta ) }{\Gamma (\alpha + 2) \Gamma (\beta)} x^{\alpha+1}(1-x)^{\beta – 1} dx $$

となる。ベータ分布の定義より、定積分が1になることから、
$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} $$
となる。

$$ V(X) = E(X^2) – \left[ E(X) \right]^2 $$
$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} – \left[ \frac{\alpha}{\alpha + \beta} \right]^2 $$
$$ = \frac{\alpha(\alpha + 1)}{(\alpha + \beta) (\alpha + \beta + 1)} – \frac{\alpha^2}{(\alpha + \beta)^2} $$
$$ = \frac{\alpha(\alpha+1)(\alpha+\beta) – \alpha^2(\alpha + \beta +1)}{(\alpha + \beta)^2(\alpha + \beta +1)} $$
$$ = \frac{\alpha \beta}{(\alpha + \beta)^2(\alpha + \beta +1)} $$

8.コーシー分布

$$ f(x;\theta) = \frac{1}{\pi} \frac{1}{1+(x-\theta)^2} , – \infty < x < \infty, – \infty < \theta < \infty $$

という確率密度に従う分布をコーシー分布と呼ぶ。裾が長く、大きい値や小さい値を取る確率がなかなか0に近づかない分布とされる。
また、コーシー分布の特徴としては平均や分散を持たないことがあげられる。身近な活用例としては、株価の分析に使われたり、ベイズ統計学における無情報事前分布として半コーシー分布というものが用いられることがある。

ここで、\( x – \theta = z \)として、コーシー分布の平均を考える。
$$ E(X) = \int_{-\infty}^{\infty} z \frac{1}{\pi} \frac{1}{1+z^2} dz $$
上限をb、下限をaとすると、
$$ = \frac{1}{\pi} \int_{a}^{b} \frac{z}{1+z^2} dz$$
$$ = \frac{1}{\pi} \left[ \frac{1}{2} \log (1+z^2) \right]^b_a $$
$$ = \frac{1}{\pi 2} \log (1+b^2) – \frac{1}{\pi 2} \log (1+a^2) $$
ここで、aとbについて極限をとると、

$$ \lim_{b \to \infty} \lim_{a \to \infty} \left( \frac{1}{\pi 2} \log (1+b^2) – \frac{1}{\pi 2} \log (1+a^2) \right) $$

となり、これは不定となり極限値は存在しない。よって、平均値は存在しないため、分散も同様に存在しない。各々の項が無限にランダムに向かうため特定の値に収束することはないというイメージをするとわかりやすい。

9.多項分布

$$P(X_1=x_1, X_2=x_2, \dots, X_k=x_k) \\
:= \frac{n!}{x_1! x_2! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} \\
x_1, x_2, \dots x_k = 0, 1, \dots , n \\
x_1 + x_2 + \dots + x_k = n
$$

\( X_i \)は事象\( A_i \)の起こった回数を表し、\( x_i \)はその実現値を表す。事象\( A_i \)は標本空間において試行した結果として生じる。その試行は独立でn回繰り返しなされる。\( \theta_i = P( A_i ) \)である。
多項分布の期待値と分散の計算のために周辺結合分布を用いる必要があるため、ここでは多項分布の周辺化を行う。

\( X_i \)と\( X_j \)の周辺化

$$ P(X_i=x_i, X_j=x_j) \\
= \sum _{ \underset{ < i,j > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

\( x_i \)と\( x_j \)の項はここでは合計の対象として除かれるので、以下のように書き表すことができる。ただし後半の表現が一見すると何も変化がないように見えるため注意されたい。

$$ = \frac{ \theta_i^{x_i} \theta_j^{x_j} }{ x_i!x_j! } \sum _{ \underset{ < i,j > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

$$ = \frac{ \theta_i^{x_i} \theta_j^{x_j} (1 – \theta_i – \theta_j)^{n- x_i – x_j} n! }{ (n – x_i – x_j)! x_i! x_j! } \sum _{ \underset{ < i,j > } { x_1,\dots , x_k } } \frac{(n – x_i – x_j )!}{x_1! \dots x_k!} \left ( \frac{\theta_1}{1-\theta_i-\theta_j} \right )^{x_1} \left ( \frac{\theta_2}{1-\theta_i-\theta_j} \right )^{x_2} \dots \left ( \frac{\theta_k}{1-\theta_i-\theta_j} \right )^{x_k} $$

後半の部分はiとjを除いた多項分布のため、1となることから、
$$ = \frac{ n!}{ (n – x_i – x_j)! x_i! x_j! } \theta_i^{x_i} \theta_j^{x_j} (1 – \theta_i – \theta_j)^{n- x_i – x_j} $$

これは3項分布に他ならない。つまり、多項分布の2つの変数の周辺分布は3項分布となる。

\( X_i \)の周辺化

先程と同様にして周辺化を行う。
$$ P(X_i=x_i) = \sum _{ \underset{ < i > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

$$ = \frac{ \theta_i^{x_i} }{ x_i!} \sum _{ \underset{ < i > } { x_1,\dots , x_k } } \frac{n!}{x_1! \dots x_k!} \theta_1^{x_1} \theta_2^{x_2} \dots \theta_k^{x_k} $$

$$ = \frac{ \theta_i^{x_i} (1 – \theta_i )^{n- x_i } n! }{ (n – x_i)! x_i!} \sum _{ \underset{ < i > } { x_1,\dots , x_k } } \frac{(n – x_i )!}{x_1! \dots x_k!} \left ( \frac{\theta_1}{1-\theta_i} \right )^{x_1} \left ( \frac{\theta_2}{1-\theta_i} \right )^{x_2} \dots \left ( \frac{\theta_k}{1-\theta_i} \right )^{x_k} $$

$$ = \frac{ n! }{ (n – x_i)! x_i!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } \\
= {}_n C _{x_i} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$

これは二項分布であることから、多項分布の一つの変数についての周辺分布は2項分布となる。

多項分布の平均

多項分布の期待値は二項分布のものと同じになる。

$$ E(X_i) = \sum _{ x_i \geq 0 , x_i \leq n } x_i \frac{ n! }{ (n – x_i)! x_i!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$

$$ = \sum _{ x_i \geq 0 , x_i \leq n } x_i \frac{ n! }{ (x_i – 1)! \\\{(n-1) – (x_i – 1)\\\}! } \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$
$$ = n \theta_i \sum _{ x_i \geq 1 , x_i \leq n } x_i \frac{ (n-1)! }{ (x_i – 1)! \\{(n-1) – (x_i – 1)\\}! } \theta_i^{x_i-1} (1 – \theta_i )^{(n-1)- (x_i-1) } $$

$$ = n \theta_i $$

多項分布の分散

多項分布の分散は二項分布のものと同じになる。

$$ V(X) = E(X^2) – \left [ E(X) \right ]^2 \\
= E(X(X-1)) + E(X) – \left [ E(X) \right ]^2
$$

$$ E(X_i(X_i-1)) \\
= \sum _{ x_i \geq 0 , x_i \leq n } x_i(x_i – 1) \frac{ n! }{ (n – x_i)! x_i!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$

$$ = \sum _{ x_i \geq 0 , x_i \leq n } \frac{ n! }{ \\{(n – 2) – (x_i -2 ) \\}! (x_i – 2)!} \theta_i^{x_i} (1 – \theta_i )^{n- x_i } $$
$$ = n(n-1) \theta_i^2 \\\
\sum _{ x_i \geq 2 , x_i \leq n } x_i \frac{ (n-2)! }{ (x_i – 2)! \\{(n-2) – (x_i – 2)\\}! } \theta_i^{x_i-2} (1 – \theta_i )^{(n-2)- (x_i-2) } $$

$$= n(n-1) \theta_i^2$$

$$ V(X) = E(X(X-1)) + E(X) – \left [ E(X) \right ]^2 \\
= n(n-1) \theta_i^2 + n \theta_i – ( n \theta_i )^2 \\
= n \theta_i (1 – \theta_i)
$$

多項分布の共分散

$$
E(X_i X_j) = \sum _{ x_i,x_j \geq 0 \\\ x_i + x_j \leq n } x_i x_j \frac{n!}{x_i! x_j!(n-x_i-x_j)!} \theta_i^{x_i}\theta_j^{x_j}(1-\theta_i – \theta_j)^{n-x_i-x_j}
$$
$$ = \sum _{ x_i,x_j \geq 0 \\\ x_i + x_j \leq n } \frac{n!}{(x_i-1)! (x_j-1)![(n-2)-(x_i-1)-(x_j-1)]!} \theta_i^{x_i}\theta_j^{x_j}(1-\theta_i – \theta_j)^{n-x_i-x_j} $$
$$ = n(n-1) \theta_i \theta_j \sum _{ x_i,x_j \geq 1 \\\ x_i + x_j \leq n } \frac{(n-2)!}{(x_i-1)! (x_j-1)![(n-2)-(x_i-1)-(x_j-1)]!} \theta_i^{x_i-1}\theta_j^{x_j-1}(1-\theta_i – \theta_j)^{(n-2)-(x_i-1)-(x_j-1)} $$

$$ = n(n-1)\theta_i \theta_j $$

$$ Cov(X_i, X_j) = E( X_i X_j ) – E(X_i) E(X_j) \\
= n(n-1)\theta_i \theta_j – n\theta_i n\theta_j \\
= – n \theta_i \theta_j
$$

多項分布の条件付き分布

\( X_j = x_j \)のときの\( X_i \)の条件付き分布は\( P(X_i = x_i | X_j = x_j) = \frac{P(X_i \cap X_j)}{P(X_j)} \)より、

$$ P(X_i = x_i | X_j = x_j) = \frac{ \frac{n!}{x_i!x_j!(n-x_i-x_j)!} \theta_i^{x_i} \theta_j^{x_j} (1-\theta_i – \theta_j)^{n-x_i-x_j} }{ \frac{n!}{x_j!(n-x_j)!} \theta_j^{x_j} (1- \theta_j)^{n-x_j} } $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \theta_i^{x_i} \left ( \frac{1}{1-\theta_j} \right )^{x_j – n} \left ( 1 – \theta_i – \theta_j \right )^{n-x_i-x_j} $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \left ( \frac{\theta_i}{1-\theta_j} \right )^{x_i} \left ( \frac{1}{1-\theta_j} \right )^{n- x_i -x_j} \left ( 1 – \theta_i – \theta_j \right )^{n-x_i-x_j} $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \left ( \frac{\theta_i}{1-\theta_j} \right )^{x_i} \left ( \frac{ 1 – \theta_i – \theta_j }{1-\theta_j} \right )^{n- x_i -x_j} $$
$$ = \frac{(n-x_j)!}{x_i!(n-x_i-x_j)!} \left ( \frac{\theta_i}{1-\theta_j} \right )^{x_i} \left ( 1 – \frac{\theta_i}{1-\theta_j} \right )^{n- x_i -x_j} $$

$$ = {}_{n-x_j} C _{x_i} \hat \theta_i ^{x_i}(1 – \hat \theta_i)^{(n-x_j)-x_i} $$

$$ \hat \theta_i = \frac{\theta_i}{1-\theta_j} $$

よってこれは二項分布であるから、多項分布の\( X_j = x_j \)のときの\( X_i \)の条件付き分布は二項分布に従うことがわかる。

10.ディリクレ分布

ディリクレ分布は\( \sum_{i=1}^{m} \theta_i = 1, \theta_i \geq 0 \)を満たす\( \theta_1, \dots , \theta_m \)の集合からなる。
ディリクレ分布の確率密度関数は以下の形をとる。

$$ f(\theta_1, \dots , \theta_m) = \frac{\Gamma (A)}{ \prod_{i=1}^{m} \Gamma (\alpha_i) } \prod_{i=1}^{m} \theta_i^{\alpha_i – 1} \\
A = \sum_{i=1}^{m} \alpha_i \\
\alpha_i > 0, \quad i = 1, \dots , m
$$

\( m = 2 \)のとき、ディリクレ分布はベータ分布の特殊型になる。

$$ f(\theta_1, \theta_2) = \frac{ \Gamma ( \alpha_1 + \alpha_2 ) }{ \Gamma (\alpha_1) \Gamma (\alpha_2) } \theta_1^{\alpha_1 – 1} \theta_2^{\alpha_2 – 1} $$

ディリクレ分布の平均

$$ E(\theta_j) = \int \dots \int \theta_j \frac{\Gamma (A)}{ \prod_{i=1}^{m} \Gamma (\alpha_i) } \prod_{i=1}^{m} \theta_i^{\alpha_i – 1} d\theta_1 \dots d\theta_{m-1} $$
$$ = \frac{\Gamma(A)}{\Gamma(A+1)} \frac{\Gamma(\alpha_j + 1)}{\Gamma(\alpha_j)} \int \dots \int \frac{\Gamma (A+1)}{ \prod_{i=1}^{m} \Gamma (\alpha_i^{\prime}) } \prod_{i=1}^{m} \theta_i^{\alpha_i^{\prime} – 1} d\theta_1 \dots d\theta_{m-1} $$

ここでは\(
\begin{cases}
\alpha_i^{\prime} = \alpha_i \quad i \neq j \\
\alpha_j^{\prime} = \alpha_j + 1 \quad i = j \\
\end{cases}
\)となる。

積分以降はディリクレ分布の定義より1なので、以下のように表される。

$$ = \frac{\Gamma(A)}{\Gamma(A+1)} \frac{\Gamma(\alpha_j + 1)}{\Gamma(\alpha_j)} $$
ガンマ関数の定義より、

$$ = \frac{\Gamma(A)}{\Gamma(A)A } \frac{\Gamma(\alpha_j)\alpha_j}{\Gamma(\alpha_j)} $$
となることから、
$$ = \frac{\alpha_j}{A} $$
となる。

ディリクレ分布の分散

分散は以下の式より求まる。
$$ V(\theta_j) = E \left( \theta_j^2 \right) – \left [ E(\theta_j) \right ]^2 $$

$$ E(\theta_j^2) = \int \dots \int \theta_j^2 \frac{\Gamma (A)}{ \prod_{i=1}^{m} \Gamma (\alpha_i) } \prod_{i=1}^{m} \theta_i^{\alpha_i – 1} d\theta_1 \dots d\theta_{m-1} $$
$$ = \frac{\Gamma(A)}{\Gamma(A+2)} \frac{\Gamma(\alpha_j + 2)}{\Gamma(\alpha_j)} \int \dots \int \frac{\Gamma (A+2)}{ \prod_{i=1}^{m} \Gamma (\alpha_i^{\prime}) } \prod_{i=1}^{m} \theta_i^{\alpha_i^{\prime} – 1} d\theta_1 \dots d\theta_{m-1} $$

$$ = \frac{\Gamma(A)}{\Gamma(A+1)(A+1)} \frac{\Gamma(\alpha_j + 1)(\alpha_j + 1)}{\Gamma(\alpha_j)} $$

$$ = \frac{\Gamma(A)}{\Gamma(A)A(A+1)} \frac{\Gamma(\alpha_j)\alpha_j(\alpha_j + 1)}{\Gamma(\alpha_j)} $$

$$ = \frac{\alpha_j (\alpha_j + 1)}{A(A+1)} $$

$$ V(\theta_j) = E \left( \theta_j^2 \right) – \left [ E(\theta_j) \right ]^2 \\
= \frac{\alpha_j (\alpha_j + 1)}{A(A+1)} – \left ( \frac{\alpha_j}{A} \right )^2 \\
= \frac{\alpha_j}{A(A+1)} – \frac{\alpha_j^2}{A^2(A+1)}
$$

多変量正規分布

追記予定

多変量正規分布の平均

追記予定

多変量正規分布の分散

追記予定

参考文献

[1] 鈴木・山田 (1996), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃
[2] 原隆 , “微分積分続論 SII-15 クラス” , 九州大学の講義資料
[3] 高校数学の美しい物語 , “ガウス積分の公式の2通りの証明”
[4] 高校数学の美しい物語 , “偶関数と奇関数の意味,性質などまとめ”
[5] 小杉考司 (2015) ,”Cauchy分布について(ベイズ塾例会資料)2015.07.26″, slideshare
[6] MAS3301 Bayesian Statistics

[数理統計学]離散型確率分布の期待値と分散の導出まとめ

はじめに

先日、知人が社労士の試験を受けていたのですが、アプリで問題を色々と解けるものがあるらしく、少しの課金で通勤電車などでの勉強が捗るそうです。統計検定に関してもそのようなものがあればいいなとは思いますが、そんなものを作る気力はないので、よくある分布とその平均と分散に関する導出についてひたすら記していきます。これで少なくとも自分の通勤時の学習が捗ると思われます。なによりも、歳をとっても最低限の数式を扱うスキルがあれば思い出せるレベルで残すことも意識したいですね。

※PC版でないと数式が切れてしまうので、SP版での閲覧はおすすめしません。

今回登場する離散分布

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

  • 1.一様分布
  • 2.ベルヌーイ分布
  • 3.二項分布
  • 4.ポアソン分布
  • 5.超幾何分布
  • 6.幾何分布
  • 7.負の二項分布

1.一様分布

確率変数が有限個の値をとり、それぞれの確率が$$\frac{1}{n}$$である分布。

一様分布の平均

$$E(X)=\frac{1}{n}\sum_{i=1}^{n}X_{i}$$

一様分布の分散

$$V(X)=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\bar X)^2$$

2.ベルヌーイ分布

等確率の原理に従うとして、あるトラックの中にある「ちくわぶ」を取り出していった場合、確率θで「ちくわぶ」は穴が空いておらず、確率(1-θ)で穴が空いているとする。穴の空いていない「ちくわぶ」の割合はベルヌーイ分布に従う。

ベルヌーイ分布の平均

$$E(X)=1 \times \theta + 0 \times (1-\theta ) \ $$
$$ = \theta $$

ベルヌーイ分布の分散

$$V(X)=(1-\theta)^2 \times \theta + (0 – \theta)^2 \times (1-\theta ) \\ $$
$$ = \left[ (1-\theta) \times \theta + \theta^2 \right] \times (1-\theta) $$

$$ = \theta(1-\theta) $$

3.二項分布

一つの事象(「ちくわぶ」に穴が空いているかどうか)に関する、確率θからなるベルヌーイ分布に従う確率変数列を
$$X_1 , X_2, \dots ,X_n$$
として、それらの和(n回試行した際の「ちくわぶ」に穴が空いていない個数)を
$$X = X_1 + X_2 + \dots + X_n$$
とする。その和の従う確率分布は二項分布と呼ばれ、
$$P(X=x):= {}_n C_x \theta^x(1-\theta)^{n-x}$$
と表される。

二項分布の平均

$$E(X)=\sum_{x=0}^{n}x \times {}_n C_x \theta^x(1-\theta)^{n-x} $$

x=0のときは0なので、

$$ = \sum_{x=1}^{n} \frac{n!}{(n-x)!(x-1)!}\theta^x(1-\theta)^{n-x}$$
$$ = n \theta \sum_{x=1}^{n} \frac{(n-1)!}{(n-x)!(x-1)!}\theta^{x-1}(1-\theta)^{n-x}$$
$$ = n \theta \sum_{x=1}^{n} \frac{(n-1)!}{\left[ (n-1) – (x-1) \right]!(x-1)!}$$
$$ \times \theta^{x-1}(1-\theta)^{ (n-1) – (x-1) } $$
ここで
$$y = x-1$$
とおくと、

$$ = n \theta \sum_{y=0}^{n-1} {}_{n-1} C_y \theta^y(1-\theta)^{n-1-y}$$
二項分布の和は1であるという性質を使えば、
$$ = n \theta $$

二項分布の分散

$$ E [ X(X-1) ] = \sum_{x=0}^{n} x(x-1) $$
$$ \times {}_n C_x \theta^x(1-\theta)^{n-x} $$

x = 1, 2のときは0なので、

$$= \sum_{x=2}^{n} \frac{n!}{(n-x)!(x-2)!} \theta^x(1-\theta)^{n-x} $$
$$= \sum_{x=2}^{n} \frac{n!}{\left[ (n-2) – (x-2) \right]!(x-2)!} \theta^x(1-\theta)^{(n-2) – (x-2) } $$
$$= n(n-1)\theta^2 \sum_{x=2}^{n} \frac{n!}{\left[ (n-2) – (x-2) \right]!(x-2)!} \theta^{x-2}(1-\theta)^{(n-2) – (x-2) } $$

ここで
$$y = x -2$$
とおくと、

$$= n(n-1)\theta^2 \sum_{y=0}^{n-2} {}_{n-2} C_y \theta^y(1-\theta)^{(n-2) – y } $$

となり、二項分布の和が1であることから、
$$= n(n-1)\theta^2$$
となる。

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

また、この式より、二項分布に関してはθの絶対値が1よりも小さいことから、平均よりも分散の方が小さくなることがわかる。

4.ポアソン分布

二項分布において、平均値がλとおく。つまり、
$$n\theta = \lambda$$
λが一定でnが非常に大きいケースにおいて、その確率分布はポアソン分布に従う。λが一定でnが非常に大きいケースというのは、θが非常に小さいということになり、その事象が滅多に起きないような対象に対して扱われることが多い。

ポアソン分布の導出は二項分布からネイピア数を用いて導出する方法や、微分方程式を用いた方法などがある。

ポアソン分布の導出(二項分布からver)

$${}_{n} C_x \theta^x(1-\theta)^{n-x}$$
$$ = \frac{n(n-1) \dots (n-x+1)}{x!} \theta^x(1-\theta)^{n-x}$$

ここで、
$$\theta = \frac{\lambda}{n}$$
とすると

$$ = \frac{1\left(1-\frac{1}{n}\right) \dots \left(1-\frac{x-1}{n}\right)}{x!} \left(\frac{\lambda}{n}\right)^x\left(1-\frac{\lambda}{n}\right)^{n-x}n^x$$
$$ = \frac{1\left(1-\frac{1}{n}\right) \dots \left(1-\frac{x-1}{n}\right)}{x!} \lambda^x\left(\left(1-\frac{\lambda}{n}\right)^{-\frac{n}{\lambda}}\right)^{-\lambda} \left(1-\frac{\lambda}{n}\right)^{-x}$$

nを無限大にすると、分母にのみnがある項は0になり、また、ネイピア数の定義から
$$ \lim_{n\to\infty}{}_{n} C_x \theta^x(1-\theta)^{n-x} = \frac{\lambda^x e^{-\lambda}}{x!}$$
となる。
$$P(X=x):= \frac{\lambda^x e^{-\lambda}}{x!}$$
この確率分布がポアソン分布と呼ばれる。(x=0,1,2,…)

実際に、確率分布の性質を確かめると、

$$ \sum_{x=0}^{\infty} P(X=x) = \sum_{x=0}^{\infty} \frac{\lambda^x e^{-\lambda}}{x!} $$
$$ = e^{-\lambda} \sum_{x=0}^{\infty} \frac{\lambda^x}{x!} $$
指数の和の公式より、
$$ = e^{-\lambda} e^{\lambda} = 1 $$
よって、確率分布であることがわかる。

ポアソン分布の導出(微分方程式ver)

abrahamcowさんのブログにまさにそれが書かれていました。図も載っていて詳しくて良いです。
他にも、英語の文献だとA Quick Way to See that the Poisson Distribution is the Appropriate Mathematical Formulation for a Counting Process with Constant Rate and Intensityなどがあります。

これらの文献を見ていただければ導出について理解できると思います。

以下、これらの文献を参考にして式展開したものを記しておきます。

まず、3つの仮定をおきます。

  • 1.互いに背反する期間では現象が起こる確率は独立。
    つまり、X(t)とX(t+h)-X(t)は独立。
  • 2.非常に小さい区間で現象が起こる確率はその区間の長さに比例する。
    $$P\{X(t+h) – X(t) = 1 \} = \lambda h + o(h)$$
    ここでo(h)は
    $$\lim_{n\to 0} \frac{o(h)}{h}=0$$
    となる関数とする。また、λ > 0とする。
  • 3.非常に小さい区間で現象が2回以上起こる確率はその区間で現象が1回起こる確率に比べて、無視できるほど小さい。
    $$P\{X(t+h) – X(t) = k\} = o(h) $$
    ただし、k > 1とする。

ここで、X(t)がrになる確率を以下の表現で表す。
$$P_r (t) = P\{ X(t) = r \} $$

事象
$$X(t+h)=r$$
は以下の3つに場合分けできる。

  • 1.$$ \{ X(t) = r,X(t+h) – X(t) = 0 \}$$
    この
    $$X(t+h) – X(t) = 0$$
    は仮定2と仮定3の式を1から差っ引いたものによって確率を計算できる。また、仮定1から互いに背反する期間では現象が起こる確率は独立のため、同時確率は各々の積で表される。

$$P\{ X(t) = r \} (1 – \lambda h – o(h) – o(h) )$$

  • 2.$$ \{ X(t) = r-1,X(t+h) – X(t) = 1 \}$$
    仮定2より、
    $$P\{ X(t) = r-1 \} (\lambda h + o(h))$$
    によって表される。

  • 3.

$$ \\{ X(t+h) = r-k,X(t+h) – X(t) = k \\} , k > 1 $$

仮定3より、
$$P\{ X(t+h) = r-k \} o(h)$$
によって表される。

これらは期間が違うため、仮定より互いに背反であることから、3つに場合分けの確率の和で表現する。
$$P_r(t+h) = P\{ X(t) = r \} $$
$$= P\{ X(t) = r \} (1 – \lambda h – o(h) – o(h) ) $$

$$+ P\{ X(t) = r-1 \} (\lambda h + o(h)) + P\{ X(t+h) = r-k \} o(h)$$

$$ = P_r(t)(1-\lambda h) + P_{r-1}(t)\lambda h + o(h) $$

$$P_r(t+h) – P_r(t) = -\lambda h P_r(t) + P_{r-1}(t)\lambda h + o(h) $$
$$\frac{P_r(t+h) – P_r(t)}{h} = -\lambda P_r(t) + \lambda P_{r-1}(t) + \frac{o(h)}{h} $$

hの極限をとると、

$$\lim_{h\to 0} \frac{P_r(t+h) – P_r(t)}{h} = -\lambda P_r(t) + \lambda P_{r-1}(t) $$

となる。

これらは
r = 0のとき
$$\frac{P_0(t)}{dt} = – \lambda P_0(t)$$

r ≠ 0のとき

$$\frac{P_r(t)}{dt} = -\lambda P_r(t) + \lambda P_{t-1}(t), r = 1,2,\dots$$

となる微分方程式として表される。

まず、r=0の場合の微分方程式を解くと、

$$\frac{P_0(t)}{dt} = – \lambda P_0(t) \Rightarrow \frac{d P_0(t)}{P_0(t)} = -\lambda dt $$
両辺積分すると、
$$ \int \frac{d P_0(t)}{P_0(t)} = -\lambda \int dt $$

$$ \log P_0(t) = -\lambda t + C$$
(Cは積分定数)

$$P_0(t) = e^{-\lambda t +C}$$
Cは任意なので、
$$e^C=C$$
として、
$$P_0(t) = Ce^{-\lambda t}$$
となる。
t = 0で1件も発生しない確率は1なので、
$$P_0(0)=C=1$$
から、C = 1であることがわかり、
$$P_0(t) = e^{-\lambda t}$$
となる。

ここで、定数係数の1階線形微分方程式の解の公式は、
$$y^{\prime} + p(x)y = q(x)$$
に対して、
$$y = e^{-p(x)x}\left[ \int q(x) e^{p(x)x}dx + C \right]$$
で表される。
この解の公式を今回のr=1の場合の微分方程式に当てはめると、
$$y^{\prime} = P_1(t) , p(x) = -\lambda , q(x) = \lambda P_0(t) $$
であるから、
$$P_1(t) = e^{-\lambda t}\left[ \int \lambda P_0(t) e^{\lambda t} dt + C \right] $$
となる。ここに、
$$P_0(t)$$
を代入すると、
$$ = e^{-\lambda t}\left[ \int \lambda e^{-\lambda t} e^{\lambda t} dt \right] $$
$$ = e^{-\lambda t}\left[ \int \lambda dt + C \right] $$
$$ = e^{-\lambda t} \lambda t + e^{-\lambda t}C $$
となる。
t=0においてX(0)=1の確率は0なので、C=0となり、
$$P_1(t) = e^{-\lambda t} \lambda t $$

同様にしてr=2の場合、
$$P_2(t) = \frac{e^{-\lambda t \left( \lambda t \right)^2}}{2} $$
r=3の場合、
$$P_3(t) = \frac{e^{-\lambda t \left( \lambda t \right)^3}}{2 \times 3} $$
となる。

ここで、任意の整数r=kにおいて、
$$P_k(t) = \frac{e^{\lambda t} \left( \lambda t \right)^k }{k!}$$
が成り立つとする。

r=k+1のとき、
$$P_{k+1}(t) = \frac{e^{\lambda t} \left( \lambda t \right)^{k+1} }{(k+1)!}$$
である。
これはポアソン分布であり、定数係数の1階線形微分方程式よりポアソン分布が導出できることが示された。

ポアソン分布の平均

$$E(X)=\sum_{x=0}^{\infty}x \frac{e^{-\lambda}\lambda^x}{x!}$$
$$ =\sum_{x=0}^{\infty} \frac{e^{-\lambda}\lambda^x}{(x-1)!}$$
$$ =\lambda \sum_{x=1}^{\infty} \frac{e^{-\lambda}\lambda^{x-1}}{(x-1)!} = \lambda $$

ポアソン分布の分散

$$E(X(X-1)) = \sum_{x=0}^{\infty}x(x-1) \frac{e^{-\lambda}\lambda^x}{x!}$$
$$ = \sum_{x=0}^{\infty} \frac{e^{-\lambda}\lambda^x}{(x-2)!}$$
$$ = \lambda^2 \sum_{x=2}^{\infty} \frac{e^{-\lambda}\lambda^{x-2}}{(x-2)!} =\lambda^2 $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ =E[X(X-1)] + E(X) -[ E(X) ]^2$$
$$ = \lambda^2 + \lambda – \lambda^2 = \lambda $$

ポアソン分布の再生性

XとYは独立にそれぞれポアソン分布λ1とλ2に従うとする。
$$P(Z=z) = \sum_{x=0}^{z} P(X=x, Y=z-x) $$
$$ = \sum_{x=0}^{z} P(X=x)P(Y=z-x) $$

$$ = \sum_{x=0}^{z} \frac{e^{-\lambda_1} \lambda_1^x }{x!} \frac{e^{-\lambda_2} \lambda_2^{z-x} }{(z-x)!} $$
$$ = \frac{e^{-(\lambda_1 + \lambda_2 )}(\lambda_1 + \lambda_2)^z}{z!} \sum_{x=0}^{z}\frac{z!}{x!(z-x)!}\left( \frac{\lambda_1}{\lambda_1+\lambda_2} \right)^x \left( \frac{\lambda_2}{\lambda_1+\lambda_2} \right)^{z-x} $$

二項分布の和は1になるので以下のように表される。
$$ = \frac{e^{-(\lambda_1 + \lambda_2 )}(\lambda_1 + \lambda_2)^z}{z!} $$
このように表されることを再生性があると呼ぶ。

5.超幾何分布

N個の製品からなる仕切りのなかで、M個の不良品があるとする。
非復元抽出でn個を抽出した際の、不良品の数をXとした際の確率分布は以下のように表される。

$$P(X=x)=\frac{{}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

xは非負の整数で

$$max \{ 0, n -(N-M) \} \leq x \leq min \{ n, M \}$$

に従うものとする。

ここで製品の数Nが標本nに比べて十分に大きい場合を考える。
$$\frac{{}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

$$ = \frac{M!}{x!\left( M-x \right)!} \times \frac{\left( N – M \right)! }{\left( n-x\right)! \left[ N-M -\left( n-x\right) \right]!} \times \frac{n! \left(N-n\right)!}{N!} $$
$$ = \frac{\left[ M(M-1)\dots (M-x+1)\right] }{x!(n-x)! \times \left[ N(N-1)\dots (N-n+1) \right]}$$
$$ \times \left[ (N-M)(N-M-1)\dots (N-M-n+x+1)\right]\times n! $$

組み合わせを使ってまとめると、分母と分子の数はn個であるから、分母分子にn個だけ1/Nを掛けると以下のようになる。
$$ = \frac{{}_n \mathrm{C} _x \times \frac{M}{N} \left( \frac{M}{N} – \frac{1}{N} \right) \dots \left( \frac{M}{N} – \frac{x-1}{N} \right) }{1 \left( 1 – \frac{1}{N} \right) \dots \left( 1 – \frac{n-1}{N} \right) } $$

$$ \times \left( \frac{N-M}{N} \right) \times \left( \frac{N-M}{N} – \frac{1}{N} \right) \dots \left( \frac{N-M}{N} – \frac{n-x-1}{N} \right) $$

ここで
$$ \theta = \frac{M}{N}$$
から、
$$ M = N\theta$$
として代入すると、

$$ = \frac{{}_n \mathrm{C} _x \times \theta \left( \theta – \frac{1}{N} \right) \dots \left( \theta – \frac{x-1}{N} \right) }{1 \left( 1 – \frac{1}{N} \right) \dots \left( 1 – \frac{n-1}{N} \right) } $$
$$ \times \left( 1-\theta \right) \times \left( (1-\theta) – \frac{1}{N} \right) \dots \left( (1-\theta) – \frac{n-x-1}{N} \right) $$

となる。
ここで、Nの極限を取ると以下のように、分母が全て1になり、以下のように表される。

$$ \lim_{N\to\infty}\frac{{}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$
$$ = {}_n \mathrm{C} _x \theta^x \left( 1 – \theta \right)^{n-x} $$

これは二項分布と同じ形となる。
つまり、超幾何分布において、製品の数が十分に大きいと二項分布と同じになる。

超幾何分布の平均

$$ E(X) = \sum_{x=0}^{n} x \frac{ {}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

$$ = \frac{1}{ {}_N \mathrm{C} _n } \sum_{x=0}^{n} x \frac{M!}{(M-x)!x!} $$
$$ \times \frac{(N-M)!}{ \left[ N-M-(n-x) \right]! (n-x)!} $$

$$ = \frac{M}{ {}_N \mathrm{C} _n} \sum_{x=0}^{n} \frac{(M-1)!}{\left[ (M-1)-(x-1) \right]!(x-1)!} $$

$$ \times \frac{\left[ (N-1)-(M-1) \right]!}{\left[ (N-1)-(M-1)-\left[(n-1)-(x-1)\right] \right]! \left[ (n-1)-(x-1) \right]!}$$

$$ = \frac{M}{ {}_N \mathrm{C} _n} $$
$$ \times \sum_{x=0}^{n} {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)} $$

$$ = \frac{M}{ \frac{N!}{(N-n)!n!} } $$
$$ \times \sum_{x=0}^{n} {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)} $$

$$ = \frac{nM}{ \frac{(N-1)!N}{\left[ (N-1)-(n-1)\right]!(n-1)!} } $$
$$ \times \sum_{x=0}^{n} {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)}$$

$$ = \frac{nM}{N} $$
$$ \sum_{x=0}^{n} \frac{ {}_{M-1} \mathrm{C} _{x-1} \times {}_{(N-1) – (M-1)} \mathrm{C} _{(n-1)-(x-1)} }{ {}_{N-1} \mathrm{C} _{n-1} } $$

超幾何分布の和は確率分布のため、1になることから、
$$ = \frac{nM}{N}$$
となる。これは二項分布の平均と同じとなる。

いちいち計算するのが面倒なので、二項係数の公式をここに載せておきます。

$$ {}_n \mathrm{C} _r = {}_n \mathrm{C} _{n-r} $$
$$r {}_n \mathrm{C} _r = n {}_{n-1} \mathrm{C} _{r-1}$$
$$ {}_n \mathrm{C} _r = {}_{n-1} \mathrm{C} _{r} + {}_{n-1} \mathrm{C} _{r-1}$$

超幾何分布の分散

$$E(X(X-1))= \sum_{x=0}^{n} x(x-1) $$
$$ \times \frac{ {}_M \mathrm{C} _x \times {}_{N-M} \mathrm{C} _{n-x}}{{}_N \mathrm{C} _n} $$

二項係数の公式より、

$$ = \sum_{x=0}^{n} (x-1)$$
$$ \times \frac{ M \times {}_{M-1} \mathrm{C} _{x-1} \times {}_{N-M} \mathrm{C} _{n-x} \times n }{ N \times {}_{N-1} \mathrm{C} _{n-1} } $$

$$ = \sum_{x=0}^{n} M(M-1) \times n(n-1) $$
$$ \times \frac{ {}_{M-2} \mathrm{C} _{x-2} \times {}_{(N-2)-(M-2)} \mathrm{C} _{(n-2)-(x-2)} }{ N(N-1) \times {}_{N-2} \mathrm{C} _{n-2}} $$

超幾何分布の和が1であることから、
$$ n(n-1) \frac{M(M-1)}{N(N-1)} $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ =E[X(X-1)] + E(X) -[ E(X) ]^2$$
$$ = n(n-1) \frac{M(M-1)}{N(N-1)}+ \frac{nM}{N} – \frac{n^2M^2}{N^2} $$
$$ = \frac{nM}{N} \left[ (n-1) \frac{M-1}{N-1} + 1 – \frac{nM}{N} \right] $$
$$ = \frac{nM}{N} \left( 1 – \frac{M}{N} \right) \left( \frac{N-n}{N-1} \right) $$

超幾何分布の平均は二項分布と同じであったが、分散に関しては、
$$ \left( \frac{N-n}{N-1} \right) $$
だけ異なる。

6.幾何分布

復元抽出で、不良品が出るまで検査した際の観測された良品の数Xの確率分布は以下の幾何分布に従う。

$$P(X=x)=\theta ( 1-\theta)^x $$

幾何分布は無記憶性という性質を持っており、それは幾何分布のみが持つとされている。無記憶性はある時点から先に時点を進めた際に、過去の時点の影響が残らないことを指している。

$$P(X=x+h | X \geq x) = \frac{P(X=x+h)}{ \sum_{y=x}^{\infty} P(X=y) }$$

$$ = \frac{ \theta ( 1-\theta)^{x+h} }{ \sum_{y=x}^{\infty} \theta ( 1-\theta)^y } $$
$$ = \frac{ \theta ( 1-\theta)^{h} }{ \sum_{y=x+1}^{\infty} \theta ( 1-\theta)^y } $$
超幾何分布の和は1であるから、
$$ = \theta ( 1-\theta)^h $$
$$ = P(X=h) $$

幾何分布の平均

$$ E(X) = \sum_{x=0}^{\infty} x \theta ( 1-\theta)^x $$
$$ = \sum_{x=1}^{\infty} x \theta ( 1-\theta)^x $$
$$ = \theta \sum_{x=1}^{\infty} x ( 1-\theta)^x $$
$$ = \theta \frac{d \sum_{x=1}^{\infty} (1-\theta)^x}{d(1-\theta)}\times (1-\theta) $$
無限等比級数の和の公式より、
$$ = \theta \frac{d \left( \frac{1-\theta}{1-(1-\theta)} \right) }{d(1-\theta)}\times (1-\theta) $$
$$ = \theta \frac{d \frac{1-\theta}{\theta}}{d \theta} \times \frac{d \theta}{d(1-\theta)}\times (1-\theta)$$
$$ = \theta \frac{-\theta – (1 – \theta)}{\theta^2} \times \frac{1}{-1}\times (1-\theta)$$
$$ = \frac{1-\theta}{\theta}$$

幾何分布に関しては、べき級数の和を使うことが多いので、毎回書くのも大変なためここで紹介しておく。

$$ 1 + r + r^2 + r^3 + \dots = \sum_{k=0}^{\infty} r^k = \frac{1}{1-r}$$

両辺をrで微分すると、

$$ 1 + 2r + 3\times r^2 + \dots = \sum_{k=1}^{\infty} kr^{k-1} = \frac{1}{(1-r)^2}$$

さらに両辺をrで微分すると、

$$ 2 + 3\times2 r + 4\times 3 r^2 + \dots = \sum_{k=2}^{\infty} (k-1)kr^{k-2} = \frac{2}{(1-r)^3}$$

となる。

幾何分布の分散

$$ E(X(X-1)) = \sum_{x=0}^{\infty} x(x-1) \theta ( 1-\theta)^x $$
$$ = \theta (1-\theta)^2 \sum_{x=2}^{\infty} x(x-1) ( 1-\theta)^x $$
べき級数の和を用いると、
$$ = \theta (1-\theta)^2 \frac{2}{\left[ 1 – (1-\theta) \right]^3 }$$
$$ = \frac{2(1-\theta)^2}{\theta^2} $$

$$V(X)=E(X^2)-[ E(X) ]^2$$
$$ = E[X(X-1)] + E(X) -[ E(X) ]^2$$

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

$$ = \frac{1-\theta}{\theta^2}$$

θの絶対値は1以下なので、幾何分布においては、分散が平均よりも大きくなることがわかる。

7.負の二項分布

r個の不良品が見つかるまでに観測した良品の数Xが従う確率分布で、以下のように表される。

$$P(X=x)={}_{r+x-1} \mathrm{C} _x \theta^r (1-\theta)^x$$
$$ x=0,1,2,\dots $$
負と呼ばれる所以は、以下のように組み合わせを負で表現できるところにある。

以下では実際に負で表現できるかを示す。
$$P(X=x)={}_{-r} \mathrm{C} _x \theta^r(\theta – 1)^x$$

$$ {}_{r+x-1} \mathrm{C} _x $$
$$ = \frac{(r+x-1)(r+x-2)\dots(r+1)r}{x!} $$

分子の数はx個あるので、おのおの-1を掛けると

$$ = (-1)^x \frac{(-r)(-r-1)\dots(-r-(x-2))(-r-(x-1))}{x!}$$

$$ = (-1)^x \frac{(-r)!}{(-r-x)!x!}$$
$$ = (-1)^x {}_{-r} \mathrm{C} _x $$
よって、
$$P(X=x)=(-1)^x {}_{-r} \mathrm{C} _x \theta^r (1-\theta)^x $$
$$ = {}_{-r} \mathrm{C} _x \theta^r (\theta-1)^x $$

負の二項分布の平均

$$E(X) = \sum_{x=0}^{\infty} x \times {}_{r+x-1} \mathrm{C} _x \theta^r (1-\theta)^x $$

$$ =\sum_{x=1}^{\infty} \frac{(r+x-1)!\theta^r (1-\theta)^x}{(x-1)!(r-1)!} $$
ここでx = y + 1 として
$$ =\sum_{y=0}^{\infty} \frac{(r+y)!\theta^r (1-\theta)^{y+1}}{(y+1-1)!(r-1)!} $$
$$ = \frac{r(1-\theta)}{\theta} \sum_{y=0}^{\infty} \frac{(r+y)! \theta^{r+1} (1-\theta)^{y}}{y!r!} $$

$$ = \frac{r(1-\theta)}{\theta} \times $$
$$ \sum_{y=0}^{\infty} \frac{(r+y)(r+y-1)\dots (r+1)}{y!} $$
$$ \times \theta^{r+1} (1-\theta)^{y} $$

$$ = \frac{r(1-\theta)}{\theta} \sum_{y=0}^{\infty} {}_{r+y} \mathrm{C} _{y} \theta^{r+1} (1-\theta)^y $$

ここで負の二項分布の和は1であることから、
$$ = \frac{r(1-\theta)}{\theta} $$
となる。

負の二項分布の分散

$$ E(X(X-1)) = \sum_{x=0}^{\infty} x(x-1) $$
$$ \times {}_{r+x-1} \mathrm{C} _x \theta^r (1-\theta)^x $$

$$ = \sum_{x=2}^{\infty} \frac{(x+r-1)(x+r-2)\dots r}{(x-2)!} \theta^r (1-\theta)^x $$

ここでx = y + 2 として

$$ = \sum_{y=0}^{\infty} \frac{(y+r+1)(y+r)\dots r}{y!} \theta^r (1-\theta)^{y+2} $$
$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} \sum_{y=0}^{\infty} \frac{(y+r+1)(y+r)\dots (r+2)}{y!} \theta^{r+2} (1-\theta)^{y} $$

$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} $$
$$ \times \sum_{y=0}^{\infty} {}_{r+y+1} \mathrm{C} _{y} \theta^{r+2} (1-\theta)^y $$

ここで負の二項分布の和は1であることから、
$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} $$

$$V(X)=E(X^2)-[ E(X) ]^2 $$
$$ = E[X(X-1)] + E(X) -[ E(X) ]^2$$

$$ = \frac{r(r+1)(1-\theta)^2}{\theta^2} + \frac{r(1-\theta)}{\theta} – \frac{r^2(1-\theta)^2}{\theta^2} $$

$$ = \frac{(1-\theta)^2r + (1-\theta)\theta r}{\theta^2} $$
$$ = \frac{(1-\theta)\left[ (1-\theta)r + \theta r \right]}{\theta^2} $$
$$ = \frac{r(1-\theta)}{\theta^2} $$

θの絶対値は1以下なので、負の二項分布においては、分散が平均よりも大きくなることがわかる。

おわりに

以上、これだけ導出を丁寧に書けば後で思い出せるだろうなというレベルで残したつもりですが、冗長的なところもあると思います。おかしなところがあったらご指摘ください。次は連続分布を扱う予定ですが、結構LaTeX書くの疲れますね。

参考文献

[1] 鈴木・山田 (1996), 『数理統計学―基礎から学ぶデータ解析』, 内田老鶴圃
[2] 緑川章一, “負の二項分布”
[3] abrahamcow (2014),”微分方程式によるポアソン分布の導出”
[4] 物理のかぎしっぽ, “定数係数1階線形微分方程式”
[5] 高校数学の美しい物語, “超幾何分布の意味と期待値の計算”