ランチ最適化問題(バンディットアルゴリズム編)

はじめに

日々の生活を営む上で、ランチをどこで食べるかということは非常に重要な問題です(2回目)。

前回はこの問題を最適停止問題と捉えて探索と活用の最適点探索を行いましたが、最適停止問題では一度見逃した店はもう選択できないという制約がありました。

飲食店で考えると一度利用した店をもう一度利用することが可能なので、 この問いに対してはバンディットアルゴリズムを用いた方が適していると考えられます。 そのため今回はバンディットアルゴリズムを用いた場合の実験を行います。

参考図書

• バンディット問題の理論とアルゴリズム

• ベイズ推論による機械学習入門

バンディットアルゴリズムについて

以下のスライドが分かりやすいです。

• https://www.slideshare.net/greenmidori83/ss-28443892

バンディットアルゴリズムは例えばカジノのスロットマシーンのように、報酬が不明な複数の選択肢がある場合にどのように選択を行えば累積報酬を最大にできるかを解くアルゴリズムで、インターネット広告配信や推薦システムで用いられています。 複数の解き方があり、ε-greedy、UCB、Thompson Samplingあたりが有名です。

ABテストとの違いが疑問として出やすいですが、ABテストはAとBのどちらが優れているかを決めるものであり、 バンディットアルゴリズムはABテストの期間中の報酬を最大になるようにAとBの出現比率を調整するものなので、 目的が異なるものになります。

実験

ベイズ統計の枠組みを利用しているThompson Samplingが面白そうだったので、今回はThompson Samplingを用いて以下の２通り場合の実験を行います。

1. 報酬がベルヌーイ分布に従う場合
2. 報酬が平均・分散未知の正規分布に従う場合

報酬がベルヌーイ分布に従う場合のThompson Sampling

仮定として、1回の試行につき各店舗で得られる利得が0（まずい）か1（うまい）の二値とし、 各店舗のもつ確率は一定であるとします。 この条件における一定の試行回数における合計利得を最大化するケースです。

今回対象店舗を10店舗とし、各店舗の期待値は以下のようになっているものとします。
[0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.45,0.5]
10番目の店舗の期待値が最も大きく0.5であるため、10番目の店舗だけを選択し続けるのが最適解となります。

実施する手順は以下になります。

1. 各店舗の利得予測値をベータ分布からランダムサンプリング
2. 最も利得の大きい店舗を選択、その店舗の利得(0か1)を観測
3. 選択店舗のベータ分布の事前確率を更新
4. 1~3を試行回数分繰り返し

コードは以下のようになります。

def thompson_bernoulli(armparams,T,ps):
    reward=0
    cumulative_reward=[]
    
    for i in range(1,T):
        #ベータ分布よりランダムサンプリング
        thetas = [np.random.beta(a=armparams[index,2], b=armparams[index,3]) for index in range(len(ps))]
        #値が最大のアームを選択
        index=max(enumerate(thetas), key=lambda x: x[1])[0]
        result = binomial(n=1, p=ps[index])
        if result == 1:
            armparams[index,0]+=1
            armparams[index,2]+=1
        else:
            armparams[index,1]+=1
            armparams[index,3]+=1
        
        reward+=result
        cumulative_reward.append(reward)
    return pd.Series(cumulative_reward)/T

#試行回数
T=10000

#アームの数
armnum=10

#それぞれのアームの正解率(当てたいもの)
ps=[0.2,0.2,0.2,0.2,0.4,0.4,0.4,0.4,0.45,0.5]

#初期値
armparams=np.zeros((armnum,4))
armparams[:,2]=1 #alpha
armparams[:,3]=1 #beta

#実行
thompson_bernoulli(armparams,T,ps).plot()

正規化した累積報酬は0.5にかなり近くなりました。

それぞれの店舗のパラメータを確認します。

resultparams=pd.DataFrame(armparams)
resultparams.columns=["win_count","lose_count","alpha","beta"]
resultparams

win_countがその店舗が当たりだった回数でlose_countがハズレだった回数です。 win_count+lose_countがその店舗が選択された回数になります。 確率値が0.2の店舗1~4は30回前後しか選択されず、確率値が0.4~0.45の店舗5~9は100~200回ほど選択されています。 最も選択したい確率値が0.5の店舗10は9000回近くと最も選択できていることがわかります。

各店舗の事後分布を確認します。

#事後分布
from scipy import stats
x = np.linspace(0,1,100)
plt.figure(figsize=(10,6))
num=1
for a,b in zip(armparams[:,2], armparams[:,3]):
    beta_pdf = stats.beta.pdf(x, a, b)
    plt.plot(x,beta_pdf,label=num)
    num+=1
plt.legend()
plt.show()

ベータ分布は更新されるごとに分布が尖っていくので、あまり選択されていない店舗1~4は分布がなだらかで、 最も選択された店舗10ではかなり尖った形になります。 この分布からランダムサンプリングした値を用いるので店舗10を選ぶ確率が非常に大きくなります。

報酬が平均・分散未知の正規分布に従う場合のThompson Sampling

店舗から得られる利得は0(まずい)か1(うまい)の2値ではなく、 その中間の値もあったり同じ店舗でもその日の食材調達状況によって変化したりする場合もあるので実際は幅をもつ連続値であると考える方が自然です。 上記を踏まえ店舗の報酬を正規分布と仮定したケースを考えます。

今回対象店舗を4店舗とし、各店舗の期待値と分散は以下のようになっているものとします。
[[0,0.3],[0,3],[2,3],[3,3]]
期待値が3である4番目の店舗を選択し続けるのが最適解となります。

さっきのベルヌーイ分布の場合は事前分布がパラメータ２つのベータ分布となり更新がそれほど大変ではなかったのですが、 平均・分散未知の正規分布の場合は事前分布がガウス・ガンマ分布となり、4つのパラメータを更新する必要が生じます。 このあたりはベイズ推論による機会学習のp94~97に詳しく書いてあり、それを参考に実装しました。 (期待通りの結果がでたのでおそらく実装はあっているはず…)

from scipy.stats import t, norm, gamma

def thompson_gauss(count,m,beta,a,b,T,ps):
    reward=0
    cumulative_reward=[]
    
    for i in range(1,T):    
        thetas=[]
        for index in range(len(ps)):
            #事前分布からランダムサンプリング
            lambda_ = gamma.rvs(a[index], loc=0, scale=1/b[index])
            sigma_ = np.sqrt(np.reciprocal(lambda_*beta[index]))
            thetas.append(norm.rvs(loc=m[index], scale=sigma_,size=1)[0])
        
        #print(thetas,max(enumerate(thetas), key=lambda x: x[1])[0])
        
        #パラメータの更新
        index=max(enumerate(thetas), key=lambda x: x[1])[0]
        result = norm.rvs(loc=ps[index][0], scale=ps[index][1],size=1)[0]
        
        bf_m=m[index]
        bf_beta=beta[index]
        
        count[index]+=1
        beta[index]=beta[index]+1
        m[index]=(result+bf_beta*m[index])/beta[index]
        a[index]=0.5+a[index]
        b[index]=(result**2+bf_beta*bf_m**2-beta[index]*m[index]**2)/2+b[index]
        
        reward+=result
        cumulative_reward.append(reward)
    
    return pd.Series(cumulative_reward)/T

#試行回数
T=10000

#アームの数
armnum=4

#それぞれのアームの平均と分散
ps=[[0,0.3],[0,3],[2,3],[3,3]]

#パラメータ初期値
count=np.zeros(armnum)
m=np.zeros(armnum)
beta=np.ones(armnum)
a=np.ones(armnum)
b=np.ones(armnum)

thompson_gauss(count,m,beta,a,b,T,ps).plot()

事前分布からランダムサンプリングする部分はガンマ分布と正規分布を分けてサンプリングしていますが、 studentのt分布を用いても導出できるはずです。

正規化した累積報酬は3に近くなり最適な選択が割とできているようです。

パラメータを確認します。

#パラメータの確認
params=pd.DataFrame()
params["count"]=count
params["m"]=m
params["beta"]=beta
params["a"]=a
params["b"]=b
params

countがその店舗が選択された回数になります。 期待値が最も大きい店舗3を選ぶことができていることが確認できます。

最後に事後分布を確認します。 今度はstudentのt分布を使ってみます。

#予測分布
x = np.linspace(-1, 4, 1000)
plt.figure(figsize=(10,6))
for i in range(len(ps)):
    mu_=m[i]
    lambda_=beta[i]*a[i]/((1+beta[i])*b[i])
    neu_=2*a[i]   
    plt.plot(x,t.pdf(x, df=neu_, loc=mu_, scale=lambda_),label=i)
plt.legend()
plt.show()

それぞれの期待値が比較的識別可能な形になっていることがわかります。 試行回数が少ない1や2の店舗の期待値は少し本来のものとはずれることがあるようです。

メモ

• 事前分布の初期値は結果に影響が出そる可能性があるので、状況に応じて変化させたり最初に数回全ての選択肢を試した方がよい場合がありそう。
• 各分布のパラメータの意味をあまりわかっていなかったので教科書上のパラメータとpythonのstatsモジュールのパラメータの対応付けに苦労した。確率分布の勉強はもう一度やった方が良さそう。
• 人間には飽きがあるので、前回選択した店舗にはマイナス値の補正をかけるという方法をとるとより選択する店舗が分散するようになり実際のケースに近づけることができそうかと思った。

今回使用したコードは以下になります。
github.com