コロナ禍で Stay home? 行動変容の数理モデルを考える

0. はじめに

今年は新型コロナウイルスに振り回された一年でした。Go ToトラベルやGo Toイートで経済を活性化したい反面、医療崩壊を避けるためには旅行や外食は避けるべきという意見もあります。私たちはWithコロナ時代の新しい生活様式へと行動変容を迫られているわけですが、人の行動はそんなに簡単には変わらないですよね。行動を変えるには、助成金とかGo Toトラベルなど何らかの補償やインセンティブが必要なのも頷けます。

ところで、人の行動決定がインセンティブやリスクに左右されることを表現する数理モデルとしてマルコフ決定過程（Markov Decision Process, MDP）が知られています。本稿では、コロナ禍における人の行動変容をMDPでモデル化するとともに、行動に伴う**報酬（または損失）**を変えることで行動がどのように変化するかをベルマン方程式を解いて確かめてみようと思います。

なお、本稿の内容は、著者が昨年に会社の仲間とともに書いた書籍「現場で使える！Python深層強化学習入門」とその付属コードを参考にしています。

1. マルコフ決定過程による行動モデル

ここでは、MDPの定義とMDPによるコロナ禍での人間行動のモデル化について解説します。

1-1. マルコフ決定過程

マルコフ決定過程（MDP）では、状態（state）およびその状態もとで選択し得る行動（action）をともに離散変数として定義します。例えば、4つの状態と2つの行動からなるMDPの場合、状態変数は4つの要素からなるカテゴリカル変数、行動は0か1の何れかの値をとるバイナリ変数で記述されます。

MDPでは、ある状態から別の状態へと直接に遷移することはできません。次の状態に遷移する前に、取るべき行動を決定する必要があります。また、次の状態へ遷移する際には、その直前で取った行動に対する**報酬（reward）**を受け取ることができます。状態 $s$ と行動 $a$ のペア $(s, a)$ のもとで次の状態 $s'$ に遷移する確率 $P(s'|s,a)$ は、それに伴う報酬 $R(s,a,s')$ と同様、モデルの定義として予め定められています。

一方、状態 $s$ において行動 $a$ を選ぶ確率 $\pi(a|s)$ は方策（policy）とよばれモデルでは定義されません。方策は、行動の主体である人（あるいはエージェント）の意思決定を担っており、与えられた環境としての MDP のもとで最適なポリシーを見つけることが問題を解くことに相当します。

1-2. 状態遷移図

MDP を定義するには、遷移確率 $P(s'|s,a)$ と報酬 $R(s,a,s')$ を設計する必要があります。それを見通しよく行うには、MDP を図として可視化するのが便利です。

状態遷移図の描き方を説明します。まず、状態変数の各要素を白丸ノード、行動変数の各要素を青丸ノードで記述します。ある状態 $s$ のもとで行動 $a$ を決定した場合、$s$ を表す白丸ノードから $a$ を表す青丸ノードへ向かう有向矢印を描きます。その結果、次の状態 $s'$ に遷移した場合、$a$ を表す青丸ノードから $s'$ を表す白丸ノードへ有向矢印を描きます。このように、状態と行動を表す2種類のノードを有効矢印で結線して得られる2部グラフのことを状態遷移図と呼びます。青丸ノード $a$ から白丸ノード $s'$ へと向かう有向矢印には、対応する遷移確率と報酬 $(p, r)$ を添書きします。

今回の問題設定では、状態変数は Home（家）、Office（職場）、Bar（居酒屋、飲食店）、Hospital（病院）の4つの状態、行動変数は stay（留まる）、move（移動する）の2つの行動からなります。MDP は以下の状態遷移図で表されます。

[図.1：状態遷移図](### 1-1.)

状態遷移図の結線の有無と有向矢印に添書きした遷移確率と報酬 $(p, r)$ が、MDP の設計を決定します。これらの情報を下表にまとめました。各行の $p$ と $r$ は、それぞれ $P(s'|s,a)$ および $R(s,a,s')$ として、同じ行の $s$, $a$, $s'$ と関係づけられています。表の最右列では、3章で説明する4種類の実験に係わる行にフラグを立てています。

状態 s	行動 a	次状態 s'	遷移確率 p	報酬 r	実験
Home	move	Office	0.8	$+$1
Home	move	Bar	0.2	$+$2	B
Home	stay	Home	0.9	$\pm$0	A
Home	stay	Hospital	0.1	$-$2
Office	move	Bar	0.5	$+$2	B
Office	move	Home	0.5	$+$1
Office	stay	Office	0.8	$+$1	D
Office	stay	Hospital	0.2	$-$2
Bar	move	Home	1.0	$+$1
Bar	stay	Bar	0.6	$-$1
Bar	stay	Hospital	0.4	$-$2
Hospital	move	Home	1.0	$\pm$0	C
Hospital	stay	Hospital	1.0	$\pm$0
[表.1：MDP の定義表](### 1-1.)

2. ベルマン方程式

所与の MDP のもとで方策を最適化するには目的関数が必要です。この目的関数は、状態 $s$ および状態と行動のペア $(s,a)$ の価値を表す関数として定義されます。これらの目的関数は、ベルマン方程式（Bellman equation）と呼ばれる再帰方程式を介してエージェントの方策と関係づけられます。ベルマン方程式を満たしつつ目的関数を最大化する方策を求めることが最適化の目的です。

2-1. ベルマン方程式

方策 $\pi$ の目的関数として2種類が定義できます。状態価値関数 $V^\pi(s)$ と行動価値関数 $Q^\pi(s,a)$ です。$V^\pi(s)$ は、方策 $\pi$ に従って行動決定する場合に、それ以後の未来にわたって得られる割り引き報酬和の期待値とし定義されます¹。状態価値関数と行動価値関数は、MDP と方策 $\pi$ を介して次式のように関係づけられます。

\begin{align}
V^\pi (s) &= \sum_a \pi(a|s)\,Q^\pi (s, a) \\
Q^\pi (s, a) &= \sum_{s'} P(s'|s, a) \Bigl[R(s,a,s')
+\gamma V^\pi(s')\Bigr]
\end{align}

この方程式をベルマン方程式と呼びます。第2式は、MDP と状態価値関数から行動価値関数を定義する式として解釈できます。第2式の行動価値関数を第1式に代入すれば、状態価値関数の再帰方程式が得られます。同様に、第1式の状態価値関数を第2式に代入すれば、行動価値関数の再帰方程式が得られます。

2-2. ベルマン方程式を価値反復法で解く

MDP が所与の場合にベルマン方程式を解く方法として、方策反復法と価値反復法があります。ここでは、方策 $\pi$ が決定論的であると仮定して価値反復法を用いて解くことにします。ベルマン方程式は再帰的なので、状態価値関数による行動価値関数の再計算と、行動価値関数を最大化による方策および状態価値関数の更新を収束するまで繰り返します。

\begin{align}
Q_{t+1} (s, a) &= \sum_{s'} P(s'|s, a)
\Bigl[R(s,a,s')+\gamma V_t(s')\Bigr] \\
V_{t+1} (s) &= \max_a Q_{t+1} (s, a) \\
\pi_{t+1}(s) &= \arg \max_a Q_{t+1}(s,a)
\end{align}

以下のコードは、上述の価値反復法を Python で実装したものです。

value_iteration.py

# 各種モジュールのインポート
import numpy as np
import copy

# 遷移確率の設定 p(s'|s, a) --> p[s, s', a]
p = np.zeros((4, 4, 2))
p[0, 1, 0] = 0.8
p[0, 2, 0] = 1 - p[0, 1, 0]
p[0, 0, 1] = 0.9
p[0, 3, 1] = 1 - p[0, 0, 1]
p[1, 0, 0] = 0.5
p[1, 2, 0] = 1 - p[1, 0, 0]
p[1, 1, 1] = 0.8
p[1, 3, 1] = 1 - p[1, 1, 1]
p[2, 0, 0] = 1.0
p[2, 1, 0] = 1 - p[2, 0, 0]
p[2, 2, 1] = 0.6
p[2, 3, 1] = 1 - p[2, 2, 1]
p[3, 0, 0] = 1.0
p[3, 3, 1] = 1.0

# 報酬値の設定 r(s, a, s') --> r[s, s', a]
r = np.zeros((4, 4, 2))
r[0, 1, 0] = 1.0
r[0, 2, 0] = 2.0  # 2.0
r[0, 0, 1] = 0.0  # 0.0
r[0, 3, 1] = -2.0
r[1, 0, 0] = 1.0
r[1, 2, 0] = 2.0  # 2.0
r[1, 1, 1] = 1.0  # 1.0
r[1, 3, 1] = -2.0
r[2, 0, 0] = 1.0
r[2, 1, 0] = 0.0
r[2, 2, 1] = -1.0
r[2, 3, 1] = -2.0
r[3, 0, 0] = 0.0  # 0.0
r[3, 3, 1] = 0.0

# 割引率の設定
gamma = 0.95

# 価値関数の初期化
v = [0, 0, 0, 0]
v_new = copy.copy(v)

# 方策分布の初期化
pi = [0.5, 0.5, 0.5, 0.5]

# 価値反復法の計算
for step in range(1000):

    # 行動価値関数を初期化
    q = np.zeros((4, 2))

    for i in range(4):

        # 行動価値関数を計算
        for j in range(4):
            q[i, 0] += p[i, j, 0] * (r[i, j, 0] + gamma * v[j])
            q[i, 1] += p[i, j, 1] * (r[i, j, 1] + gamma * v[j])

        # 行動価値関数のもとで greedy に方策を改善
        if q[i, 0] > q[i, 1]:
            pi[i] = 1
        elif q[i, 0] == q[i, 1]:
            pi[i] = 0.5
        else:
            pi[i] = 0

    # 改善された方策のもとで価値関数を計算
    v_new = np.max(q, axis=-1)

    # 計算された価値関数 v_new が前ステップの値 v を改善しなければ終了
    if np.sum(v_new - v) <= 0:
        break

    # 価値関数を更新
    v = copy.copy(v_new)

    # 現ステップの価値関数と方策を表示
    print('step:', step, ' policy:', pi, ' value:', v)

3. 実験結果

ここでは、デフォルト設定に対して報酬設定を変更した4通りの場合について、最適方策 $\pi(s)$ がどのように変更されるかを計算により調べました。

デフォルト設定

デフォルト設定では、遷移確率、報酬とも表.1 の通りです。この場合、最適方策は $\pi =(1,1,1,1)$ であり、すべての状態において move を選択する方策となります。家にこもらず、毎日、職場に出勤して定時で仕事を切り上げ、居酒屋にも通うけど長居はせずに家に帰る、というスタイルです。ちなみに病院において move を選ぶのは、病院で完治したら直ぐに家に帰って社会復帰に備えることだと解釈できます。

A : 在宅勤務を支援

これは、表.1 の3行目において報酬を $\pm 0$ から $+2$ に増やすことを意味します。家に居ても仕事をして稼ぐことができるので出社する必要がなくなります。実際、最適方策を計算すると $\pi = (0, 1, 1, 1)$ であり、家では stay、職場や居酒屋、飲食店では move を選ぶというスタイルです。仮に出社したり外食しても長居はしません。

B : リスク回避を徹底

これは、表.1 の2行目と5行目において、報酬を $+2$ から $-2$ に減らすことを意味します。デフォルトでは、居酒屋での飲み会や飲食店での外食は、職場や家庭のコミュニケーションを促進するという意味でプラスに働きますが、コロナ禍では感染リスクを高めるのでマイナスと評価して避けることになります。この場合、最適方策は $\pi = (1, 0, 1, 1)$ であり、職場では stay、それ以外では move を選びます。つまり、職場に出勤したら食事や飲みには行かずに仕事に集中し、外食する場合でも長居はしないというスタイルです。

C : 感染したら治療と社会復帰を支援

これは、表.1 の下から2行目において、報酬を $\pm 0$ から $+6$ に増やすことを意味します。コロナに感染した場合には、治療費を補助したり有効な薬やワクチンを投与するなどして治療と社会復帰を支援します。この場合も、最適方策は $\pi = (1, 0, 1, 1)$ となります。たとえコロナに感染しても治癒して復帰できるという安心感のもと、リスクに注意しながら仕事に励むスタイルと解釈できそうです。

D : 残業を増やす

これは、表.1 の7行目において、報酬を $+1$ から $+3$ に増やすことを意味します。残業が増えて食事や飲みに行く時間がないため、結果的にリスク回避を促します。実際、最適方策を計算すると $\pi = (1, 0, 1, 1)$ となり B、C と同じ結果が得られます。

4. まとめ

実験結果を表.2 にまとめておきました。B、C、D は最適方策は同じですが状態価値の平均値で見ると異なります。リスク回避の B は他に比べて状態価値が最も低くリスクを恐れるあまりコミュニケーション機会が減って幸福度が下がっていると解釈できるかもしれません。残業を増やせば所得も増えるので D の状態価値が高いのは頷けます。次いで在宅勤務を表す A の状態価値が高いのは、家で効率よく働けるからでしょうか。コロナ感染後に手厚くサポートする C は、状態価値が既定値よりも高く、Withコロナ時代の行動変容への示唆を与えていると言えるかもしれません。

実験	最適方策	Home	Office	Bar	Hospital	状態価値
既定値	(1,1,1,1)	move	move	move	move	24.7
A	(0,1,1,1)	stay	move	move	move	28.7
B	(1,0,1,1)	move	stay	move	move	7.5
C	(1,0,1,1)	move	stay	move	move	26.4
D	(1,0,1,1)	move	stay	move	move	30.8
[表.2：実験結果](### 1-1.)

以上、マルコフ決定過程としてコロナ禍での行動変容について計算実験をもとに考察しましたが、ここでの結論は著者による主観的解釈が多分に含まれています。その点、あしからずご容赦ください。

1. 定義式については、冒頭で紹介した書籍の2.2節を参照してください。https://www.shoeisha.co.jp/book/detail/9784798159928 ↩