チョコボール統計

チョコボールの秘密を統計解析で明らかにしていく。おもちゃのカンヅメ欲しい。

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

【概要】

  • 正味重量の分布を眺めてみたら、ちょっと気になることが出てきたのでデータを分析してみた
  • 賞味期限が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年以降の商品は、平均値が下がっていることから、単純にしきい値を変えただけではないことは明らかです。また、分散が低くなっていることからも、製造ラインに何らかの変化が起きたことは明らかであろうと考えています。

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

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

  • 変化位置を目視で決めている → 次回は変化点の推定まで行ってみようと思います
  • 計測器の劣化である可能性もあるが未検証
  • 正規分布のパラメータ推論の問題にするのはデータ傾向からちょっと乱暴

次回以降、これらの点について検証を進めていこうと思います。

【トップに戻る】

実装コード

【トップに戻る】

参考文献

  1. 須山敦志, 機械学習スタートアップシリーズ ベイズ推論による機械学習, 講談社, 2019

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

  • 作者:須山 敦志
  • 出版社/メーカー: 講談社
  • 発売日: 2017/10/21
  • メディア: 単行本(ソフトカバー)

  1. C.M.ビショップ(著), パターン認識機械学習 下 (Pattern Recognition and Machine Learning), Springer, 2007

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

  1. Pythonによるベイズ統計モデリング pymc3を使った入門書。 理論面も多少書かれており、こちらで手を動かしながら学習するのが良いと思う。

【トップに戻る】

広告

Amazonの欲しいものリスト作りました。 チョコボールのカンパ募集中です。
チョコボールをカンパする

森永製菓  チョコボール<ピーナッツ>  28g×20箱

森永製菓 チョコボール<ピーナッツ> 28g×20箱

  • 発売日: 2016/03/01
  • メディア: 食品&飲料

*1:ピーナツ味のデータ量が多いので、これだけを対象にします

*2:賞味期限の記録がないなど

*3:2019年12月以前の最小値を下回るデータが増えたので、この傾向は目でみて明らかな気はします

*4:精度(分散)は0より大きい実数ですので、そのようなパラメータを表現するためガンマ分布を採用しました

*5:計20,000サンプル