チョコボールの内容量を調べてみる

【概要】

正味重量の分布を眺めてみたら、ちょっと気になることが出てきたのでデータを分析してみた
賞味期限が2020年の商品から重量傾向が変わっているようだ（重量が低下している）
製造装置の精度が向上して余分なマージンを減らしても良くなったとか？

【目次】

はじめに
データの確認
内容量の推定
おわりに
実装コード
参考文献
広告

はじめに

前回、計測記事のナンバリングが200回を超えたことを記念して、これまでに計測してきたチョコボールについてざっくりと集計してみました。

chocolate-ball.hatenablog.com

この記事の最後に、一箱の重量の月毎の遷移を掲載しました。この図、とても興味深いですね。2020年以降が賞味期限の商品の重量が減っているように見えます。
そこで今回は、この重量の低下傾向が統計的に確かなのかを調査してみたいと思います。

【トップに戻る】

データの確認

まずはデータを可視化してみましょう。対象データは2019年12月29日時点までに計測したピーナツ味です*1。

計測数は全体で481個です。

項目	数
計測データ数	481
銀のエンゼル出現数	19
金のエンゼル出現数	1

このうち、計測ミス*2を除外した477個が分析の対象です。

賞味期限月毎の計測数は以下のとおりです。なお、チョコボールは製造年月日の記載が無く、賞味期限月の記載しかないため、賞味期限月の単位で集計をしていきます。

f:id:hippy-hikky:20200102222354p:plain — 賞味期限月毎の計測データ数．横軸は賞味期限月，縦軸は計測データ数．

次に、賞味期限月毎の重量を箱ひげ図でプロットしたものが以下の図です。

f:id:hippy-hikky:20200102222549p:plain — 正味重量の分布．横軸は賞味期限月，縦軸は正味の重量(g)．

この図から2020年2月を境に傾向が変化したように見えるので、2020年2月より前と以降の2つのグループにデータを分けてヒストグラムを描いてみます。

f:id:hippy-hikky:20200102223511p:plain — 正味重量のヒストグラム．水色の領域は全データのヒストグラム．青線は2020年1月より前の期間のデータ．緑線は2020年1月以降の期間のデータ．

重量が低下していそうなことがはっきりとわかりますね。

なお、これらの図に描かれている赤線はパッケージに記載されている内容量（28g）のラインです。これまでの計測では、内容量を下回る商品は確認されていません。

【トップに戻る】

内容量の推定

前章でのデータ可視化結果から、2020年2月以降が賞味期限の商品は重量の傾向が変わっているように見えます。しかし、各月のデータ数は10個程度なので、データの偏りである可能性も否定はできません*3。そこで、正味重量を確率モデルで表現し、パラメータの推論をして、この傾向が明らかなのかを確認してみます。

モデル

重量 $Y=\left\{y_1, y_2, \cdots, y_N\right\}$ はガウス分布 $\mathcal{N}(y_n | \mu, \lambda)$ に従うと仮定することにします。この分布の平均パラメータ $\mu$ もガウス分布 $\mathcal{N}(\mu | \mu_m, \lambda_m)$ に従うと仮定します。 $\lambda$ は精度パラメータで、分散の逆数となります。この精度パラメータはガンマ分布 $\mathrm{Gam}(\lambda | a, b)$ に従うと仮定します*4。

ここでは単純なモデルとして、先に確認したデータに基づき、2020年1月を境にデータを2つのグループに分割し、それぞれパラメータ( $\mu, \lambda$ )を推論することにします。変化点の推論については後で別にブログを書こうと思います。

このようにするとモデル（同時分布）は次のようになります。

${ \displaystyle \begin{eqnarray} p(Y, \mu_k, \lambda_k) = p(Y | \mu_k, \lambda_k)p(\mu_k)p(\lambda_k) = \prod^N_n \mathcal{N}(y_n | \mu_k, \lambda_k) \mathcal(\mu_k | \mu_{mk}, \lambda_{mk}) \mathrm{Gam}(\lambda_k | a_k, b_k) \end{eqnarray} }$

ここで、kはデータのグループを示すインデックスで、今回は $k=\{0,1\}$ です。

パラメータ推論

推論対象のパラメータは、平均パラメータ $\mu$ と精度パラメータ $\lambda$ です。

ガウス分布の平均と精度パラメータの推論は、共役事前分布であるガウス・ガンマ分布を利用することで解析的に求めることができます。しかしここでは、MCMC アルゴリズムを利用した近似解を求めます。解析解については、「ベイズ推論による機械学習入門」などを参照ください。また、MCMC アルゴリズムについては、「PRML(下)」などを参照ください。MCMC アルゴリズムを利用した推論の実装については、今回はPyMC3というフレームワークを利用します。PyMC3については、公式ドキュメントや「Pythonによるベイズ統計モデリング」などを参照ください。

今回実装したコード全体は「実装コード」にnotebookを貼っているので参照ください。重要なところは以下の部分です。これで上述のモデルを定義しています。 $\mu$ の事前分布はパッケージ記載の内容量である28gを平均値としたガウス分布です。精度パラメータ $\lambda$ の事前分布は $a=1, b=5$ のガンマ分布です。

spec = 28.
tau_m = 0.1
a, b = 1., 5.

with pm.Model() as model:
    mus = pm.Normal('mu', mu=spec, tau=tau_m, shape=2)
    taus = pm.Gamma('tau', alpha=a, beta=b, shape=2)
    weights = pm.Normal('weights', mu=mus[idx], tau=taus[idx], observed=data.net_weight.values)

このようにモデルを定義すれば、あとはMCMC アルゴリズムを利用して事後分布からサンプリングをします。PyMC3では以下のように書きます。

with model:
    trace = pm.sample(5000, chains=4)

サンプルは5,000点で、このサンプルを4セット繰り返します*5。

結果

推論結果は次のようになりました。

f:id:hippy-hikky:20200103212358p:plain — サンプル点の傾向．上段が平均パラメータ、下段が精度パラメータ．左がサンプルのヒストグラム(KDE)，右がサンプル系列のプロット．

右のヒストグラムを見ると4セットのサンプルがほぼ重なっています。右のサンプル系列のプロットを見ると、トレンドのないランダムなサンプルのように見えます。ということで、求めたい事後分布からサンプルを得られていると見て良さそうです。

（上図とほとんど変わらないですが）平均パラメータ $\mu_{0}, \mu_{1}$ の事後分布の推論結果は以下のとおりです。

f:id:hippy-hikky:20200103213205p:plain — 平均パラメータの推論結果．青の分布が2020年1月より前のデータ群，オレンジの分布が2020年1月以降のデータ群．赤線が仕様上の内容量．

この結果から、明らかに平均値が下がったことがわかります。どのくらい下がったのかについては以下の分布となります。

f:id:hippy-hikky:20200103213446p:plain — 2つのデータ群の平均値の差の分布．横軸は2群の差の重量(g)．

ということで、2020年1月以降の商品はそれ以前の商品と比較して0.6g~0.75g程度軽くなっているという推論結果となりました。

【トップに戻る】

おわりに

以上の推論結果より、2020年以降が賞味期限の商品の平均重量が下がったことはデータとしてほぼ確信できるということがわかりました。パッケージ記載の内容量に近い商品が多くなったということです（内容量を下回る商品は確認されていません）。

公表されている内容量に実際の量が近づいている、また、精度も向上(分散が低下)していることから、製造装置の精度が向上し、余分なマージンを取らなくても良くなったということが考えられます。
2019年までの内容量傾向を見ると歪んだ正規分布状です。これは、29gくらいの量でしきい値を設定し、下回る商品は出荷しないとしているように見えます。2020年以降の商品は、平均値が下がっていることから、単純にしきい値を変えただけではないことは明らかです。また、分散が低くなっていることからも、製造ラインに何らかの変化が起きたことは明らかであろうと考えています。

何が起きているのかの真相は不明ではありますが、データからこのような事情を推察するのは楽しいものですね。

しかしながら、検証が甘い点がいくつかあります。