ボージョレ・ヌヴォーのコピペをMCMCで分析してみた

#はじめに
Qiitaには初投稿となります。よろしくお願いします。
毎年この時期になると、ボージョレ・ヌヴォーがニュースに取り上げられ、その度に例のコピペが話題になります。あれを元に、何年産がおいしい年なのかというランキングを作ろうとした人は、多いのではないかと思います。しかし、難しいのですね。1位は分かるにしても、それ以下の順位を出そうとすると情報の不足に突き当たります。

一応、グラフ化している人もいるようですが、もう少し定量的にやりたい。しかし全ての情報が与えられているわけではない。となれば、こういうときはベイズ推定でしょう、というわけでやってみました。

#分析の準備
##使用したコピペ
まずは、参考にしたコピペです。

• 95年「ここ数年で一番出来が良い」

• 96年「10年に1度の逸品」
• 97年「1976年以来の品質」
• 98年「10年に1度の当たり年」
• 99年「品質は昨年より良い」
• 00年「出来は上々で申し分の無い仕上がり」
• 01年「ここ10年で最高」
• 02年「過去10年で最高と言われた01年を上回る出来栄え」「1995年以来の出来」
• 03年「100年に1度の出来」「近年にない良い出来」
• 04年「香りが強く中々の出来栄え」
• 05年「ここ数年で最高」
• 06年「昨年同様良い出来栄え」
• 07年「柔らかく果実味が豊かで上質な味わい」
• 08年「豊かな果実味と程よい酸味が調和した味」
• 09年「50年に1度の出来栄え」
• 10年「2009年と同等の出来」「今年は天候が良かった為、昨年並みの仕上がり。爽やかでバランスが良い」
• 11年「2009年より果実味に富んだリッチなワイン」「出来が良く、豊満で絹のように滑らかな味わい」
• 12年「ボジョレー史上最悪の不作」「糖度と酸度のバランスが良く、フルーティーな味わい」
• 13年「みずみずしさが感じられる素晴らしい品質」
• 14年「太陽に恵まれ、グラスに注ぐとラズベリーのような香りがあふれる、果実味豊かな味わい」
• 15年「過去にグレートヴィンテージと言われた2009年を思い起こさせます」
• 16年「採れたて感あふれる新酒」
• 17年「今世紀最高と言われた2015年を彷彿とさせる芳醇さ」←NEW!

##前提事項
分析では、次のような仮定・手法を利用しました。

・ボージョレ・ヌヴォーの味の年ごとの統計は、標準正規分布に従うものとする。これを事前分布とする。
・記述が複数ある年は、具体的な文章の方を選択する。
・「○○年に一度の出来」などと書かれた年は、標準正規分布の累積分布関数を使って定量化させる。（MCMCで振らない。）
・「過去○○年で最高」「○○年より上」などと書かれた年は、その条件に合わない場合に尤度を大きく下げる罰則を設ける。
・「○○年と同等」の年は、○○年の結果をそのまま利用する。
・その他「中々の出来栄え」「申し分の無い仕上がり」などの年は定量化せず、無情報として扱う。（結果的にはほぼ事前分布の通りとなる。）

「10年に1度」と「ここ10年で最高」は明確に区別していますが、12年「史上最悪の不作」だけ少し手を抜いてしまいました。後者または両方のハイブリッドでやるべき所を、前者だけで計算しています。結果にはほとんど影響しないと思いますが。

これらの条件の下で、MCMCで分析してみます。
##年と年の拘束関係
年同士に拘束関係があるところに線を引くと、次のようなグラフになります。完全に独立しているのは16年のみですね。

#コーディング
筆者は、最近はMathematicaしか使えない体になってしまったので、MCMCもMathematicaです。Qiitaユーザーで使っている人は少ないのでしょうね。ごめんなさい。

Mathematicaで使えるMCMCのツールと言えば、https://github.com/joshburkart/mathematica-mcmcにあるものぐらいです。筆者は、これを少しだけ改造して使っています。
##「○○年に一度の出来」「○○年と同様」

f$1996 = InverseCDF[NormalDistribution[], .9]; (* 96年「10年に1度の逸品」*)
f$1998 = InverseCDF[NormalDistribution[], .9]; (* 98年「10年に1度の当たり年」*)
f$2003 = InverseCDF[NormalDistribution[], .99]; (* 03年「100年に1度の出来」*)
f$2006 = f[2005]; (* 06年「昨年同様良い出来栄え」*)
f$2009 = InverseCDF[NormalDistribution[], .98]; (* 09年「50年に1度の出来栄え」*)
f$2010 = f[2009]; (* 10年「2009年と同等の出来」*)
f$2017 = f[2015]; (* 17年「今世紀最高と言われた2015年を彷彿とさせる芳醇さ」*)

変数名がf[○○○○]だったりf$○○○○だったりとぶれていますが、これは本当は前者を使いたいのですが、MCMCのツールが前者を変数として認識してくれないので、とりあえず前者で書いて後で変換する、変換できないところは最初から後者で書くという複雑な事情によります。

##尤度の式

frule = f[x_] :> ToExpression[ToString[f] <> "$" <> ToString[x]];
expr = {
	(* 95年「ここ数年で一番出来が良い」*)
	-1000 Boole[Max[Table[f[n], {n, 1991, 1994}]] > f[1995]], 
	(* 97年「1976年以来の品質」(正規分布15年分と91-96年より上)*)
	Log[PDF[OrderDistribution[{NormalDistribution[], 16}, 16],f[1997]]],
	-1000 Boole[Max[Table[f[n], {n, 1991, 1996}]] > f[1997]],
	(* 99年「品質は昨年より良い」*)
	-1000 Boole[f[1999] <= f[1998]],
	(* 01年「ここ10年で最高」*)
	 -1000 Boole[Max[Table[f[n], {n, 1991, 2000}]] > f[2001]],
	(* 02年「過去10年で最高と言われた01年を上回る出来栄え」*)
	-1000 Boole[f[2002] <= f[2001]],
	(* 05年「ここ数年で最高」*)
	-1000 Boole[Max[Table[f[n], {n, 2001, 2004}]] > f[2005]],
	(* 11年「2009年より果実味に富んだリッチなワイン」*)
	-1000 Boole[f[2011] <= f[2009]],
	(* 12年「ボジョレー史上最悪の不作」*)
	-1000 Boole[InverseCDF[NormalDistribution[], 1/(2012 - 1968)] < f[2012]],
	(* 17年「今世紀最高と言われた2015年を彷彿とさせる芳醇さ」*)
	-1000 Boole[Max[Table[f[n], {n, 2001, 2014}]] > f[2015]]
    } // Total;
(* 高速化のため、正規分布の密度関数の対数を自前で用意した。*)
logNormalDistributionPDF[mu_, sigma_, x_] := -((-mu + x)^2/(2 sigma^2)) - 1/2 Log[2 棠・棧] - Log[sigma];
(* 事前分布を加える。*)
expr = expr + Total[Map[logNormalDistributionPDF[0, 1, #] &, vars]] //. frule;

2012年「ボジョレー史上最悪の不作」は、輸出が本格的に始まった1968年を起点として、それ以降で最悪であると見なして正規分布の累積分布関数から値を求め、それ以下の味であるという式を作っています。また分析の起点を1991年として、95年「ここ数年で一番出来が良い」は91～94年をいずれも上回るという式になっています。（なので、本当は91年から計算しています。）

##MCMCを回す

<< mcmc.m
vars = f /@ 
  Flatten[{Range[1991, 1995], 1997, Range[1999, 2002], 2004, 2005, 
    2007, 2008, 2011, 2012, 2013, 2014, 2015, 2016}]
mcmc = MCMC[expr, Map[{#, 0, .1, Reals} &, vars /. frule], 80000000]

8000万回もやれば十分でしょう。

#実行
いよいよ実行です。途中経過はこんな感じで、上下に薄い色で示しているのは90%と95%の範囲です。採択率34.3%ということで、中々良い感じだと思います。ここ以降、味を表す度合いを、一般的に馴染みがある偏差値に変換して表示しています。

それにしても、MCMCを回すのが楽しくて仕方ありません。どこかに１日中MCMCを回してお金をもらえる仕事はないですかね。

#実行結果
ジャーン！

箱ひげ図で表現してみました。黄色いボックスは50%の確率で味の偏差値がこの範囲に入ることを、上下の髭はこの範囲外の値はほぼありえないことを示します。

##ランキング
事後分布の平均値で見た、ランキングは次の通りです。
1, 2位：2015年, 2017年(78.9)
3, 4位：2005年, 2006年(75.4)
5位：2003年(73.3)
6位：2011年(73.2)
7, 8位：2008年, 2009年(70.5)
9位：2002年(70.4)
10位：2001年(67.7)
今年は15年とタイであり(そういう仮定でした)、95年以降で最高の年ということになります。まあ、今世紀最高というのが入っているとそうなりますよね。5位と6位は微妙で、2003年＞2011年となる確率は58.7%です。

#さいごに
尤度関数には、まだまだ改善の余地があるかもしれません。それによって微妙な年の順位が入れ替わるかもしれないので、今後も研究して行きたいと思っています。無情報とした年も、何らかの方法で定量化したいところです。AHPとかに頼るしかないですかね。