【Pythonで学ぶデータ分析】母平均のベイズ推定と予測サンプルの作成 〜 規格外製品の廃棄コストを見積もる：やさしい推測統計（ベイズ統計編）
「不良品で幾ら損する？」をベイズ統計で見積もってみましょう。製品のサイズを測ったデータから、平均やばらつきを推定し、さらに「これから作る製品が規格外になる確率」までをPythonを使って予測します。『社会人1年生から学ぶ、やさしいデータ分析』ベイズ統計編の第4回です。
この記事は会員限定です。会員登録（無料）すると全てご覧いただけます。
規格外製品の廃棄コストはどれぐらい？ 連載『社会人1年生から学ぶ、やさしいデータ分析』のベイズ統計編。今回は第4回です。製作された製品のサイズを基に、母平均と母標準偏差を推定し、将来製作される製品のサイズを予測します。それにより、規格外となる確率を求め、規格外製品の廃棄コストを見積もろうというわけです。前回は母平均と母標準偏差の推定や、ベイズ検定を扱いましたが、今回はさらに予測サンプルを作成します（事後予測分布が作成できます）。
この連載では、簡単な事例を見ながら、ベイズ統計の考え方と分析の進め方を解説します。新しい用語や考え方が幾つも出てきますが、全てを理解しなくても大丈夫です。登場するたびに、分かりやすく丁寧に説明するので、その時点で分からないことがあっても気にせず、先に進んでください。回を追うごとに、少しずつ理解が深まります。どんな考え方で、どんな手順で分析を進めるのか、といった大きな流れを捉えてください。
連載：
この連載では、データをさまざまな角度から分析し、その背後にある有益な情報を取り出す方法を学ぶ『社会人1年生から学ぶ、やさしいデータ分析』シリーズの「記述統計と回帰分析編」「確率分布編」「推測統計（区間推定編・仮説検定編）」に続く「ベイズ統計編」です。
これまでの推測統計を土台に、近年活用が広がっているベイズ統計の考え方と分析手順を、古典的な手法との違いを整理しながら解説します。初めての方でも無理なく理解できるよう具体例を通して進めるとともに、ベイズ的なアプローチの特徴やメリットを実感できる構成とし、「どのように考え、どう使い分けるのか」に重点を置いて解説していきます。
羽山博筆者紹介： IT系ライターの傍ら、かつて、非常勤講師として東大で情報・プログラミング関連の授業を、一橋大でAI関連の授業を担当。まだ発展途上のようだが、生成AIへの質問と回答を分野ごとにまとめ、体系的に整理、補足しつつ、さらに個別の学びを促す記述を加えた資料を一貫して作成できるツールが実用レベルで本格的に使えるようになると、書籍の概念が変わるのではないかと思ったりしている。一斉授業だったものが、ユーザーに合った個別指導のできる書籍といったイメージ。そうなると、出版社のビジネスモデルも大きく変わり、私の仕事もなくなってしまいそうだけど、まあ、そのときは、縁側でお茶でもすすりながら、庭の景色を眺めて暮らすことにしようと思う。ウチには縁側も庭もないけど。
規格外製品の廃棄コストはどの程度か
今回も簡単な事例を使って、ベイズ推定の考え方や手順を見ていきます。今回は、少し応用に重点を置きます。図1で今回やることのイメージを確認しておきましょう。
図1 ある工作機械での規格外製品の廃棄コストはどれぐらいか今回のテーマは2つ。観測されたデータを基に、母平均と母標準偏差を推定することと、推定された母平均と母標準偏差を使って予測サンプルを作り、規格外製品となる確率を求め、規格外製品の廃棄コストを見積もること。
まず、工作機械によって製作した製品の母平均と母標準偏差をベイズ推定します。これは前回のおさらいです。続いて、それを基に予測サンプルを作り、規格外製品となる確率を求めます。その確率が求められれば「製作した個数×規格外製品となる確率×1個当たりの廃棄コスト」で全体の廃棄コストが求められるというわけです。
この事例では、製品のサイズが15±0.5（ミリ）の範囲に収まっていれば、規格通りに製作されているものと見なします。では、母平均と母標準偏差のベイズ推定から始めましょう。前回の内容は十分理解できた、という方は、以下のお話については斜め読み程度でも構いません。事例と値は前回のものと異なりますが、問題の構造は全く同じです。
製品の母平均と母標準偏差をベイズ推定する
工作機械で製作された製品を30個抜き出してサイズを測定してみたところ、以下のような値であったものとします（これらのデータは架空のものです）。
data = [14.95, 14.94, 15.00, 15.14, 14.90, 15.22, 15.05, 15.08, 15.08, 14.92,
14.85, 15.02, 15.06, 14.74, 15.01, 15.05, 14.82, 14.80, 15.12, 15.13,
15.00, 14.94, 14.86, 14.92, 15.07, 14.95, 15.02, 14.99, 15.16, 14.86]
母集団は正規分布に従っているものとして話を進めます。正規性の確認については、前回やった通りなので、ここでは省略します。こちらに置いたサンプルプログラムを実行すると、正規性を確認するためのヒストグラムとQ-Qプロットが表示されます。リンクを開くとPythonのプログラムが表示されるので、まず、最初のコードセルをクリックし、［Shift］＋［Enter］キーを押してコードを実行してください（グラフに日本語を表示できるようにするための準備です）。続いて、2番目のコードセルを同様に実行するとグラフが表示されます。
ベイズの定理により事後分布を求めるための、おなじみの式を以下に記しておきます。私たちが知りたいのは、左辺の事後分布でしたね。ここでは、母数をまとめてθと表記していますが、今回の例であれば、母数は平均μと標準偏差σです。
この式の意味については前々回と前回で詳しく説明しました。「忘れてしまった」「よく分からなかった」という方もご心配なく。具体的な操作を行う中で少しずつ理解が深まるので、焦ることはありません。このまま先に進んでもらっても構いませんが、一応、それぞれの項について簡単にまとめておきます。
- 左辺
- 事後分布 π(θ | x)： データxが得られたときの母数θの分布←これを求めたい
- 右辺
- 尤度関数 L(x | θ)： ある母数θを決めたとき、データxが得られることがどの程度もっともらしい（尤もらしい）かを表す関数。xが変数ではなく、θが変数であることに注意
- 事前分布 π(θ)： データが得られる前に想定される母数θの分布
- 周辺尤度 m(x)： 事前分布で重み付けした尤度の合計（加重平均）。定数となる ← 今のところこれは気にしなくてもいい
MCMC（マルコフ連鎖モンテカルロ法）のメトロポリスヘイスティングス法などでは、周辺尤度を計算しなくても、事後分布が求められます。従って、尤度関数と事前分布さえ分かれば、事後分布が求められるということでした。
この事例では、前回同様、平均μと標準偏差σの推定を行いたいので、以下のような事前分布を設定することにします。
- 平均の事前分布： 正規分布とする。事前の平均は15、標準偏差は1とする
- 標準偏差の事前分布： 半正規分布とする。事前のスケール（広がりを決める値）は1.2とする
平均については、機械の仕様がそもそも平均15ミリなので、それについてはある程度信ぴょう性があるものと考えて事前の平均を15としています。ただし、性能のばらつきについては不明なので、広く取るために事前の標準偏差を1としています（サンプルの標本標準偏差は0.11なので、かなり広めに取ってあります）。
標準偏差の事前分布は、母標準偏差そのものにどの程度の幅があると考えられるかということです。いわば母標準偏差の標準偏差です（ベイズ統計では、母数を変数として考えるので、母標準偏差も「動き」ます。従って、母標準偏差の標準偏差を考えるわけです）。標準偏差の値は負になることはないので、半正規分布（負の部分を折り返して正の値だけを取るようにした正規分布）としています。その幅についても不明なので、広く取るために1.2としています。
コードを見てみましょう。サンプルファイルの3番目のコードセルに入力されています。考え方は前回のコードと全く同じですが、今回はトレースプロットを表示することにします。トレースプロットとは、チェーン（サンプリングされた一連のデータ）を先頭から順に、つまり時系列でプロットしたグラフです。ただし、以下のコードはpymcのバージョン5.28.5、arvizのバージョン0.22.0で作成したものです。バージョンが異なると正しく実行できないことがあります。サンプルファイルには、pymc6.0.1、arviz1.2.0で実行できるコードも4番目のセルに含めてあります（arviz1.2.0では信用区間の既定値もhdi《最高密度区間》からeti《等尾区間》に代わり、その確率も89％になっています）。
# リスト1：平均と標準偏差のベイズ推定（PyMCを使った例）
import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt
# データ
data = [14.95, 14.94, 15.00, 15.14, 14.90, 15.22, 15.05, 15.08, 15.08, 14.92,
14.85, 15.02, 15.06, 14.74, 15.01, 15.05, 14.82, 14.80, 15.12, 15.13,
15.00, 14.94, 14.86, 14.92, 15.07, 14.95, 15.02, 14.99, 15.16, 14.86]
with pm.Model() as model:
# 事前分布（値は異なるが前回のコードと考え方は同じ）
mu = pm.Normal("mu", mu=15, sigma=1) # 平均は正規分布
sigma = pm.HalfNormal("sigma", sigma=1.2) # 標準偏差は半正規分布
# 尤度 (正規分布)
obs = pm.Normal("obs", mu=mu, sigma=sigma, observed=data)
# MCMCによるサンプリング
trace = pm.sample(draws=2000, tune=1000, chains=2, random_seed=42)
# トレースプロットの表示（前回のコードと異なる部分）
az.plot_trace(trace, chain_prop={"color": ["C0", "C1"]}, figsize=(12,6))
plt.show()
# サマリーの表示
az.summary(trace)
データと事前分布が異なるだけで、それ以外は前回のコードと全く同じ。ただし、結果として、トレースプロットとサマリーを表示するようにした。トレースプロットはチェーンごとに色を変えてある。
前回は、ArviZと呼ばれるモジュールを利用して、plot_posterior関数でサンプリング結果を可視化しましたが、今回は、plot_trace関数を使ってトレースプロットを表示しました。plot_trace関数に指定したchain_prop引数では、チェーンの表示方法などを辞書形式で指定できます。ここでは、色として、既定の色を表すC0（通常はブルー）とC1（通常はオレンジ）に指定しています。図2にトレースプロットと、summary関数で得られた要約とを併せて示します。
図2 トレースプロットと結果の要約左の2つのグラフがサンプリングされた平均の分布と標準偏差の分布。ヒストグラムの形ではなく、滑らかな曲線で表示されている。右の2つのグラフが、サンプリングの過程を可視化するトレースプロット。偏りがなく「毛虫」のように表示されていれば、うまくサンプリングができている。下の表がsummary関数の出力（結果の見方は以下で解説）。
先に、summary関数の出力を見ておきます。「mean」の列が平均（期待値）の点推定値です。母平均の点推定値が14.988、母標準偏差の点推定値が0.121となっていることが分かります。信用区間は「hdi_3％」の列の値から「hdi_97％」の列の値です。母平均の94％信用区間は14.944 ≤ μ ≤ 15.028であることが分かります。
上で述べたようにトレースプロットは、サンプリングされた値を順に（時系列的に）プロットしたものです。中心付近に値が集中して塗りつぶされ、そこから外れた値も幾らかあるような「毛虫」のように表示されていれば、サンプリングがうまくいっていることが分かります。このことは、summary関数の出力を基に数値で確認できます。
- mcse_mean： 平均の標準誤差。小さいほどサンプリングが安定している
- ess_bulk： 分布の中心付近の有効サンプルサイズ。大きいほどサンプリングが信頼できる
- r_hat： 収束診断値。1に近ければ分布が収束している（安定している）。1.1より大きい場合は十分に収束していないので、drawsやtuneの値を増やす必要がある
ここに来て、ちゃぶ台返しのような話ですが、実は、正規分布の場合、母分散の事前分布を逆ガンマ分布とし、母平均の事前分布をその分散に依存する正規分布とすれば、母平均の事後分布も正規分布になり、PyMCによってサンプリングを行わなくても解析的に母平均の推定ができます。連載の第2回でも触れましたが、このように、事後分布が同じ形になる事前分布を共役事前分布と呼びます。今回は、前回に引き続きPyMCによるサンプリングの流れでお話を進めているので、共役事前分布を利用する例については割愛します。……が、その方がはるかに簡単なので、次回の番外編で補足として触れることとします。
コラム トレースプロットを自分で描いてみる
サンプリングされた値を自分で順にプロットすれば、トレースプロットが描けます。実際にやってみると、トレースプロットがどのようなものであるかが実感できると思います。サンプリングされた値を実際に表示するためのコードは前回紹介したので、それを少し応用すれば簡単です。コードはサンプルファイルの5番目のコードセルに入力されています。
# リスト2：サンプリングされた平均と標準偏差を表示して、トレースプロットを自分で描いてみる
mean_samples = trace.posterior['mu'].values.flatten() # サンプリングされた平均
mean_stds = trace.posterior['sigma'].values.flatten() # サンプリングされた標準偏差
# サンプリングされた値を表示する
print(mean_samples.round(3)) # 小数点以下3桁まで
print(mean_stds.round(3))
# トレースプロット（平均のみ）を自分で作成する
plt.figure(figsize=(6, 3))
plt.title('mu')
plt.xlim(0, 1999)
plt.plot(mean_samples[0:2000], color="C0", alpha=0.3) # 最初のチェーン
plt.plot(mean_samples[2000:4000], color="C1", alpha=0.3) # 2番目のチェーン
plt.show()
# 出力例（サンプリングされた値）：
# [14.997 14.993 14.996 ... 14.986 15.006 15.021]
# [0.13 0.145 0.098 ... 0.113 0.128 0.128]
サンプリングされた平均は14.997, 14.993, 14.996 ... という値。flattenメソッドで1次元に並べたので、前半の2000個が最初のチェーンの値、後半の2000個が2番目のチェーンの値。それらを折れ線グラフとしてプロットすれば、トレースプロットになる（図2の右上とほぼ同じグラフが表示される）。また、サンプリングされた値をヒストグラムにし、それを基にスムーズな折れ線を作成すれば、図2の左上と同様のグラフになる。
規格外製品の廃棄コストを計算する 〜 予測サンプルを基に規格外の確率を求める
続いて、規格外製品の廃棄コストを見積もってみます。これが今回の主要な目標でしたね。この例では、許容できる誤差を0.5（ミリ）とするので、製作される製品のサイズが14.5より小さいものと15.5より大きいものが規格外製品となります。
規格外製品となる確率が求められれば、廃棄コストも計算できます。では、その確率をどのようにして求めればいいのでしょうか。単純に考えると、事後分布の平均と標準偏差を基に、14.5より小さい確率と15.5より大きい確率を求めれば、規格外になる確率が求められると思われるかもしれません。
しかし、そうはいきません。そのようにして求めた確率は、あくまでも母平均が14.5より小さい確率と15.5より大きい確率です。やるべきことは、これから製作される製品のサイズを予測し、それが14.5より小さい確率と15.5より大きい確率を求めることです。
そこで、まず、サンプリングされた2×2000個の平均と標準偏差（リスト1のプログラムで求められた値＝コラムのリスト2で実際に表示した値）を使って、製品のサイズとして予測される値（予測サンプル）を作成します。といっても、自分で予測サンプルを作る必要はありません。sample_posterior_predictive関数を使えば簡単に求められます。その際、予測値は観測されたデータと同じ個数、つまり30個ずつ作られます。つまり、
- 0番目の（14.997, 0.13）という平均と標準偏差を使って、30個のサンプルを作成する
- 1番目の（14.993, 0.145）という平均と標準偏差を使って、30個のサンプルを作成する
- 2番目の（14.996, 0.098）という平均と標準偏差を使って、30個のサンプルを作成する
……
- 3999番目の（15.021, 0.128）という平均と標準偏差を使って、30個のサンプルを作成する
というように、「2チェイン」×「2000個の平均と標準偏差」×「30個の予測サンプル」＝合計120,000個の予測サンプル（将来製作されるであろう製品のサイズを予測した値のサンプル）が作成されるわけです。次に、それらの値が規格外になる確率を求めます。全ての値について14.5より小さいか、15.5より大きいかを調べるだけでできます。従って、以下のようなコードを使えば、廃棄コストが得られます（リスト3：サンプルファイルの6番目のコードセルに入力されています）。
# リスト3：廃棄コストを計算する：平均のベイズ推定を実行してから、このコードを実行する
m = 15 # 基準値
delta = 0.5 # 許容される誤差
low = m - delta # 基準値の下限
high = m + delta # 基準値の上限
# 予測サンプルを作成する
with model:
ppc = pm.sample_posterior_predictive(trace, random_seed=42, progressbar=False)
# 予測された個々の製品のサイズ
y_pred = ppc.posterior_predictive["obs"].values.flatten()
# 製品のサイズが規格外になる確率
prob = np.mean((y_pred < low) | (y_pred > high))
total_number = 100000 # 全体の個数
cost_per_unit = 1000 # 1個あたりの廃棄コスト
# 全体の個数×規格外の確率×廃棄コストで、全体の廃棄コストを求める
disposal_cost = total_number * prob * cost_per_unit
# 結果の表示
print(f"規格外の確率：{prob:.6%}") # 0.029167%という結果が表示される（実行環境により多少変わる場合がある）
print(f"廃棄コスト：{disposal_cost:.0f}" ) # 29167という結果が表示される
pymcモジュールのsample_posterior_predictiveが予測サンプルを作成するための関数。ppc.posterior_predictive["obs"].valuesにより予測サンプルの値を取り出し、flattenメソッドにより、2チェイン×2000サンプル×30個の予測サンプルを一次元に並べる。後は、それらの予測サンプルが規格外になる確率を求めて、廃棄コストを計算するだけ。
予測サンプルを作成するためには、ベイズ推定によって求めた平均と標準偏差を使うので、その推定を行ったモデルを利用する必要があります。リスト1で、with pm.Model() as model:と記述して、モデルにmodelという名前を付けたので、それを続けて利用します。with model:がそのための記述です。その中で（字下げをした範囲で）、sample_posterior_predictive関数を呼び出して予測サンプルを作成するというわけです。
返り値は、予測サンプルを含むオブジェクトです。上の例ではppcという変数でそのオブジェクトを参照できるようにしています。後は、flattenメソッドを使って、予測サンプルの値を一次元の配列y_predに変換し、14.5未満または、15.5より大きくなる確率を求めて、廃棄コストを計算するだけです。
規格外になる確率を求めるnp.mean((y_pred < low) | (y_pred > high))というコードの意味が少し分かりにくいかもしれませんね。本筋からは外れる話ですが、mean関数の引数として指定した部分について簡単に説明しておきます。以下、話を簡単にするために、ちょっと違った値を使います。例えば、y_predの値が[14.9, 15.1, 14.3, 15.7, 15.2]であったものとしましょう。表1で計算を追いかけてみます。「|」はORを求める演算子です。つまり、どちらかがTrueであれば、Trueとなります。
| ypred | 14.9 | 15.1 | 14.3 | 15.7 | 15.2 |
|---|---|---|---|---|---|
| ypred < low | False | False | True | False | False |
| ypred > high | False | False | False | True | False |
| (y_pred < low) | (y_pred > high) | False | False | True | True | False |
この例であれば、14.3と15.7が規格外なので、その部分だけがTrueとなり、規格内の値である場合はFalseとなる。
結果は、[False, False, True, True, False]となりました。mean関数により、その平均を求めるわけですが、Pythonのブーリアン値は数値として計算される場合、Trueは1、Falseは0として扱われるので、この結果は2/5＝0.4となります。これは、Trueである（つまり、規格外である）確率にほかなりません。リスト3でも、同じ方法を使って、サンプリングされた予測値が規格外である確率を求め、それを基に廃棄コストを計算しているというわけです。
というわけで、規格外となる確率は0.029167％（100,000個中29.167個）で、10万個製作した場合の廃棄コストは29,167円となりました。
コラム arviz.InferenceData形式のデータには何が含まれているか
※本コラムの内容は、pymcのバージョン5.28.5、arvizのバージョン0.22.0での動作に基づきます。
pm.sample関数やpm.sample_posterior_predictive関数の返り値などのarviz.InferenceData形式のデータには数多くの情報が含まれています。その中には幾つかのグループがあり、いずれもxarrayと呼ばれる多次元配列のデータ（xarray.Dataset）となっています。例えば、pm.sampleの返り値には、posterior、sample_stats、observed_dataという3つのグループが含まれています。リスト4のようにグループ名を指定すれば、その内容を一覧表示できます。
trace["posterior"]
pm.sample関数の返り値であるtraceに含まれる"posteror"グループの内容を表示するコード。trace.posteriorというコードでも同じ結果が得られる（実行例は図3）。同様に、ppc["posterior_predictive"]またはppc.posterior_predictiveというコードを入力・実行すれば、予測サンプルの内容を見ることができる。
リスト4のコードを実行すると、図3のような結果が表示されます。サンプリングされた値などのデータは「Data variables:」の下に表示されています。さらに、右端のディスクのアイコンをクリックすると、詳細な値を表示できます。図3は、muの詳細（サンプリングされた平均の一覧）を表示した例です。
図3 traceに含まれるデータのうち、muについての詳細を表示したリスト4の実行結果のうち、muの右に表示されているディスクのアイコンをクリックして詳細を表示したところ。サンプリングされた値が表示されている（もう一度ディスクのアイコンをクリックすると、表示が閉じる）。muやsigmaはリスト1の事前分布の定義で指定した名前。
なお、ArviZ 1.0以降では、ここで説明したInferenceDataは廃止され、xarrayのDataTreeという構造に置き換えられました。そのため最新環境では型名やtrace["posterior"]の戻り値などが本コラムの説明と異なる場合があります。ArviZの変更点などについては、こちらを参照してください。
今回は、正規分布の母数（平均と標準偏差）のベイズ推定についておさらいし、続いて、予測サンプルを作成して廃棄コストを見積もりました。繰り返しをいとわず、ゆっくりとお話ししているので、ベイズ統計の基礎と応用について一歩ずつ理解が深めていけることと思います（分からない部分があっても気にせず、のんびりと進めましょう）。それでも、プログラミングはちょっと苦手だという方に朗報です！
次回は番外編として、グラフィカルなインターフェースで対話的にベイズ推定やベイズ検定ができる神ツールJASPを紹介します。しかも、これまでの古典的な（頻度論による）統計とベイズ統計の両方に対応しています。どうぞ、お楽しみに！
今回使った新出の関数（主なもの）
引数も主なものだけを掲載します。引数に「=値」と書かれたもの（例：=True）は既定値を表します。
plot_trace： トレースプロットを描く関数
- モジュール： arviz
- 形式： plot_trace(trace, var_names=None, chain_prop=None, figsize=None)
- 引数：
- trace： arviz.InferenceData形式のデータ（通常、pymc.sampleの返り値を指定する）
- var_names： グラフを描きたい変数名のリスト
- chain_prop： matplotlib.pyplotのplot関数に指定できるオプションを辞書形式で指定する
- figsize： グラフのサイズ（インチ単位）
※ arviz1.2.0ではchain_propとfigsizeの指定はできない。変更点については、こちらを参照。
sample_posterior_predictive： 予測サンプルを作成する関数
- モジュール： pymc
- 形式： sample_posterior_predictive(trace, random_seed=None, progressbar=True)
- 引数：
- trace： arviz.InferenceData形式のデータ（通常、pymc.sampleの返り値を指定する）
- random_seed： 乱数の初期設定。値を指定しない場合は毎回異なる乱数となる
- progressbar： サンプリングの進行状況を表示する場合はTrueを、表示しない場合はFalseを指定する
- 返り値： 予測サンプルを含むデータ（arviz.InferenceData形式。ただし、PyMC 6以降では、xarray.DataTree形式）。予測サンプルはposterior_predictive属性により得られる
Copyright© Digital Advantage Corp. All Rights Reserved.
アイティメディアからのお知らせ
注目のテーマ
編集部からのお知らせ
ITmediaはアイティメディア株式会社の登録商標です。