【Pythonで学ぶデータ分析】母平均と母標準偏差をベイズ推定する～シュークリームの重さは100gと異なるか？：やさしい推測統計（ベイズ統計編）

ある喫茶店のシュークリームの重さを例に、ベイズ統計で「平均」や「ばらつき」をどう推定するのか、さらに基準値と違いがあるかどうかをどう確かめるのかを解説します。『社会人1年生から学ぶ、やさしいデータ分析』ベイズ統計編の第3回です。

» 2026年06月17日 05時00分公開

この記事は会員限定です。会員登録（無料）すると全てご覧いただけます。

　ある喫茶店のシュークリームが、見本よりもはるかに大きい？！　いわゆる「逆詐欺」は本当なのか。連載『社会人1年生から学ぶ、やさしいデータ分析』のベイズ統計編、第3回となる今回は、この商品の重さを手掛かりに、母平均や母標準偏差を推定し、平均が基準値と異なるかをベイズ検定で確かめます。前回は二項分布で比率を扱いましたが、今回は連続的な量を扱う点が新しいところです。

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

連載：

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

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

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

羽山博

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

母平均と母標準偏差をベイズ推定する

　今回も簡単な事例を使って、ベイズ推定／検定の考え方や手順を見ていきます。今回の目標は幾つかのサンプルを基に、商品の母平均と母標準偏差のベイズ推定を行うとともに、母平均がある値と等しいかどうかを検定するというものです。ここでは、ある喫茶店のシュークリームの重量を例に取って見てみます。図1が今回やることのイメージです。データは架空のものです。

図1　シュークリームの重量はどれぐらい？
観測されたデータを基に、平均と標準偏差をMCMC法（マルコフ連鎖モンテカルロ法）によりサンプリングし、母平均μと母標準偏差σを推定するとともに、母平均が100と異なるかどうかを検定する。

　前回までのお話は、クーポンの当選確率や医療におけるタッチングによる癒やしの確率のような「比率」に関するものでした。今回の例は、シュークリームの重量についてのお話なので、「定量的な値」に関するものです。

　例えば、皆さんが、ある喫茶チェーンの監察官で、店舗を視察していたものとしましょう。そのチェーン店ではシュークリームの重量は平均100グラムで、標準偏差は5グラムという規定になっています。しかし、作られたものを手にすると、なんとなくずっしりと重いような気がします。そこで、シュークリームを30個取り出し、重さを測ってみました。観測されたデータは以下の通りです。

data = [106.0, 92.4, 105.9, 100.1, 97.1, 107.6, 108.1, 101.8, 113.3, 107.5,
101.1, 112.0, 106.2, 104.7, 99.6, 104.4, 96.1, 102.0, 110.2, 103.6,
101.8, 106.7, 111.7, 94.2, 103.9, 99.7, 107.3, 103.5, 98.6, 106.6]

　これらのデータを基に、母平均と母標準偏差をベイズ推定してみたいと思います。また、ベイズ因子を求め、モデル比較を行うことにより、平均がある値（ここでは100）と等しくないというモデルがどの程度支持されるかを評価してみます。

　というわけで、今回は以下のような手順で進めていきます。なお、今回の事例（データ）は架空のものです。

（準備として）サンプルデータを基に正規性を確認しておく
測定された値を基に、母平均と母標準偏差をベイズ推定する
Savage-Dickey法によりベイズ因子を求め、母平均が100と等しいかどうかを検定する

サンプルデータを基に正規性を確認しておく

　今回の事例では、シュークリームの重量が正規分布に従っているという前提で、その母平均や母標準偏差を推定します。従って、母集団の正規性を確認しておく必要があります。ここでは、簡易的な方法ですが、ヒストグラムやQ-Qプロットを使って可視化してみましょう（図2）。

図2　ヒストグラムとQ-Qプロット
ヒストグラムはやや右に山が寄っているが、極端な歪（ゆが）みや大きな外れ値はなさそう。曲線はKDE（Kernel Density Estimation：カーネル密度推定）と呼ばれる方法により、分布を滑らかに表したもの。また、Q-Qプロットもデータがほぼ直線上に並んでいるので、正規分布に近いと考えられる。

　ヒストグラムとQ-Qプロットと併せて見ると、ほぼ正規分布と考えてよさそうだということが分かります。ヒストグラムについては、階級の分け方によって見え方がずいぶんと変わってくることに注意してください。KDEによる曲線を表示しておけば、見た目の違いが「緩和」され、分布の形が分かりやすくなります。

　ただし、これらのグラフを描画するためのコードについては、ここでは解説しません（本筋のお話から外れるので）。サンプルプログラムをこちらに置いておきます。リンクを開くとPythonのプログラムが表示されるので、まず、最初のコードセルをクリックし、［Shift］＋［Enter］キーを押してコードを実行してください（グラフに日本語を表示できるようにするための準備です）。続いて、2番目のコードセルを同様に実行するとグラフが表示されます。これらのグラフ作成方法に興味のある方はコード内のコメントに記した簡単な解説を参照してください。

　母集団から取り出したサンプルが正規分布している場合でも、母集団が正規分布であるとは限りません。

　ただ、サンプルサイズがある程度の大きさ（中心極限定理に基づく一般的な目安ではn>30とされることが多い）で、分布の大きな歪みや外れ値がなければ、ある程度頑健であるといわれています（※この30は、今回のシュークリーム30個とは関係なく、数値がたまたま一致しているだけです）。母集団に外れ値があり、裾が重くなっている（中心から離れた値がある程度多い）場合には母集団をt分布と考えて分析したり、母集団に偏りがある場合は対数を取って正規分布に近い形にしてから分析を行うこともよくあります。

　ここでは、母集団を正規分布と見なして分析を進めます。

母平均と母標準偏差を推定する

　さて、ここからが本題です。平均と標準偏差のベイズ推定に移ります。前回、あまり深く考えなくてもいいとお話しした式ですが、念のため掲載しておきます。

　深く考えなくてもいいのは確かなのですが、この式はベイズ統計の「道しるべ」となるので、形だけは気にしておいてください。要するに、事前分布と尤度関数が決まれば、事後分布が求められるということでした。日常の言葉で言うと、それまではこう思っていた（事前分布）が、データが得られた（尤度関数）ので、考えがこう変わった（事後分布）ということです。

　分母の周辺尤度は、全体の大きさを1にするための定数です。前回は、解析的に（数式を計算して）事後分布を求める例を紹介しましたが、今回は、周辺尤度を解析的に求めるのが難しい場合に、シミュレーションを行って事後分布を求める例を紹介します。その場合、周辺尤度の計算をうまく回避できるので、事前分布と尤度関数さえ定義できればいい、ということになります。

　周辺尤度の計算を回避できる理由をごく簡単に説明すると、シミュレーションのためのMCMC（マルコフ連鎖モンテカルロ）法の一つであるメトロポリス・ヘイスティングス法などでは、発生させたサンプルと、次に発生させたサンプルの比を取るので、周辺尤度が約分されてうまく消えるからです。

事前分布を定義する

　この事例では、母平均μと母標準偏差σを推定したいので、以下のような事前分布とします。（1）式では母数をθという文字で代表的に表していますが、この例では、母数はμとσの2つとなっています。

平均の事前分布： 正規分布とする。事前の平均は100、標準偏差は10とする。
標準偏差の事前分布： 半正規分布とする。事前の標準偏差は5とする。

　後でコードの全体像を見ますが、先行して事前分布を定義する部分だけを抜き出してみると、以下のようになります。なんとなく、意味が分かりますね。

mu = pm.Normal("mu", mu=100, sigma=10) # 平均は正規分布
sigma = pm.HalfNormal("sigma", sigma=5) # 標準偏差は半正規分布

　簡単にコードの意味を解説しますが、細かいことは気にせず、写経する（そのまま書き写す）ぐらいの気持ちで眺めてください。pm.Normalはpymcモジュールに含まれる関数で、正規分布をモデル化するのに使われます。同様にpm.HalfNormalは、半正規分布をモデル化するのに使われます。標準偏差については、正の値しか取らないので、事前分布を半正規分布としたわけです（後のコラムでもう少し詳しく説明します）。

　引数について、簡単に説明しておきます。

pymcモジュールのNormal関数：
- "mu"：モデル内部で確率変数を識別するための名前。後で結果を参照するときにも使う
- mu=100：事前分布の平均
- sigma=10：事前分布の標準偏差（平均のばらつきを表す）

pymcモジュールのHalfNormal関数：
- "sigma"：モデル内部で確率変数を識別するための名前。後で結果を参照するときにも使う
- sigma=5：事前分布の標準偏差（標準偏差のばらつきを表す）

　ここで、注意すべき点は、1行目のsigma=10というのは、母平均がmu=100からどの程度ばらつくかという事前の信念を表す値です。データがどの程度ばらつくかということではありません。母平均そのものがどれぐらいの範囲に分布すると思われるか、ということです。ベイズ統計では、母数は変数でしたね。従って、母平均は決まった値ではなく、「動く」ものと考えられます。それがどの程度か、ということです。

　一方、母集団のデータがどの程度ばらつくかという、母標準偏差についての事前の信念は、2行目のsigmaで表されます。そして、その引数のsigma=5というのは母標準偏差そのものがどの程度ばらつくかということです。ベイズ統計では母標準偏差も変数なので、母標準偏差そのものがどの程度ばらつくか、ということです。つまり、母標準偏差の標準偏差ですね。

コラム　標準偏差の事前分布として使われる半正規分布とは？

　半正規分布は、図3のように、正規分布を0で半分に折り返して積み上げたようなものです。

図3　正規分布と半正規分布
半正規分布は正規分布の負の部分を0を軸にして折り返して重ねたようなイメージ。図にはσ＝1の正規分布と半正規分布が含まれているので、それらを見比べると分かりやすい。σが大きくなると、分布が平らになっていく（事前情報が弱くなっていく）ことが分かる。このグラフを描くプログラムもサンプルファイルに含まれている（3番目のコードセル）。

　なお、母数を分散σ²とするなら、逆ガンマ分布が事前分布としてよく使われます。

　話を本筋に戻しましょう。この例では、事前の情報がほとんどありません。そのような場合には、σの値をやや大きくしておき、事前分布の幅を広く取っておくといいでしょう。一方、過去の経験などから、確信を持って見積もられる値があれば、それを指定します。

　なお、ベイズ統計では事前分布を、事前の「信念」と表現することがよくあります。そして、その信念がどの程度強いかを「確信度」と表現することがあります。このように「分析者の判断を反映してよい」という立場から、「主観」という言葉もよく使われます。

　もちろん、これらの用語は個人の勝手な思い込みを表すものではありません。先行研究やこれまでに収集したデータなどに基づく根拠を持って「信じるに値する」といった意味です。

尤度関数を定義する

　上記の事前分布と次に見る尤度関数を指定すれば、事後分布が得られます。尤度関数は、データを観測したときの母数の尤（もっと）もらしさの分布でした（前回お話ししました）。尤度関数の定義についても、その部分だけ抜き出してみます。

obs = pm.Normal("obs", mu=mu, sigma=sigma, observed=data)

　ここでは、尤度関数として、母集団の分布である正規分布を指定します（もし、母集団をt分布とみなすなら、pm.StudentT関数を使いますが、今回はそれについては取り扱いません）。尤度関数の定義では、事前分布の定義と同じpm.Normal関数を使っていますが、事前分布はあくまで母平均の分布を表すもので、尤度関数は観測データを基にした（母数の尤もらしさの）分布であることに注意してください。

　前回の例では、母数の事前分布はベータ分布で、尤度関数は母集団の分布である二項分布（と同じ形の分布）だったので、混同することはなかったかと思いますが、今回は、たまたま平均の事前分布と尤度関数がどちらも正規分布だっただけです。

　尤度関数に含まれる変数を識別するための名前としては"obs"という文字列を指定しています（後で使います）。ここでは、尤度関数にobserved引数が指定されている点に注目してください。observed=に観測データを指定すれば、その観測データを基に尤度関数が定義されることになります。引数の意味を箇条書きにしておきます。

pymcモジュールのNormal関数：
- "obs"：モデル内部で確率変数を識別するための名前。後で結果を参照するときにも使う
- mu=mu：尤度関数の平均
- sigma=sigma：尤度関数の標準偏差
- observed=data：観測されたデータの指定（この指定により、尤度関数とみなされる）

　mu=mu, sigma=sigmaの部分では、引数の右辺の部分に事前分布の返り値を指定していますね。つまり、以下の図4のような関係になっています。

図4　事前分布の定義と尤度関数の定義
引数muや引数sigmaには事前分布の返り値を指定しているが、単なる値を表しているわけではなく、事前分布で定義された「計算のルール」のようなものを表しており、その中の平均や標準偏差を利用できるようになると考えればよい。

　muやsigmaは、プログラミング言語的にはTensorVariableと呼ばれるデータ型なのですが、詳細についてはあらためてお話しします。とりあえずは、mu=には、平均の事前分布を参照する変数（たまたまmuという名前でした）を指定し、sigma=には、標準偏差の事前分布を参照する変数（こちらも、たまたまsigmaという名前でした）を指定する、と覚えておいてください。

ベイズ推定のためのコードを書いてみよう

　事前分布と尤度関数が定義できたので、あとはシミュレーション（サンプリング）を行うだけです。ここからは、コード全体を見た方が分かりやすいでしょう。マルコフ連鎖モンテカルロ法（MCMC）によるサンプリングなどを簡単に行うためのモジュールpymcを使います。以下のコードがサンプルファイルの4番目のコードセルに入力されています。

# リスト1：平均のベイズ推定（PyMCを使った例）
import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt

# データ
data = [106.0, 92.4, 105.9, 100.1, 97.1, 107.6, 108.1, 101.8, 113.3, 107.5,
        101.1, 112.0, 106.2, 104.7, 99.6, 104.4, 96.1, 102.0, 110.2, 103.6,
        101.8, 106.7, 111.7, 94.2, 103.9, 99.7, 107.3, 103.5, 98.6, 106.6]

with pm.Model() as model:

    # 事前分布
    mu = pm.Normal("mu", mu=100, sigma=10) # 平均は正規分布
    sigma = pm.HalfNormal("sigma", sigma=5) # 標準偏差は半正規分布

    # 尤度 (正規分布)
    obs = pm.Normal("obs", mu=mu, sigma=sigma, observed=data)

    # MCMCによるサンプリング
    trace = pm.sample(draws=2000, tune=1000, random_seed=42)

# サンプリング結果（グラフ）の表示
az.plot_posterior(trace)
plt.show()

# サマリーの表示
az.summary(trace)

リスト1　母平均（と母標準偏差）をベイズ推定するためのコード
pm.Model()は、プログラミング言語的には「Modelクラスのインスタンス（実体）を作成している」ということだが、意味としては、サンプリングや結果の記録などに必要な道具を用意する、というイメージ。as modelによりmodelという名前を付け、withの範囲内（字下げされている範囲内）で、この「道具」を利用できるようにしている。実行結果はこの後に掲載する（図5）。

　事前分布と尤度関数の定義についてはすでに見た通りですが、サンプリングや結果の記録を行うために、with pm.Model() as model:の範囲内に書きます。これらを用意しておけば、後はpymcモジュールのsample関数を呼び出してサンプリングを行うだけです。sample関数を抜き出して見てみましょう。

trace = pm.sample(draws=2000, tune=1000, random_seed=42)

　この例であれば、sample関数によりサンプリングが行われ、平均と標準偏差のペアが2000個×2回分作られます（※2回になる理由は後述します）。後で具体例を見ますが、サンプリングされる値は個々のデータではなく、平均や標準偏差です。このことはとても重要なので、絶対に忘れないようにしてください。引数の意味は以下の通りです。

pymcモジュールのsample関数：
- draws=2000： 2000個のサンプルを作成するという意味（既定値は1000個）
- tune=1000：取り出したサンプルの最初の1000個を捨てる（既定値は1000個）
- random_seed=42：乱数の初期設定のための指定（ここでは、42とした）

　draws=2000, tune=1000なら、実際には3000個のサンプルが作られます。ただし、最初の方のサンプルはまだ安定しないので、tune引数で指定した1000個が捨てられ、draws引数で指定した、残りの2000個が使われます。このような操作によって作られた一連のサンプルのことをチェーンと呼びます。

　上のコードには書かれていませんが、チェーンを幾つ作るかはchains引数で指定できます。例えば、チェーンを4回作るなら、chains=4と指定します。詳細については、「今回使った新出のクラス・関数・メソッド」で説明しますが、Google Colaboratoryの、筆者が利用している環境では既定値は2です。従って、今回の例であれば、draws引数で指定した2000個のサンプルが、2回作られることになります（全部で4000個のサンプルとなります）。

　サンプリングには乱数が使われるので、毎回結果が異なります。ここでは、乱数の初期設定を決めておき、毎回同じ結果が出るようにしています（動作確認などに便利です）。

　乱数の初期設定に使った42という値には、特に意味はありません（どんな値でも構いません）。有名な話ですが、42という値はアダムス著の『銀河ヒッチハイク・ガイド』（1979年）というSF小説の第27章に登場する興味深い（?）数値にちなんだものです（ネタバレになるので、これ以上は説明しません）。この原稿の執筆時点では、Amazonプライムで吹き替え版の映画を観ることもできます（1時間7分あたりにそのシーンがあります）。

母平均と母標準偏差の推定結果を見てみよう

　サンプリング結果としては、サンプリングされた平均と標準偏差以外にも、サンプリングの履歴などに関するさまざまな情報が含まれます。上のコードでは、結果をtraceという変数で参照できるようにしているので、それを使えば可視化や結果の要約などができます。その際、可視化や結果の要約のためにArviZと呼ばれるモジュールが利用できます。例えば、plot_posterior関数でサンプリング結果を可視化でき、summary関数で結果の要約が表示できます（図5）。

図5　サンプリング結果の可視化と要約
上の2つのグラフがサンプリングされた平均の分布と標準偏差の分布。ヒストグラムの形ではなく、滑らかな曲線で表示されている。下の表がsummary関数の出力。結果の見方は以下で解説する。

　plot_posterior関数にサンプリングされた結果（trace）を指定すれば、事後分布のグラフが表示できます。グラフには平均の点推定値と信用区間も表示されていますが、詳細な値については、summary関数の出力で確認できます。なお、特定の変数のみのグラフだけを描くのであれば、az.plot_posterior(trace, var_names=["mu"])のように、var_names引数に変数名のリストを指定します。変数名は、事前分布を定義する際に最初の引数に指定した文字列（"mu"や"sigma"でしたね）で識別されます。

　summary関数の出力を見てみましょう。「mean」の列が点推定値（平均）となっています。母平均の点推定値が103.766、母標準偏差の点推定値が5.251となっていることが分かります。「sd」の列は作成されたサンプルの標準偏差を表します。サンプリングされた平均そのものや標準偏差そのもののばらつきですが、とりあえず、それほど気にする必要はありません。

　信用区間は、「hdi_3%」の列の値から「hdi_97%」の列の値です。母平均については、101.999 ≤ μ ≤ 105.506であることが分かります。既定の設定では94％信用区間となっていることにも注目してください（※グラフ中央の「94% HDI」がこれに当たります。HDIは信用区間を表す指標）。これは、頻度論で一般的に使われる95％信頼区間とは異なる、という意味合いを強調するためです。95％信用区間を表示したい場合は、summary関数の引数にhdi_prob=0.95を追加します。

　ただし、この信用区間の意味は要注意です。例えば、101.999 ≤ μ ≤ 105.506というのは、データがその範囲に94％の確率で収まっているということではなく、母平均が94％の確率でその範囲にあるということです。同様に、σの信用区間が4.028 ≤ σ ≤ 6.596となっているのも、母標準偏差が94％の確率でその範囲にあるということです。

　上記以外で、summary関数の出力で注目すべき列の簡単な解説を箇条書きにしておきます。

mcse_mean：平均の標準誤差。小さいほどサンプリングが安定している
ess_bulk：分布の中心付近の有効サンプルサイズ。大きいほどサンプリングが信頼できる
r_hat：収束診断値。1に近ければ分布が収束している（安定している）。1.1より大きい場合は十分に収束していないので、drawsやtuneの値を増やす必要がある

サンプリングされた値を見てみよう

　結果が得られたところで、いったん立ち止まって、サンプリングされた平均と標準偏差の値を実際に見てみましょう。すでにお話しした通り、サンプリングされた値は、個々のデータ（個々のシュークリームの重さ）を表す値ではありません。平均や標準偏差の値が幾つも作られているということです。それを実際に確認しておこうというわけです。以下のコードを実行すると、実際にサンプリングされた平均と標準偏差の値が表示できます（サンプルファイルの5番目のコードセルに入力されています）。

print(trace.posterior['mu'].values.flatten()) # サンプリングされた平均
print(trace.posterior['sigma'].values.flatten()) # サンプリングされた標準偏差

リスト2　サンプリングされた事後分布のデータを表示する
trace.posterior['変数名']と書けば、その変数の事後分布に関する情報が得られる。valuesは、その中に含まれる値そのもので、チェーン数×サンプル数の配列になっている。今回の例だと、2000個のサンプルを2チェーン作成したので、2×2000の配列。それをflattenメソッドによって一次元に並べて表示する。

　以下に示すのは、リスト2を実行した結果を見やすく書き直したものです。上段が平均、下段が標準偏差のサンプリング結果です。


変数＼インデックス	0	1	2	……	3998	3999
平均（"mu"）	103.78	103.75	104.05	……	103.46	103.17
標準偏差（"sigma"）	5.356	5.223	4.841	……	5.818	4.753
表1　サンプリングされた平均（mu）と標準偏差（sigma）の値

　しつこいようですが、これらの値は個々のデータを表す値ではなく、いわば、母数の候補としてサンプリングされた幾つもの平均と標準偏差の値です。図5のサマリーに表示されていた母平均の推定値103.766は、これらのサンプル（上段の103.78, 103.75, 104.05, ..., 103.46, 103.17）の平均です。同様に、図5のサマリーに表示されていた母標準偏差の推定値5.251は、下段のサンプルの平均です。

　次は、母平均が100と等しいかどうかを調べるためにベイズ検定を行うのですが、もうちょっと立ち止まっておさらいをしておきましょう。以下のコラムは、事前分布、尤度関数、周辺尤度についてもう少し実感を持っていただくための簡単な例です。一刻も早くベイズ検定を行いたいという方は読み飛ばしていただいても構いません。

コラム　ベイズの定理をもっと簡単な例で確認する

　以下の式に含まれる各項の詳しい意味は前回のコラムで解説しました。しかし、まだ実感が湧かない、という方も多いのではないかと思います。

　そこで、極めて簡単な例でこの式の計算を行ってみましょう。例えば、事象（確率変数）が、A、B、Cの3つしかないものとします。例えば、Aが負け、Bが勝ち、Cが引き分け、などです。

　データを観測する前は、これらの値が0.2, 0.6, 0.2という確率だと見積もられたものとします。これが事前分布ですね。事前分布は確率分布なので合計すると1になります。

　次に実際にデータを観測したところ、尤度が0.1, 0.8, 0.1だったとします。これが尤度関数です。実際にはB（勝ち）がよく起こったわけですね。

　表にまとめながら見ていきます。表中の「事前分布」と「尤度関数」の値はすでに示した通りです。

図6　事前分布・尤度関数・周辺尤度・事後分布を理解する
ごく簡単な例で、事前分布、尤度関数から周辺尤度を計算し、さらに事後分布を求める。具体的な計算方法については、この後説明する。

　（1）式の分子を求めてみます。事前分布と尤度関数のそれぞれの値を掛けると、

という分布になります。つまり、分子は0.02, 0.48, 0.02という分布です。表ではブルーの枠で囲んだ（X）の部分に当たります。ただし、これは確率分布ではないので、合計は1になりません（0.02＋0.48＋0.02＝0.52ですね）。

　次に分母です。周辺尤度は、事前分布で重み付けした尤度の合計です。つまり、

となります。表ではオレンジの枠で囲んだ（Y）の部分です。この値、つまり周辺尤度は母数（パラメーター）の全体的な尤もらしさを表す定数です。

　最後に、事後分布です。これは簡単ですね。分子／分母なので、

となります。A（負け）の確率が0.038、B（勝ち）の確率が0.923、C(引き分け)の確率が0.038となりました。B（勝ち）の確信度がより高くなったというわけですね。なお、事後分布も確率分布なので、全体では0.038＋0.923＋0.038＝0.999 ≈ 1となります（ぴったり1にならないのは小数点以下の誤差によります）。

　実際には、尤度関数や事前分布が、今説明したような離散型の分布ではなく、連続型の分布になっていることが多いので、周辺尤度は合計ではなく、積分した値になっています。先ほど事後分布の合計がほぼ1になったのは、分子を周辺尤度（0.52）で割ったからです。このように、周辺尤度は事後分布全体を1にするための定数だと分かりますね。

　「周辺」の意味についても、表の中の値を集計して、その外側（＝周辺）に書き出した合計値、といったイメージで理解できると思います。先ほど各候補（A・B・C）の「事前分布×尤度」を全部足して0.52を求めましたが、この合計が「周辺」尤度です。候補ごとの値を合計して表の端に出すので「周辺」と呼ばれる、と考えればよいでしょう。

母平均が100と等しいかどうかをベイズ的に検定する

　前回お話ししたように、古典的な検定に対応するものとして、ベイズ因子（ベイズファクター）によるモデル比較が利用できます。ベイズ因子は2つの仮説（モデル）の周辺尤度の比で求められます。上のコラムで、周辺尤度は「母数（パラメーター）の全体的な尤もらしさを表す定数」であるというお話をしましたが、複数のモデルを設定し、周辺尤度の比を取れば、全体としてどちらがより尤もらしいか、ということが分かるわけです。

　しかし、周辺尤度を求めるのはちょっと面倒です。そこで、Savage-Dickey法を利用して、ある点（例えば、100）における事前分布と事後分布の確率密度関数の値の比を求めます（それが、周辺尤度の比と等しくなることの証明はかなり難しいので、ここでは触れません）。

　ここでの帰無仮説と対立仮説は以下のようになります。

帰無仮説H₀： μ ＝ 100
対立仮説H₁： μ ≠ 100

　対立仮説は、古典的な仮説検定の書き方に合わせて上のように記されるのが一般的ですが、正確には「μが100という固定値とは異なる」ということではなく「μは平均100、標準偏差10の正規分布に従う」ということです（100と10は、μの事前分布を定義したときの値で、具体的には前半の「事前分布を定義する」でmu = pm.Normal("mu", mu=100, sigma=10)と設定した値です）。対立仮説をわざわざ分布の形で表すことにちょっと違和感があるかもしれませんが、これはμを「100ぴったり」と決めつけず、平均100あたりで標準偏差10程度の不確実性を持つ（＝μに幅を持たせて表している）と考えているということです。ただし、事後分布ではμの平均は103.766、標準偏差は0.946（summaryのsd欄に表示されている値）に更新されます。

　以下のコードがサンプルファイルの5番目のセルに入力されています。リスト1に続けて実行すれば、ベイズ因子BF₁₀の値が表示されます。

# Savage-Dickey法による1群の平均のベイズ検定（μ=100の位置での確率密度関数の比を求める）
from scipy.stats import gaussian_kde, norm

# 事後サンプル
mu_samples = trace.posterior["mu"].values.flatten()

# KDEによる事後分布の確率密度
kde = gaussian_kde(mu_samples) # 滑らかな曲線にする
posterior_density = kde.evaluate([100])[0] # 100の位置での確率密度

# 事前分布の確率密度
prior_density = norm.pdf(100, loc=100, scale=10) # 100の位置での確率密度

# ベイズ因子
BF10 = prior_density / posterior_density # 事前/事後でBF10を求める

# 結果はBF10の値が68.6となる。Jeffreysの評価尺度によると
# 対立仮説H1：μ≠100を支持する「非常に強い証拠」となっている。
print("BF10:", BF10) # 68.60016058862766という結果が表示される

リスト3　ベイズ因子を求めて、ベイズ仮説検定を行う
事後サンプルの取り出し方は、リスト2で見た通り。続いて、scipy.statsモジュールのgaussian_kdeクラスを利用し、度数分布を滑らかな曲線にする。その100の位置にある値が事後分布の確率密度関数の値。gaussian_kdeクラスのevaluate関数を使って取り出す。引数にリストを指定すれば複数の値を取り出せるので（ここでは[100]しか指定していないが）、返り値の[0]番目だけを取り出す。事前分布の確率密度関数の値は、おなじみの正規分布（scipy.statsモジュールのnorm関数）に100を指定し、平均100、標準偏差10の場合の値を取り出す。後は割り算をするだけ。

　結果はBF₁₀＝68.60となりました。前回掲載したJeffreysの評価尺度によると対立仮説H₁を帰無仮説H₀よりも支持する「非常に強い証拠」となっていることが分かります。

　ここまで、いったい何の話だったか忘れてしまった方もおられるかもしれませんが、喫茶店でのシュークリームの話でしたね。上の結果から、シュークリームの重さの平均は100グラムとは異なる（100グラムよりは重い）といえることが分かりました。

　もちろん、監察官としての対応は「そのままにしてお客さまに喜んでもらおう。いや、もっと盛っていい」となるか「コストダウンのために、平均がきちんと100グラムになるようにしよう」となるか、（ここでは分析を行っていませんが）「標準偏差を小さくして不公平感をなくそう」となるかは、考え方次第です。

　Savage-Dickey法が適用できる前提については前回お話ししました。仮説がネストしており、対立仮説の事前分布が連続型で、パラメーターが1つであるといった制約でした。ここでは、パラメーターはμとσですが、検証したいパラメーターがμで、事前分布が独立して定義されていれば問題ありません。

　「仮説がネストしている」というのは、帰無仮説H₀が対立仮説H₁の特殊な値であるということでしたね。対立仮説は、正確には「μは平均100、標準偏差10の正規分布に従う」でした。帰無仮説の「平均が100である」というのはその中の1点なので、前提を満たしています。

　ただし、KDEによる事後分布の確率密度については近似値なので、結果もあくまで近似的なものです。そのため、サンプリングの個数を増やして評価するのが一般的です。例えば、リスト1のpm.sample関数でdraws=8000としてサンプリングを行い、リスト3を実行すると、BF₁₀＝79.45という値が得られます（得られる値は異なりますが、分布によほどの歪みがない限り、大勢に影響することはありません）。

　Savage-Dickey法が適用できない場合にはSMC法（Sequential Monte Carlo法：逐次モンテカルロ法）を使って周辺尤度の比を求めるなどの方法もありますが、今回は解説しません（サンプルファイルにはその例も含めてあります）。

　なお、ここでは、平均がどの程度かという推測は行いましたが、実際にシュークリームを作成したときに、商品の重さがどの範囲に入っているか（例えば、95～105グラムに入る確率はいくらか）ということは、この段階では求められません。ここで求めたのはあくまでも平均がどのように分布しているかということで、個々の商品の重さを予測しているわけではありません。例えば、平均の94％信用区間は母平均がその中にある確率が94％であるということです。個々の商品の重さを測ったとき、94％の確率で信用区間に入っているということではありません。

　確かに、平均の点推定値が103.766、標準偏差の点推定値が5.251なので、これらの値を使って、これから作成するシュークリームが95～105グラムである確率を簡易的に求めることもできそうではあります。しかし、これらの点推定値を使うと、平均や標準偏差の不確実さが十分に反映されません。そういった不確実さも反映した予測を行うためには、サンプリングされた数多くの平均と標準偏差を基に事後予測分布を作成し、それを基にある範囲内に入っている確率を求める、という方法を使います。……が、今回はそろそろお腹いっぱいかと思います。その話は次回に回すこととします。

　今回は、正規分布の母数（平均と標準偏差）のベイズ推定を行い、平均のベイズ検定を行いました。今のところは、理屈を全て理解することよりも、慣れることの方が大事なので、単にサンプルコードを実行するだけでなく、ぜひ写経して実際に実行してみてください。

　さて、次回は、サンプリングされた平均と標準偏差を基に事後予測分布を作成し、工作機械によって製作される製品が仕様を満たしているかどうかを判定し、規格外品の廃棄コストを見積もります。どうぞ、お楽しみに！

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

　引数も主なものだけを掲載します。引数に「=値」と書かれたもの（例：=True）は既定値を表します。

Model：確率変数、事前分布、尤度関数、観測値などを管理するためのクラス

モジュール： pymc（実際にはpymc.model.coreに含まれる）
形式： Model()

　※通常は、以下のように使う

import pymc as pm
with pm.Model() as モデル名:
    # 事前分布の定義
    # 尤度関数の定義（観測データを指定）
    # サンプリングの実行

Normal：正規分布に従う確率変数を定義するクラス

モジュール： pymc
形式： Normal(name, mu, sigma, observed)
引数：
- name：確率変数に付ける名前を文字列で指定する。この名前でModelに登録される
- mu：平均
- sigma：標準偏差
- observed：観測データ。この引数の指定がない場合は事前分布として扱われ、この引数の指定がある場合は尤度関数として扱われる

HalfNormal：半正規分布に従う確率変数を定義するクラス

モジュール：pymc
形式： HalfNormal(name, sigma, observed)
引数：
- name：確率変数に付ける名前を文字列で指定する。この名前でModelに登録される
- sigma：標準偏差
- observed：観測データ。この引数の指定がない場合は事前分布として扱われ、この引数の指定がある場合は尤度関数として扱われる

sample： MCMCサンプリングを行う関数

モジュール： pymc
形式： sample(draws=1000, tune=1000, chains=None, cores=None, random_seed=None, progressbar=True)
引数：
- draws：最終的に取り出すサンプルの個数
- tune：ウォームアップ（バーンイン）の個数。ここで指定した個数の最初のサンプルはまだ安定していないので捨てる
- chains： MCMCチェーンの数。drawsとtuneで指定したサンプリングを何回行うか。None（既定値）の場合は、cores引数で指定された値と2のうち大きい方が自動的に設定される
- cores：並列的に実行されるチェーンの数。None（既定値）の場合は、CPUのコア数と4の小さい方が自動的に設定される
- random_seed：乱数の初期設定。値を指定しない場合は毎回異なる乱数となる
- progressbar：サンプリングの進行状況を表示する場合はTrueを、表示しない場合はFalseを指定する
返り値：サンプリングの履歴（サンプリングされた値や診断結果など）。通常、arvizによって分析したり要約したりする

plot_posterior：事後分布のグラフを描く

モジュール： arviz
形式： plot_posterior(data, var_names=None)
引数：
- data： arviz.InferenceData形式のデータ（通常、pymc.sampleの返り値を指定する）
- var_names：グラフを描きたい変数名のリスト

summary：分析結果についての要約統計量を求める

モジュール： arviz
形式： summary(data, var_names=None, hdi_prob=None)
引数：
- data： arviz.InferenceData形式のデータ（通常、pymc.sampleの返り値を指定する）
- var_names：要約統計量を表示したい変数名のリスト
- hdi_prob：信用区間の範囲を指定する。既定値は0.94