検索
連載

【Pythonで学ぶデータ分析】母平均に差があるかどうかをベイズt検定で調べる 〜 運動部と非運動部の体力差はあるのか？：やさしい推測統計（ベイズ統計編）

運動部と非運動部の生徒の体力テストを例に、母平均に差があるかどうかをベイズ統計により検定します。古典的なt検定のp値に代わるものとしてベイズ因子を利用します。『社会人1年生から学ぶやさしいデータ分析』ベイズ統計編の第6回です。

Share
Tweet
LINE
Hatena
「やさしい推測統計（ベイズ統計編）」のインデックス

連載目次

　母平均の差の検定をベイズ統計で行うには？　連載『社会人1年生から学ぶ、やさしいデータ分析』のベイズ統計編、第6回となる今回は、古典的なt検定に代わる方法として、ベイズt検定を行う方法を解説します。第3回では、1群の母平均が基準値と等しいかどうかをベイズ検定により確かめましたが、今回は独立した2群の母平均に差があるかどうかを確かめようというわけです。最後に、対立仮説や帰無仮説が正しい確率を求めます。これは古典的な検定ではどうしてもできなかったことです。ベイズ統計の圧倒的なメリットが実感を持って理解できると思います。

　この連載では、簡単な事例を通してベイズ統計の考え方と分析の進め方を解説します。新しい用語や考え方が幾つも出てきますが、全てを理解しなくても大丈夫です。登場するたびに、分かりやすく丁寧に説明するので、その時点で分からないことがあっても気にせず、先に進んでください。回を追うごとに、少しずつ理解が深まります。どんな考え方で、どんな手順で分析を進めるのか、といった大きな流れを捉えてください。

連載：

『社会人1年生から学ぶ、やさしい推測統計（ベイズ統計編）』

社会人1年生から学ぶ、やさしい推測統計（ベイズ統計編）

　この連載では、データをさまざまな角度から分析し、その背後にある有益な情報を取り出す方法を学ぶ『社会人1年生から学ぶ、やさしいデータ分析』シリーズの「記述統計と回帰分析編」「確率分布編」「推測統計（区間推定編・仮説検定編）」に続く「ベイズ統計編」です。

　これまでの推測統計を土台に、近年活用が広がっているベイズ統計の考え方と分析手順を、古典的な手法との違いを整理しながら解説します。初めての方でも無理なく理解できるよう具体例を通して進めるとともに、ベイズ的なアプローチの特徴やメリットを実感できる構成とし、「どのように考え、どう使い分けるのか」に重点を置いて解説していきます。

羽山博
羽山博

筆者紹介： IT系ライターの傍ら、かつて、非常勤講師として東大で情報・プログラミング関連の授業を、一橋大でAI関連の授業を担当。まだ発展途上のようだが、生成AIへの質問と回答を分野ごとにまとめ、体系的に整理、補足しつつ、さらに個別の学びを促す記述を加えた資料を一貫して作成できるツールが実用レベルで本格的に使えるようになると、書籍の概念が変わるのではないかと思ったりしている。一斉授業だったものが、ユーザーに合った個別指導のできる書籍といったイメージ。そうなると、出版社のビジネスモデルも大きく変わり、私の仕事もなくなってしまいそうだけど、まあ、そのときは、縁側でお茶でもすすりながら、庭の景色を眺めて暮らすことにしようと思う。ウチには縁側も庭もないけど。


母平均の差をベイズt検定によって分析する

　第1回〜第5回までは、1群のサンプルについて母数をベイズ推定したり、基準値との差についてベイズ検定を行ってきたりしました。今回は2群の母平均について、差があるかどうかをベイズt検定により確かめます。今回は、両側検定と片側検定、標準偏差が等しい場合／等しいとは考えられない場合の検定、さらには仮説が正しい確率を求める……というように、話がかなり広がりますが、一歩ずつ進めるために、標準偏差が等しい場合の両側検定からスタートします。それが分かれば、少しばかり補正や計算を加えるだけで他の分析もできます。

　最初に、図1で今回の事例と今回やることを確認しておきましょう。図には効果量δを推定したり、それを基にベイズ因子を求めるといったことが書かれています。この例での「平均に差がある」という結果は、正確には「効果量に差がある」ということですが、詳細についてはこの後少しずつお話ししていきます。

今回やること
図1　独立した2群の母平均の差をベイズt検定により確かめる
運動部に属している生徒と運動部に属していない生徒の体力テストのデータ（14歳男子の架空データ）。成績は握力、反復横跳び、持久走、50m走、ハンドボール投げなどの得点（1〜10点）を合計したもの。このデータを基に、成績の母平均に差があるかどうかをベイズt検定によって確認する。最後に、対立仮説が正しい確率を求める方法も解説する（古典的なt検定では仮説が正しい確率を求めることはできなかったが、ベイズ統計ではそれができる！）。

　今回のサンプルデータは、スポーツ庁による体力・運動能力調査の結果（運動部・スポーツクラブ所属の有無と体格測定・テストの結果・14、15歳男・女子）を参考にして作成した架空のデータです。実際の値は整数ですが、この架空データでは、小数点第一位までの値としました。

　それぞれの群は独立で、正規分布に従っており、標準偏差は等しいものとします。「独立」とは、簡単に言えば関係がないということです。具体的には、運動部の群と非運動部の群は別々の人だということです。例えば、先頭の値の44.8点と45.8点は別の人の成績です。つまり、ここでの前提は古典的なt検定での「ステューデントのt検定」と同じです。ちなみに、同じ人が2回体力テストを受けて、1回目が44.8点で2回目が45.8点だとすると、独立ではなく「対応のある」データとなります。


AI博士

　実際のところ、2群の標準偏差が等しいと考えることにはやや無理がありますが、ここでは手順を学ぶことを主旨とし、このまま進めていきます。なお、標準偏差が等しいと考えられない場合（「ウェルチのt検定」に対応する場合）は、平均の標準偏差を使います。これについては、後ほど補足します。また、対応のあるデータについて、平均の差を検定する例は、次回解説します。


　さて、今回は効果量δに注目して分析を行います。その理由を簡単に説明しましょう。ここでの効果量とは、サンプルから計算されたいわゆるCohen's dと呼ばれる値ではなく、母集団の効果量です。つまり、以下の式で表されます。

　分析の方法としては、平均そのものの差を考えても構わないのですが、効果量を使うと次のようなメリットがあります。

  • 母数がδだけで済む
    • 帰無仮説H0はμ1 ＝ μ2でなく、δ ＝ 0と表せる
    • 対立仮説H1はμ1 ≠ μ2でなく、δ ≠ 0と表せる
  • 単位が異なっても統一的に分析ができる
    • 例えば、δ ＝ 1ということは、元のデータの単位がいくらであっても、μ1とμ2の差が標準偏差1つ分であるということを表す

　つまり、帰無仮説がδ ＝ 0のように、ある1つの点（実際にはその点にピークのあるスパイク）で表せるというわけです。また、対立仮説についても、母数が1つだけなら、事前分布が比較的シンプルに表せます。対立仮説のδ ≠ 0というのは、あくまで古典的な検定に合わせた表現で、正確には、δが何らかの分布に従うということでしたね。

　効果量の事前分布としては、コーシー分布が一般的に使われます。コーシー分布のパラメーターは中心を表すαと、スケール（幅）を表すβの2つです。図2で確率密度関数の形を見ておき、その後で具体的な手順に移りましょう。

コーシー分布
図2　コーシー分布の確率密度関数
線が重なってやや見づらいが、中心αとスケールβの値を変えてコーシー分布の確率密度関数をプロットした例。スケールが大きくなると、グラフが平坦になる（弱情報な事前分布となる）。比較のために標準正規分布（standard normal）の確率密度関数も併せて表示してある。コーシー分布は標準正規分布よりも裾が重い（中心から離れた位置にもある程度の高さがある）ことが分かる。効果量の事前分布としてはα＝0、β＝0.707のコーシー分布が一般的に使われる

　図2のグラフを描くプログラムは、サンプルファイルに含めてあります。リンクを開くと、説明やPythonのプログラムが表示されるので、最初のコードセルをクリックし、［Shift］＋［Enter］キーを押してコードを実行すればグラフが表示されます。

事前分布と尤度関数を定義して、サンプリングを行う

　では、少しずつベイズ推定と検定のための手順を見ていきます。データについては、そのまま使っても構わないのですが、あらかじめ平均を0、標準偏差を1にしておく（標準化しておく）と、事前分布の定義が簡潔になります。コードの全体像は後でまとめて見ることとして、標準化を行う部分を抜き出しておきます。これは簡単ですね。サンプルファイルの2番目のコードセルをクリックし、［Shift］＋［Enter］キーを押して、実行しておいてください。

import numpy as np

x1 = [44.8, 44.5, 36.4, 50.4, 40.8, 48.5, 53.4, 53.2, 47.1, 70.5,
      63.7, 67.3, 68.2, 56.4, 55.0, 54.6, 52.7, 67.0, 57.1, 35.9,
      53.4, 49.7, 44.1, 50.7, 56.4, 44.0, 60.6, 52.2, 55.8, 59.8]
x2 = [45.8, 39.1, 43.6, 61.6, 46.1, 47.8, 49.6, 47.6, 40.8, 44.3,
      30.0, 43.0, 44.3, 49.3, 35.9, 58.5, 42.7, 37.2, 39.7, 55.8,
      24.1, 53.6, 43.8, 33.7, 25.9, 38.9, 40.1, 52.9, 44.3, 41.1]

# x1とx2を標準化する
x_all = np.concatenate((x1, x2))  # x1とx2をつないで1つにまとめる
x_mean = x_all.mean()  # 全体の平均
x_std = x_all.std()  # 全体の標準偏差
x1_standardized = (x1 - x_mean) / x_std  # 標準化されたx1
x2_standardized = (x2 - x_mean) / x_std  # 標準化されたx2

リスト1　データを標準化する
標準化するには、各データから全体の平均を引き、全体の標準偏差で割るとよい。これで全体の平均が0、全体の標準偏差が1にそろえられる。以降の分析は標準化されたデータを使って行うことになる。

　続いて、事前分布と尤度関数の定義です。おさらいからスタートしましょう。毎回掲載しているベイズの定理は以下の通りでした。数式の意味が理解できなくても問題ありません。式の形だけが分かれば十分です。

　周辺尤度は定数でしたね。今のところ、周辺尤度を計算する必要はないので、比例するという意味の∝という記号を使って、以降、必要に応じて以下のように表すことにします。

　要するに、事前分布と尤度関数が定義できれば、事後分布が求められるというわけです。というわけで、事前分布と尤度関数を定義したいのですが、今回は、少し段階を踏む必要があるので、最初に全体像を図3で見ておきます。流れだけをざっと見ておいてください（なんとなく眺める程度で十分です）。

事前分布と尤度関数
図3　事前分布と尤度関数の考え方
注目しているパラメーターは効果量δ。その事前分布についてはコーシー分布とする。また、全体の平均と標準偏差の分布についても定義しておく。これらは、尤度関数のパラメーターμ1とμ2を求めるために使う補助的なパラメーターとなっている。μ1とμ2は、全体の平均μと効果量δ、共通の標準偏差σから求めることができる。計算方法については後述する。

　（2）式や図3では、事前分布や尤度関数のパラメーターをθで代表させて表してありますが、複数のパラメーターの場合もあります。例えば、上の事前分布π(θ)のθを個別に表すと、π(δ, μ, σ)となります。ちなみに、それぞれのパラメーターが独立であれば（この例では独立です）、積で表せるので、π(δ, μ, σ)＝π(δ)π(μ)π(σ)となります。

事前分布を定義する

　では、図3の流れに沿って、事前分布を定義しましょう。効果量δの事前分布としては、一般的に使われる中心が0、スケールが0.707のコーシー分布とします。全体の平均μについては、データが標準化されているので、平均0、標準偏差1の正規分布とすればいいでしょう。共通の標準偏差σについては、元の正規分布の標準偏差を1とした半正規分布とします。μとσは上でも少し触れたように、尤度関数のパラメーターμ1とμ2を求めるための補助的なパラメーターです。μ1とμ2については尤度関数の定義と併せて後述することとして、ここまでをまとめると、以下のようになります。

  • 効果量δ： α＝0、β＝1/√2＝0.707のコーシー分布
  • 共通の標準偏差σ： 元の正規分布の標準偏差がσ＝1である半正規分布
  • 全体の平均μ： 平均μ＝0、標準偏差σ＝1の正規分布

　これらの定義を書くと以下のようになります（リスト2）。

import pymc as pm

with pm.Model() as model:
    # 事前分布
    delta = pm.Cauchy("delta", alpha=0, beta=0.707)  # 効果量
    sigma = pm.HalfNormal("sigma", sigma=1)  # 共通の標準偏差（半正規分布とする）
    mu = pm.Normal("mu", mu=0, sigma=1)  # 全体の平均（標準化されているので、平均は0、標準偏差は1としておく）

リスト2　事前分布の定義
ここまで見てきた事前分布の定義を素直に書いたもの。コーシー分布はpymcモジュールのCauchyクラスを利用して定義する。引数alphaに分布の中心を、引数betaにスケールを指定すればよい。スケールとしては0.707という値が使われることが多い。標準偏差は正の値しか取らないので半正規分布（HalfNormalクラス）とした。全体の平均については、データが標準化されているので、平均0、標準偏差1の標準正規分布を想定する（Normalクラス）。

尤度関数を定義する

　続いて、尤度関数の定義です。尤度関数としては、母集団が正規分布しているという前提なので、正規分布とします。パラメーターは、それぞれの群の平均μ1とμ2です。これらについては、効果量δと全体の平均μ、全体の標準偏差σから計算して求めます。図4に示したような方法です。

尤度関数のパラメーターを計算によって求める
図4　尤度関数のパラメーターを計算によって求める
図の下に示したように、μ1−μ2の幅が、δ × σとなるので、μ1は全体の平均μから、＋(δ × σ)/ 2だけ離れていることが分かる（図の右側の式）。同様に、μ2は、全体の平均μから、−(δ × σ)/ 2だけ離れている（図の左側の式）。

　このように、計算によって新たな変数を定義する場合、pymcモジュールのDeterministic関数を使って記述します（リスト3）。

    mu1 = pm.Deterministic("mu1", mu + delta * sigma / 2)
    mu2 = pm.Deterministic("mu2", mu - delta * sigma / 2)

リスト3　尤度関数で使われるパラメーターを計算で求める
Deterministic関数は、実際にはpymc.model.core.Deterministic関数だが、pymc.Deterministicと記述すればよい。mu1は全体の平均（mu）に効果量（delta）×共通の標準偏差（sigma）の半分を加えたものとなる。mu2は、逆に、全体の平均（mu）から効果量（delta）×共通の標準偏差（sigma）の半分を引いたものとなる。

　Deterministicは「決定論的」と訳されますが、「計算すればいつも同じ値が求められる」といった意味です。といっても、単に計算を行うのではなく、計算の方法を定義したもの（計算グラフと呼ばれます）を作るといった意味です。引数として指定したmu、delta、sigmaや返り値のmu1、mu2は、TensorVariableと呼ばれる型の変数で、いずれも値そのものではなく計算の方法を定義したものです。


AI博士

　ベイズ「t検定」なら、尤度関数としてt分布を使うのでは？　と思われるかもしれませんが、単に、古典的なt検定に対応する手法であるという意味で、ベイズt検定と呼んでいるだけです。分散（標準偏差）の不確実さについては事前分布で定義されているので、母集団が正規分布に従うという前提であれば、尤度関数としても正規分布を素直に使えばいい、ということです。ただし、外れ値がある場合にも、ある程度頑健な分析を行いたいときには尤度関数としてt分布を使うこともあります。


　尤度関数は簡単です。上で作成したmu1とmu2を使って、以下のように記述します。

    pm.Normal("x1_obs", mu=mu1, sigma=sigma, observed=x1_standardized)
    pm.Normal("x2_obs", mu=mu2, sigma=sigma, observed=x2_standardized)

リスト4　尤度関数の定義
運動部の成績の分布は、平均mu1、標準偏差sigmaに従う。観測されたデータ（observed引数）には、標準化されたデータx1_standardizedを指定する。一方の非運動部の成績の分布は、平均mu2、標準偏差sigmaに従う。観測されたデータ（observed引数）には、標準化されたデータx2_standardizedを指定する。

サンプリングを行う

　サンプリングのためのコードはこれまでと同じです。つまり、pymcモジュールのsample関数を使うだけです。これまで断片的に記したコードも含め、結果の要約を表示するところまで、コードの全体を掲載します。サンプルファイルの3番目のコードセルに入力されています。

import pymc as pm
import arviz as az

with pm.Model() as model:
    # 事前分布
    delta = pm.Cauchy("delta", alpha=0, beta=0.707)  # 効果量
    sigma = pm.HalfNormal("sigma", sigma=1)  # 共通の標準偏差（半正規分布とする）
    mu = pm.Normal("mu", mu=0, sigma=1)  # 全体の平均（標準化されているので、平均は0、標準偏差は1としておく）

    # それぞれの平均（効果量の分だけ差があるものとする）
    mu1 = pm.Deterministic("mu1", mu + delta * sigma / 2)
    mu2 = pm.Deterministic("mu2", mu - delta * sigma / 2)

    # オリジナルの平均
    # mu1_orig = pm.Deterministic("mu1_orig", mu1 * x_std + x_mean)
    # mu2_orig = pm.Deterministic("mu2_orig", mu2 * x_std + x_mean)

    # 尤度
    pm.Normal("x1_obs", mu=mu1, sigma=sigma, observed=x1_standardized)
    pm.Normal("x2_obs", mu=mu2, sigma=sigma, observed=x2_standardized)

    # サンプリング
    trace = pm.sample(draws=8000, tune=1000, chains=2, random_seed=42)

# 事後分布の要約
az.summary(trace, var_names=["delta", "mu1", "mu2"])

リスト5　事後分布を求めるためにサンプリングを行う（コード全体）
「# サンプリング」というコメント以降が、これまで断片的に記したコードの続き。後でベイズ因子を求めるために、やや多めに（8000個の）サンプリングを行った。事後分布の要約として、効果量deltaと運動部の平均mu1、非運動部の平均mu2の推定値などを表示する。なお、この平均は標準化された値であることに注意。

　pymc.sample関数によってサンプリングが実行されると、各変数のサンプルが自動的に作られます。上のコードでは、それらをarvizモジュールのsummary関数で要約しています。実行すると、図5のような結果が表示されます。

事後分布の要約
図5 事後分布の要約
効果量δの点推定値は、deltaの行のmean列の値、つまり1.017となる。94％信用区間はhdi_3％〜hdi_97％なので、0.490〜1.547となる。μ1やμ2についても同様に見ていけばよい。ただし、μ1やμ2の推定値は標準化された値によるものであることに注意。r_hatが1.0である（＝収束している）ことも確認しておこう。

　今回はベイズ因子を求めることが目的なので、元の成績について母平均などの推定を行っていません。元の成績の推定値が必要であれば、リスト5の「#オリジナルの平均」の下にある2行のコードのコメント（#）を外した上で、az.summaryのvar_names引数として指定したリストに"mu1_orig"と"mu2_orig"を追加してください。

ベイズ因子を求める 〜 独立した2群のベイズt検定を行う

　今回の主な目的はベイズ因子を求めることです。Savage-Dickey法を使うなら、事前分布（コーシー分布）での、0に対する確率密度と、事後分布をスムーズな曲線にした場合の0に対する確率密度の比を取れば、ベイズ因子BF10の値が簡単に求められます。この連載の第3回で見たコードとほとんど同じですね。異なるのは事前分布としてコーシー分布を使っていることだけです。以下のようなコードになります（サンプルファイルの4番目のコードセルです）。

from scipy.stats import gaussian_kde, cauchy

# 事後分布
delta_posterior = trace.posterior['delta'].values.flatten()
kde = gaussian_kde(delta_posterior)  # 滑らかな曲線にする
posterior = kde.evaluate(0)[0]  # δ=0の位置の確率密度関数の値

# 事前分布（コーシー分布）
prior = cauchy.pdf(0, loc=0, scale=0.707)  # δ=0の位置の確率密度関数の値

# ベイズ因子
bf10 = prior / posterior
print(f"BF10: {bf10:.3f}")

# 出力例：
# BF10: 233.456

リスト6　Savage-Dickey法によりベイズ因子BF10を求める
事後分布をscipy.statsモジュールのgaussian_kde関数を使って、サンプリングされたデータのヒストグラムを滑らかな曲線にする。事後分布と事前分布のδ＝0に対する確率密度関数の値を求め、その比を取ればベイズ因子が求められる。結果は233.456となる。

　上の結果をJeffreysの評価尺度に照らし合わせると、H1を支持する「決定的な証拠となっている」ことが分かります。従って、対立仮説δ ≠ 0（効果量は0ではない）が支持されます。実質的に、μ1の方がμ2よりも大きいと言えるので、やはり、運動部の生徒の方が体力テストの成績がいいということになります。

　参考として、Savage-Dickey法を視覚的に理解できるようなグラフを図6に示しておきます。0.0の位置で、事前分布（prior）の確率密度が、事後分布（posterior）ではぐっと小さくなっていることが分かります。つまり、δ＝0近くの値が得られる傾向が減った、つまり、δ＝0であるとは考えにくい、ということがよく分かります。

ベイズ因子の可視化
図6　ベイズ因子を可視化した例
このグラフを描くためのコードもサンプルファイルに含めてある。0.0の位置では、事前分布の確率密度よりも事後分布の確率密度の方が小さい。データを観測した後では、0.0近くの値が得られる傾向が小さくなったことが分かる。

片側検定を行いたい場合は？

　これまでは、対立仮説をH1:δ ≠ 0としてきました。δの値は正である可能性も負である可能性もあるので、これはいわゆる両側検定に当たります。しかし、μ1>μ2であると考えられる根拠があれば、対立仮説をH1:δ > 0と置くこともできます。これは、片側検定に当たりますね。そのような場合は、以下のような補正を行います。

　BF＋0は、H1：δ > 0という仮説（対立仮説）を支持するベイズ因子を表しています。逆に、BF0＋は、H0：δ ＝ 0という仮説（帰無仮説）を支持するベイズ因子を表し、BF0＋ ＝ 1 / BF＋0という関係になります。


AI博士

　対立仮説がH1：δ < 0である場合は、BF−0がH1：δ < 0という仮説（対立仮説）を支持するベイズ因子を表し、BF0−がH0：δ ＝ 0という仮説（帰無仮説）を支持するベイズ因子を表します。やはり、BF0− ＝ 1 / BF−0です。


　P(δ >0|D)は、データ（D）が得られたときに（つまり事後分布で）δが0より大きくなる確率です。これは、gaussian_kdeクラスのintegrate_box_1dメソッドによって求められます。このメソッドは、gaussian_kdeクラスの曲線を積分するためのメソッドで、引数には下限と上限を指定します。従って、コードは以下のようになります。

p_delta_gt_0 = kde.integrate_box_1d(0, np.inf)
bfplus0 = bf10 * 2 * p_delta_gt_0
print(f"BF+0: {bfplus0:.3f}")

# 出力例：
# BF+0: 466.843

リスト7　片側検定のための補正を行う
リスト6に続けて上のコードを実行すれば、片側検定の結果が得られる。integrate_box_1dメソッドの引数には積分の下限として0、上限としてnp.inf（無限大）を指定している。これで、δ>0である確率が求められる。後は、補正のための式に当てはめるだけ。

　片側仮説の場合、負の部分を考えなくていいので、上に示したような補正を行わず、効果量の事前分布として半コーシー分布を使うという方法もあります。ただし、実際にやってみると、結果がかなり大きな値になってしまいます（一般に、片側検定のベイズ因子は、両側検定の2倍ぐらいの値になります）。これは、サンプリング時に0に近い値があまり生成されないので、事後分布をスムーズな曲線にしたときに、0に対する確率密度関数の値が極めて小さくなるからです。このように、Savage-Dickey法では、効果量がかなり大きい場合には、誤差が大きくなる可能性があります（結論の方向は変わりません）。

標準偏差が異なる場合のベイズt検定を行うには

　これまでは、2群の標準偏差が等しいものとして、平均の差を検定する方法を見てきました。それは、古典的な検定では、ステューデントのt検定に当たるものです。しかし、2群の標準偏差が等しいとは考えられない場合もあるはずです（むしろ、その場合の方が多いかと思います）。こちらは、古典的な検定の「ウェルチのt検定」に当たるものです。その場合、平均的な標準偏差の値を使います。具体的には、以下の式のように標準偏差を調整します。

　それぞれの平均mu1とmu2の計算には、この計算で求めたσを使い、尤度関数には、それぞれの標準偏差sigma1とsigma2を使います。コードの一部分を修正するだけなので、具体的な書き方についてはサンプルファイルのコードを参照してください。解説も加えてあるので、違いが確認できます。

コラム　ベイズ因子が周辺尤度の比で求められる理由

　ベイズ因子とは、データを観測することによって、事前オッズがどれだけ大きく／小さくなったかを表す係数です。ここでの「オッズ」は確率の比を表します。従って、このことは以下のように表されます。

　左辺の分子は、データDが得られたときに仮説H1が正しい確率ですね。左辺の分母は、データDが得られたときに仮説H0が正しい確率です。H1に注目してみると、（3）式の左辺はデータが得られたときにH1がH0に対してどの程度有利であるかを表しています。

　一方、右辺の事前オッズは、データが得られる前にH1がH0に対してどの程度有利であると考えられるかということですね。通常は五分五分であるとするのが普通なので、1:1と考えます。従って、BF10は、データDが得られたときに、このオッズがどの程度拡大するか、つまり、H1がどの程度有利になるかを表す係数になるというわけです。

　ところで、P(H | D)はベイズの定理により以下のように表されます。

　従って、（3）式の左辺は、以下のようにも書けます。

　右辺はP(D)で約分できます。

　（3）式と（6）式を見比べると、

であることが分かります。（7）式の右辺の分母と分子はそれぞれの仮説における周辺尤度です。ベイズ因子が周辺尤度の比であることが分かりますね。ただ、それらがなぜ周辺尤度を表しているのか、という理由はちょっと難しいのですが、尤度関数を事前分布で重み付けして積分し、母数θを式から消した値だからです（このことを周辺化といいます）。一応、式を掲載しておくと、

によって、θが消されるからということです。が、これ以上の深入りは避けることにしましょう。ともあれ、ベイズ因子が周辺尤度の比で表されることについては、（7）式に示した通りです。

　なお、実際に周辺尤度の比を求めて、ベイズ因子を計算するプログラムも、参考としてサンプルファイルに含めてあります。また、第5回で紹介した、数値積分による方法も参考としてサンプルファイルに含めてあります。いずれもサンプルファイル中に簡単な説明を加えてあるので、興味のある方はご参照ください。


ベイズ統計では仮説が正しい確率が求められる！

　ベイズ因子は、仮説が正しい確率ではなく、どちらの仮説がより強く支持されるかという相対的な値でした。例えば、BF10は、仮説H1の方が、仮説H0よりもどの程度強く支持されるかということになります。さらに、ベイズ因子を基に、仮説が正しい確率を求めることもできます。コラムの最初に示した（3）式を使います（式を再掲します。コラムの内容を読んでいなくても大丈夫です）。

　ここで、事前オッズが分かれば、P(H1| D)、つまり「データDが得られたときに仮説H1が正しい確率」を求めることができます。一般に、データを観測する前は、どちらの仮説も確からしさが等しいと考えられるので、事前オッズは1（つまり1：1）とするのが一般的です。その場合、（3）式は、

と書けます。仮説は2つで排反なので（同時には起こらないので）、P(H0 | D)＝1−P(H1 | D)です。これを（9）式に代入すると、

となります。この式をP(H1 | D)について解くと（途中の計算は省きますが）、

となります。リスト6で求めたBF10＝233.456を代入すると、

という極めて高い確率になります。ちなみに、P(H0| D)は、

です。古典的な検定では帰無仮説や対立仮説が正しい確率を求めることはできませんが、ベイズ統計では、帰無仮説が正しい確率や対立仮説が正しい確率も求められるのです！

　今回は、効果量に注目して、独立した2群の平均に差があるかどうかをベイズt検定により分析しました。その中で、Deterministic関数による新たな変数の計算方法の定義や、両側検定と片側検定の取り扱いを見ました。また、コラムではベイズ因子が周辺尤度の比であることを紹介し、最後に、ベイズ因子を基に、仮説が正しい確率を求める方法を解説しました。盛りだくさんな内容でしたが、サンプルプログラムを単に実行するだけでなく、説明と併せて追いかけていくと少しずつ理解が深まっていくと思います。

　次回は、対応のあるデータでの平均の差についてベイズt検定を行う方法を見ていきます。手に汗握るホラー映画を観る前よりも、観た後の方が握力が向上しているのか、という事例を取り扱います。どうぞ、お楽しみに！

今回使った新出のクラス・関数・メソッド（主なもの）

　引数も主なものだけを掲載します。

Cauchy： コーシー分布に従う確率変数を定義するクラス

  • モジュール： pymc
  • 形式： Cauchy(name, alpha, beta)
  • 引数：
    • name： 確率変数に付ける名前を文字列で指定する。この名前でModelに登録される
    • alpha： 分布の中心
    • beta： スケール（分布の広がり）

Deterministic： 計算により新たな確率変数を定義する関数

  • モジュール： pymc
  • 形式： Deterministic(name, var)
  • 引数：
    • name： 確率変数に付ける名前を文字列で指定する。この名前でModelに登録される
    • var： 変数の計算方法
  • 返り値： 作成された確率変数（TensorVariable型の変数）

cauchy.pdf： コーシー分布の確率密度関数の値を求める関数

  • モジュール： scipy.stats
  • 形式： cauchy.pdf(x, loc, scale)
  • 引数：
    • x： 確率変数の値。リストや配列を指定すると複数の値が求められる
    • loc： 分布の中心
    • scale： スケール（分布の広がり）
  • 返り値： xに対する確率密度関数の値

integrate_box_1d： gaussian_kdeクラスの積分値を求めるメソッド

  • モジュール： scipy.stats
  • 形式： gaussian_kde.integrate_box_1d(low, high)
  • 引数
    • low： 積分の下限
    • high： 積分の上限
  • 返り値： 積分値
「やさしい推測統計（ベイズ統計編）」のインデックス

「やさしい推測統計（ベイズ統計編）」

Copyright© Digital Advantage Corp. All Rights Reserved.

[an error occurred while processing this directive]
ページトップに戻る