即売会釣り銭問題

　文学フリマでヘルツシュタルク氏へのインタビュー本を頒布いたしました。ぴっちぶれんどのブースに来ていただいた方々、どうもありがとうございました。
　さて、今回の文学フリマでは、ある統計を取っていました。それは、みなさんがどの券種のお金で購入しようとするかです。即売会のたびに釣り銭、特に大量に必要な100円玉をどの程度用意しておけば良いのか迷うので、その消費量を調査し必要量を数理的に決定して、次回に生かそうと考えたわけです。

　次の即売会では700円の本のみを100冊販売すると仮定します。そのために、お客さんが700円の本を買う場合にどのような券種を差し出すかを調査したわけですが、今回は300円で併売した「発笑探検隊」さんのクルタマニュアル本と一緒に購入していただいた方が多く、またメモの取り忘れなどがあり、700円の本をのみを購入する場合の動向データは19人についてのみ得られました。結果は
・700円ちょうど（100円玉2枚+500円玉）を出した人は7人
・1,000円札または10,000円札の人は10人
・1,200円（100円玉2枚+1000円札）を出して500円のお釣りを受け取った人が2人
でした。それぞれの場合で手元の100円玉は、+2枚、-3枚、±0枚増加します。その確率を順にα、β、1-α-βとすると、釣り銭問題は次のように言い換えられます。

　ある人が、原点から出発して1秒ごとに確率αで前に2歩進み、確率βで後ろに3歩下がり、その他の場合は動かないとする。時刻0秒から100秒までの位置101点について、もっとも後ろ側に存在する位置の確率的分布はどのようになるか。

　これはランダムウォーク問題になります。歩数が対称ならばカタラン数云々で解析的に解けると思うのですが、そう簡単に行きそうにありません。そこで、乱数を振ってモンテカルロ・シミュレーションを実行してみます。Mathematicaを使えば簡単です。
　とりあえずはデータを元にα=7/19、β=10/19と置いて10,000ケース（100万人分）を試行してみました。

alpha=7/19;beta=10/19;
dist1=Table[Min[0,Accumulate[
RandomVariate[
ProbabilityDistribution[
Piecewise[{
{alpha,x==2},
{beta,x==-3},
{1-alpha-beta,x==0},
{0,True}}],
{x,-3,2,1}],100]
]],{10000}]//EmpiricalDistribution;
DiscretePlot[CDF[dist1,x],{x,-250,10}]
Mean[dist1]//N
Quantile[dist1,0.05]

　これが、釣り銭が最大で何枚減るかの累積確率分布です。グラフのX軸の数字に対して、Y軸の値はその枚数以上釣り銭が減る確率がどの程度かを表しています。例えば釣り銭が足りなくなる危険率を0.05にしたいならば、あらかじめ124枚を持って行けば良いことになります。ちなみに平均値は86.7枚でした。
　ところが今回のデータサンプル数が19と少ないことから、α、βの推定値自体の信頼性に問題があります。そこで、ベイズ推定をやってみたいと思います。事前分布を一様分布だと仮定すると、α、βの事後確率は次元2、形状母数(7+1,10+1,2+1)のディリクレ分布に従います。グラフにするとこんな感じです。

Plot3D[PDF[DirichletDistribution[{7+1,10+1,2+1}],{x,y}],
{x,0,1},{y,0,1},PlotRange->All]

　これを元に、先ほどのプログラムを修正してみました。α、βの値を乱数で決定し、さらに100人分の試行で乱数を振るわけです。結果はこうなりました。

dist2=Map[Min[0,Accumulate[
RandomVariate[
ProbabilityDistribution[
Piecewise[{
{#[[1]],x==2},
{#[[2]],x==-3},
{1-#[[1]]-#[[2]],x==0},
{0,True}}],
{x,-3,2,1}],100]
]]&,RandomVariate[DirichletDistribution[{7+1,10+1,2+1}]
,10000]]//EmpiricalDistribution;
DiscretePlot[CDF[dist2,x],{x,-250,10}]
Mean[dist2]//N
Quantile[dist2,0.05]

　この場合、平均値は82.9、危険率5%で164枚となりました。つまり100円玉は164枚用意しておけば95%の確率で釣り銭不足は発生しないということになります。次の即売会では、この計算をもとにして準備しようと思います。