チョコボールの大きさの変化点を捉える

概要
はじめに
データ
アプローチ
モデル
- チョコボールの個数の変化点予測
- 1個あたり平均重量の変化点予測
変化点の予測結果
- チョコボールの個数の変化点予測結果
- 1個あたり平均重量の変化点予測結果
まとめ
付録
参考文献
広告

概要

ここ最近、チョコボールの粒が大きくなったような気がする
ベイズ推定により、変化点の検出を試行
どうやら、賞味期限が2018年6月から2018年7月に変わるタイミング（計測時期的には2017年12月中旬）のあたりで変化したらしいと予測できた
個数が1~2個減り、チョコボールの粒の平均重量が約0.2g（11%）重くなっているらしい

はじめに

去年の11月から半年弱、チョコボールの開封を不定期にやっているのですが、最近のチョコボールは粒が揃っていて、大きい気がします。
例えば、このように。

f:id:hippy-hikky:20180415220437j:plain:w200 — 2017年11月26日計測

f:id:hippy-hikky:20180415220453j:plain:w300 — 2017年4月14日

普通はこのように写真を並べて目視で評価するかもしれませんが、目視では確認できる量に限界があります。また、主観に左右されるので信頼できないですよね。（上の写真じゃ撮影位置が違うし全然変化がわからない。。。）
そこで、差があるのか、いつから差が出たのかを統計的に評価してみます。

この問題は、時系列データの変化点検出に応用できます。

【トップに戻る】

データ

今回使うデータは、ピーナツ味のデータのみを使います。他の味は機械的に作るのでバラつきがほとんど無いでしょうし、開封している数もピーナツ味が圧倒的に多く、他は少ないので。

ピーナツ味の全体の集計は下表の通りです。

項目	数
データ数	１２９箱開封
平均内包個数	16.5個
１個あたりの平均重量の平均	1.79g

データの散布図を以下に示します。

f:id:hippy-hikky:20180415232110p:plain — 左：全体重量の散布図。中央：個数の散布図。右：１個あたり平均重量の散布図。横軸は全て開封数。

データは本来であれば、製造年月日でデータを並べるのが正しいのですが、製造日がわからないので、賞味期限でデータを並べています。（ただ、賞味期限は1月単位でしか表示がないので、その間の並びは不定です）

上図を見ると、全体重量（左図）は変化していないのに、チョコボールの数（中央図）は６０箱目くらいで少なくなっているように見えます。この傾向は１個あたり平均重量（右図）が顕著で、一箱あたりのチョコボール数が減っているので、１粒あたりの平均重量が重くなっている（＝大きくなっている）ことが見えてきます。

では、実際に統計的に変化が起きていると言えるのか検証してみます。

【トップに戻る】

アプローチ

変化点の検出には、毎度おなじみのベイズ推定を用います。

変化が１回だと仮定し、変化前後での母分布のパラメータを推定します。また、変化の起こった日（何個目のチョコボールから変化したのか）の予測も合わせて行います。ちなみに、本当に変化が起きていない場合には、予測される変化点前後のパラメータには差異が出て来ないはずです。

実はこの問題、参考文献に載っている例とほぼほぼ同じです。詳しくは書籍を読んでいただけたらと。

【トップに戻る】

モデル

変化点の有無はチョコボールの個数についての分析と、１個あたり平均重量の分析の２種類の分析を行うことで評価してみます。

チョコボールの個数の変化点予測

はじめに１箱あたりに入っているチョコボールの数が変化しているのかを検証するためのモデルを立てます。

i箱目のチョコボールの個数 $N_i$ はパラメータ $\lambda$ のポアソン分布に従うと仮定します。

${ \displaystyle N_i \sim Poi(\lambda) }$

変化点 $\tau$ でポアソン分布のパラメータが $\lambda_1$ から $\lambda_2$ に変化したとします。ちなみに、ポアソン分布の期待値はパラメータ $\lambda$ になりますので、 $\tau$ 箱目から、チョコボール数の平均が $\lambda_1$ から $\lambda_2$ に変化したということをモデル化したものになります。

$\lambda$ の事前分布には、指数分布を用います。

${ \displaystyle \lambda_i \sim Exp(\alpha) }$

ここで、 $\alpha$ は事前分布のパラメータであり、あまり過敏にならなくても良いものらしいです。データ数が少なすぎなければ、事前分布がずれていてもリーズナブルな事後分布に収束してくれるはずです。

$\tau$ も予測対象なので、事前分布を決める必要があります。しかし、いつ変化したかという情報は何も持っていないので、一様分布を事前分布とします。

以上をまとめてpymcのコードにしたものは、付録を参照ください。

1個あたり平均重量の変化点予測

次に、一粒あたりの重量が変化しているのかを検証するためのモデルを立てます。これは、１粒あたりの重量を予測する問題を応用します（過去の記事参照）。

i箱目の１個あたり平均重量 $w_i$ は $N(\mu, \sigma^{2})$ の正規分布に従うと仮定します。

${ \displaystyle ｗ_i \sim N(\mu, \sigma^{2}) }$

上記個数の変化点予測の場合と同様に、変化点 $\tau$ で正規分布のパラメータが $(\mu_1, \sigma_1^{2})$ から $(\mu_2, \sigma^{2}_{2})$ に変化したと仮定します。

今回は、 $\mu$ の事前分布には中心が観測データの平均、分散を100とした正規分布を使ってみます。 $\sigma^{2}$ の事前分布には一様分布を使うことにします。以前の記事と同じ設定です。

変化点 $\tau$ は上記と同じで一様分布を用います。

こちらのコードについても付録を参照ください。

【トップに戻る】

変化点の予測結果

変化点を予測した結果を見てみます。

チョコボールの個数の変化点予測結果

個数の変化点を予測した結果が以下の図です。

f:id:hippy-hikky:20180418000422p:plain — 上：lamba_1の事後分布, 中央：lambda_2の事後分布, 下：tauの事後分布

$\lambda_1$ と $\lambda_2$ の差異があるかを確認したところ、 95%以上の確信度で差異があると予測しています。変化点の前後で、約17個から15~16個に減っているという予測結果になりました。総重量が変化していないので、個数が減ったということは、粒が大きくなったということです。

変化点については、 60箱目くらいから80箱目くらいで変化していると予測しています。先に示した散布図からも同様のことが読み取れますね。

しかし、変化点の分布は広く、変化点を正確には予測できないようです。実際、上記散布図の中央の個数の散布図を見てみると、目で見てもあまりはっきりと変化点を読み取ることが難しいですね。

1個あたり平均重量の変化点予測結果

平均重量の変化点を予測した結果が以下の図です。分散の事後分布については、ここでは掲示していません。

f:id:hippy-hikky:20180418000533p:plain — 上：mu_1の事後分布, 中央：mu_2の事後分布, 下：tauの事後分布

結果は、個数の変化点予測と同じように、60箱目から80箱目付近でチョコボールの傾向に差が出ているということを示しています。同じものを見方の異なるデータで表現したに過ぎないので、結果が同じになるのは当然ですね。
しかし、下段の $\tau$ の予測分布を比較してみると、明らかにこの場合の $\tau$ の分布が特定の値（60箱目、74箱目）に集中していることがわかります。これは、60箱目か74箱目で変化しているということに強い確信を持っていることを示しています。 2山できているので、変化点が一つであるとの仮定が間違っていた可能性が高いです。

変化点前後の事後分布から、 1個あたり約0.2g重くなっているということがわかりました。全体の1個あたり重量の平均が約1.8gなので、11%も重くなった（≒粒が大きくなった）ということが言えそうです。

【トップに戻る】