In [1]:
import numpy as np
from scipy.optimize import linprog
import quantecon.game_theory as gt
じゃんけん・ゲームを例に計算してみる．
あとあと便利なので，NumPy array としてプレイヤー0の利得行列を定義しておく:
In [2]:
U = np.array(
[[0, -1, 1],
[1, 0, -1],
[-1, 1, 0]]
)
In [3]:
p0 = gt.Player(U)
p1 = gt.Player(-U.T)
プレイヤー1の行列は -U の転置 (.T) であることに注意．
In [4]:
g = gt.NormalFormGame((p0, p1))
In [5]:
print(g)
In [6]:
gt.lemke_howson(g)
Out[6]:
In [7]:
gt.support_enumeration(g)
Out[7]:
主問題: $$ \min u $$ subject to $$ U y - \mathbf{1} u \leq \mathbf{0},\quad \mathbf{1}' y = 1,\quad y \geq \mathbf{0} $$
これを
scipy.optimize.linprog
の形式
$$
\max c' x
$$
subject to
$$
A_{\mathit{ub}} x \leq b_{\mathit{ub}},\quad
A_{\mathit{eq}} x = b_{\mathit{eq}},\quad
l \leq x \leq u
$$
に合わせると，
$$
x =
\begin{pmatrix}
y \\ u
\end{pmatrix},\ %
c =
\begin{pmatrix}
\mathbf{0} \\ 1
\end{pmatrix},\ %
A_{\mathit{ub}} =
\begin{pmatrix}
U & -\mathbf{1}
\end{pmatrix},\ %
b_{\mathit{ub}} = \mathbf{0},\ %
A_{\mathit{eq}} =
\begin{pmatrix}
\mathbf{1}' & 0
\end{pmatrix},\ %
b_{\mathit{eq}} = \begin{pmatrix}1\end{pmatrix}
$$
$$
l =
\begin{pmatrix}
0 & \cdots & 0 & -\infty
\end{pmatrix}',\ %
u =
\begin{pmatrix}
\infty & \cdots & \infty
\end{pmatrix}'
$$
両プレイヤーの戦略の数をそれぞれ $m$, $n$ とする:
In [8]:
m, n = U.shape
各入力を定義する:
In [9]:
c = np.zeros(n+1)
c[-1] = 1
c
Out[9]:
In [10]:
A_ub = np.hstack((U, np.full((m, 1), -1)))
A_ub
Out[10]:
In [11]:
b_ub = np.zeros(m)
b_ub
Out[11]:
In [12]:
A_eq = np.ones((1, n+1))
A_eq[0, -1] = 0
A_eq
Out[12]:
In [13]:
b_eq = np.ones(1)
b_eq
Out[13]:
In [14]:
bounds = [(0, None)] * n + [(None, None)]
bounds
Out[14]:
scipy.optimize.linprog に渡して解かせる:
In [15]:
res_p = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds)
結果:
In [16]:
res_p
Out[16]:
プレイヤー1の均衡戦略:
In [17]:
res_p.x[:-1]
Out[17]:
ゲームの値:
In [18]:
res_p.x[-1]
Out[18]:
scipy.optimize.linprog は双対解を返してくれないようなので，双対問題も定式化して解かせる:
$$
\min -v
$$
subject to
$$
-U' x + \mathbf{1} v \leq \mathbf{0},\quad
\mathbf{1}' x = 1,\quad
x \geq \mathbf{0}
$$
In [19]:
c = np.zeros(m+1)
c[-1] = -1
A_ub = np.hstack((-U.T, np.full((n, 1), 1)))
b_ub = np.zeros(n)
A_eq = np.ones((1, m+1))
A_eq[0, -1] = 0
b_eq = np.ones(1)
bounds = [(0, None)] * m + [(None, None)]
In [20]:
res_d = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds)
res_d
Out[20]:
プレイヤー0の均衡戦略:
In [21]:
res_d.x[:-1]
Out[21]:
ゲームの値:
In [22]:
res_d.x[-1]
Out[22]:
微妙は誤差はある:
In [23]:
res_p.x[-1] - res_d.x[-1]
Out[23]:
In [24]:
def solve_zerosum_lemke_howson(U):
g = gt.NormalFormGame((gt.Player(U), gt.Player(-U.T)))
NE = gt.lemke_howson(g)
return NE
In [25]:
def solve_zerosum_linprog(U, method='revised simplex'):
U = np.asarray(U)
m, n = U.shape
# Primal problem
c = np.zeros(n+1)
c[-1] = 1
A_ub = np.hstack((U, np.full((m, 1), -1)))
b_ub = np.zeros(m)
A_eq = np.ones((1, n+1))
A_eq[0, -1] = 0
b_eq = np.ones(1)
bounds = [(0, None)] * n + [(None, None)]
res_p = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method=method)
# Dual problem
c = np.zeros(m+1)
c[-1] = -1
A_ub = np.hstack((-U.T, np.full((n, 1), 1)))
b_ub = np.zeros(n)
A_eq = np.ones((1, m+1))
A_eq[0, -1] = 0
b_eq = np.ones(1)
bounds = [(0, None)] * m + [(None, None)]
res_d = linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method=method)
NE = (res_d.x[:-1], res_p.x[:-1])
return NE
In [26]:
solve_zerosum_lemke_howson(U)
Out[26]:
In [27]:
solve_zerosum_linprog(U)
Out[27]:
In [28]:
m, n = 4, 3
U = np.random.randn(m, n)
U
Out[28]:
In [29]:
solve_zerosum_lemke_howson(U)
Out[29]:
In [30]:
solve_zerosum_linprog(U)
Out[30]:
In [ ]: