2020-06-03

【Python】path の確認

Python

ライブラリや他のディレクトリのモジュールをimport するときにどこを見に行っているのかを確認する

import pprint
import sys

pprint.pprint(sys.path)

これで見ることができます．
sys, pprint ともに標準ライブラリです．
sys は基本的な設定や実行環境の変数がまとめられたライブラリで，sys.path で格納されたパスを見ることができます．
pprint はデータを見やすい形にして出力してくれるライブラリで，pprint.pprint でリストの要素ごとに改行して表示してくれます．

2020-05-23

運搬経路問題（配送最適化問題，Vehicle Routing Problem) をPuLPで解く

最適化

Qiitaの記事

qiita.com

2020-05-22

【強化学習】モデルフリー（モンテカルロ法，Q学習，SARSA，Actor-Critic）

強化学習

強化学習について軽くおさらいします．

強化学習は，長期的な報酬が最大化するような方策を学習することが目的でした．

学習対象である良い方策を求めるためには，長期的な報酬（価値）を求める必要がありました．

環境（モデル：遷移関数，報酬関数）がわかっていればモデルベースの手法（価値反復法，方策反復法）で価値を推定し方策を求めることができました．

しかし，環境が未知の場合は，エージェントと環境の相互作用に得られたデータをもとに価値，方策推定をします．モデルを陽に使わずに方策を求める方法を「モデルフリー」な手法と呼びます．

モデルフリーな手法では，推定した価値と経験して得られた結果の差を小さくするように価値を更新し推定していきます．

モンテカルロ法

更新を1エピソードで行う手法をモンテカルロ法といいます．更新式は次式です．

$V^{\pi}(s_t) \leftarrow V^{\pi}(s_t) + \alpha(G_t - V^{\pi}(s_t))$

$G_t$ は割引報酬和です．実際に得られた即時報酬を割り引いて足し合わせた値と推定価値の差を修正します． $\alpha$ はモンテカルロ法では状態 $s$ に遷移した回数の逆数とし，平均値を出しています．

TD法

例えば価値の更新を1ステップごとに行うときの更新式の一例は次式のようになります．

$V^{\pi}(s_t) \leftarrow V^{\pi}(s_t) + \alpha \bigl(R(s_t, a_t, s_{s+1}) + \gamma V^{\pi}(s_{t+1})- V^{\pi}(s_t) \bigr)$

この $R(s_t, a_t, s_{s+1}) + \gamma V^{\pi}(s_{t+1})- V^{\pi}(s_t)$ をTD誤差 (TD; Temporal Difference) といい，1ステップでの価値の更新分を示しています．

更新を1ステップごとに行う方法をTD(0)法，数ステップで行う方法をMulti-step learning，数ステップを組み合わせて更新する方法をTD( $\lambda$ )法と呼びます．

ここでは，TD(0)法の代表的な手法である「Q学習，SARSA, Actor-Critic」について説明します．

Q学習

ここでは行動価値関数 $Q^{\pi}(s,a)$ の推定を考えます．状態価値関数を行動価値関数に書き換え，価値が最大となる行動を選択すると仮定すると，上記の式は次式に書き換えることができます．

$Q^{\pi}(s_t,a_t) \leftarrow Q^{\pi}(s_t,a_t) + \alpha \bigl(R(s_t, a_t, s_{s+1}) + \gamma \max_a Q^{\pi}(s_{t+1}, a_{t+1}) - Q^{\pi}(s_t, a_t) \bigr)$

これがQ学習におけるQ値の更新式です．「価値が最大となる行動を選択する」という仮定を置き，方策を用いていないため「Off-policy（方策オフ型）」な手法と呼ばれます．

SARSA

Q学習とは異なり，行動は方策に従って選択されます．そのため更新式は次式のようになります．

$Q^{\pi}(s_{t},a_{t}) \leftarrow Q^{\pi}(s_{t},a_{t}) + \alpha \bigl( R(s_t, a_t, s_{s+1}) + \gamma Q^{\pi}(s_{t+1},a_{t+1}) - Q^{\pi}(s_t, a_t) \bigr)$

SARSA は，state-action-reward-(next)state-(next)actionの頭文字をとったものです．

方策を用いているので「On-policy（方策オン型）」な手法と呼ばれています．

（Q学習とSARSAの違い）

Q学習もSARSAもTD誤差を用いて価値を更新するので「value ベース」の手法と捉えることができます．（SARSAでは価値を更新することで方策を更新することにも繋がりますが，value ベースと言ってよいでしょう．）

Q学習は「Off-policy」，SARSAは「On-policy」です．Q学習は，常に価値が最大の行動をとるため，目的に向かって楽観的に最短経路で進みます．一方，SARSAは，方策に従って行動するため，慎重に安全な経路で進みます．

またQ学習では，更新している方策を使うという制約がないため，過去のデータを使ったり，他の方策を使ったりすることができ，サンプル数を増やすことができます．これはOn-policyなSARSAではできないメリットです．

Actor-Critic

Q学習やSARSAは「value ベース」の手法でした．これに対し，Actor-Critic は「value ベース」と「policy ベース」を組み合わせた手法です．

方策反復法でいう方策改善がActor，方策評価がCritic にあたります．方策と価値関数が別々に更新される（異なる関数を扱える）ため，例えば方策のほうをニューラルネットワークで近似するとすると，連続値を入力できるため，連続行動空間の問題にも使えるというメリットがあります．

2020-05-21

【強化学習】方策勾配法，方策勾配定理

強化学習

方策勾配法

強化学習の目的は「長期的な報酬を最大化する方策を求めること」でした．

方策勾配法は，方策をパラメータを持った関数として定義し，方策の価値が最大となるパラメータを勾配法で求める方法です．

パラメータ $\theta$ を持つ方策 $\pi_{\theta}$ の価値を $J(\theta)$ とします．方策の価値の勾配を $\nabla_{\theta} J(\theta)$ ，学習率を $\alpha$ とすると，パラメータの更新式は次式のようになります．

$\theta'=\theta+\alpha \nabla_{\theta} J(\theta)$

次に $\nabla_{\theta} J(\theta)$ の求め方について説明します．結果から示すと， $\nabla_{\theta} J(\theta)$ はこうなります．

$\nabla_{\theta} J(\theta) = \mathbb{E} _{d,\pi} [ \nabla_{\theta} \ln \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a)$ ]

これを方策勾配定理と呼びます．([Sutton+, 1999] Policy Gradient Methods for Reinforcement Learning with Function Approximation)

方策勾配定理の導出

方策の価値を状態価値関数とします．

$V^{\pi_{\theta}}(s) = \sum_{a} \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a)$

状態価値関数を $\theta$ で偏微分すると，

$\nabla_{\theta} V^{\pi_{\theta}}(s) = \sum_{a} \nabla_{\theta} \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a) + \sum_{a} \pi_{\theta}(a|s) \nabla_{\theta} Q^{\pi_{\theta}}(s,a)$

ここで，

$Q^{\pi_{\theta}}(s,a)=\sum_{s'} P(s'|s,a) \left( R(s,a,s') + \gamma V^{\pi_{\theta}}(s') \right)$

より

$\nabla_{\theta} Q^{\pi_{\theta}}(s,a) = \gamma \sum_{s'} P(s'|s,a) \nabla_{\theta} V^{\pi_{\theta}}(s')$

を代入して

$\nabla_{\theta} V^{\pi_{\theta}}(s) = \sum_{a} \nabla_{\theta} \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a) + \gamma \sum_{a} \pi_{\theta}(a|s) \sum_{s'} P(s'|s,a) \nabla_{\theta} V^{\pi_{\theta}}(s')$

$f(s) = \sum_{a} \nabla_{\theta} \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a)$ とすると，

$\nabla_{\theta} V^{\pi_{\theta}}(s) = f(s) + \gamma \sum_{a} \pi_{\theta}(a|s) \sum_{s'} P(s'|s,a) \nabla_{\theta} V^{\pi_{\theta}}(s')$

となり， $\nabla_{\theta} V^{\pi_{\theta}}(s)$ が再帰的な形になります．もう1ステップ進めると，

$\nabla_{\theta} V^{\pi_{\theta}}(s) = f(s) + \gamma \sum_{a} \pi_{\theta}(a|s) \sum_{s'} P(s'|s,a) \left(f(s') + \gamma \sum_{a'} \pi_{\theta}(a'|s') \sum_{s''} P(s''|s',a') \nabla_{\theta} V^{\pi_{\theta}}(s'') \right)$

ここで， $P(s \rightarrow s', t)$ を $t$ ステップで $s$ から $s'$ に遷移する確率とすると，

$\nabla_{\theta} V^{\pi_{\theta}}(s) = f(s)+\gamma P(s \rightarrow s', 1) f(s') + \gamma^2 P(s \rightarrow s'', 2)f(s'')+\cdots$

$= \sum_{s} \sum_{t=0}^{\infty} \gamma^t P(s \rightarrow s', t) f(s')$

$= \sum_{s'} \sum_{t=0}^{\infty} \gamma^t P(s \rightarrow s', t) \sum_{a} \nabla_{\theta} \pi_{\theta}(a|s') Q^{\pi_{\theta}}(s',a)$

これで，価値の勾配を示せることができました．もう少し正確には，方策の価値は「初期状態 $s_0$ で方策 $\pi_{\theta}$ に従って行動したときの価値」なので $J(\theta)=V^{\pi_{\theta}}(s_0)$ となります．

また， $\sum_{s'} \sum_{t=0}^{\infty} \gamma^t P(s \rightarrow s', t)$ を初期状態 $s_0$ から状態 $s$ へと遷移する確率に書き換えると， $\sum_{s} \sum_{t=0}^{\infty} \gamma^t P(s_0 \rightarrow s, t)$ となり，これを $\sum_s d^{\pi_{\theta}}(s)$ とします．

これより，方策の価値の勾配は次式で表されます．

$\nabla_{\theta} J(\theta)=\sum_s d^{\pi_{\theta}}(s) \sum_a \nabla_{\theta} \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a)$

$= \sum_s d^{\pi_{\theta}}(s) \sum_a \pi_{\theta}(a|s) \frac{\nabla_{\theta} \pi_{\theta}(a|s)}{\pi_{\theta}(a|s)} Q^{\pi_{\theta}}(s,a)$

$= \sum_s d^{\pi_{\theta}}(s) \sum_a \pi_{\theta}(a|s) \nabla_{\theta} \ln \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a)$

$= \mathbb{E} _{d,\pi} [ \nabla_{\theta} \ln \pi_{\theta}(a|s) Q^{\pi_{\theta}}(s,a)$ ]

これで，方策勾配定理が導出できました．

2020-05-09

【強化学習】モデルベース（価値反復法，方策反復法）

強化学習

強化学習の目的は「長期的な報酬を最大化する方策を求めること」であり，長期的な報酬は価値関数で表現されます．

ある状態 $s$ で方策 $\pi$ に従って行動したときの価値関数を状態価値関数と呼び，次の式で定義します．

$V^{\pi}(s) = \sum_a \pi(a|s) \sum_{s'} P(s'|s,a) \bigl(R(s,a,s') + \gamma V^{\pi}(s')\bigr)$

価値反復法

価値反復法 (Vlaue Iteration) とは，方策は「価値が最大となる行動を選択する」という仮定のもとで，価値を繰り返し更新しながら価値を推定する方法です．価値をベースに考えるので「Value ベース」の手法と呼ばれます．

方策は「価値が最大となる行動を選択する」ことなので，状態価値関数を以下の式に書き換えます．

$V^{\pi}(s) = \max_a \sum_{s'} P(s'|s,a) \bigl(R(s,a,s') + \gamma V^{\pi}(s')\bigr)$

この価値関数の形を用いて更新を行います．

具体的なアルゴリズムは以下の通りです．

Initialize $V(s)$ （例えば全ての状態に対して $V(s)=0$ ）

Repeat
　　 $\Delta$ ← 0
　　For each $s \in S$
　　　　 $v$ ← $V(s)$
　　　　 $V(s)$ ← $\max_a \sum_{s'} P(s'|s,a) \bigl(R(s,a,s') + \gamma V(s')\bigr)$
　　　　 $\Delta$ ← $\max(\Delta, |v-V(s)|)$
　　until $\Delta$ < $\theta$ （ $\theta$ は小さい正数）

方策反復法

方策反復法 (Policy Iteration) とは，価値反復法とは異なり，行動は方策に基づき選択されます．方策を更新していくことで方策を求めます．具体的には，方策を更新すると価値も更新されるため，更新された方策の下で価値を推定する「方策評価」と，推定した価値により方策を更新する「方策改善」を繰り返すことで方策を求めます．方策をベースに考えるので「policy ベース」と呼ばれます．

具体的なアルゴリズムは以下の通りです．

Initialize $V(s), \pi(a|s)$ （例えば全ての状態に対して $V(s)=0$ ， $\pi(a|s)$ は一様分布）

Policy Evaluation
　　Repeat
　　 $\Delta$ ← 0
　　For each $s \in S$
　　　　 $v$ ← $V(s)$
　　　　 $V(s)$ ← $\sum_a \pi(a|s) \sum_{s'} P(s'|s,a) \bigl(R(s,a,s') + \gamma V(s')\bigr)$
　　　　 $\Delta$ ← $\max(\Delta, |v-V(s)|)$
　　until $\Delta$ ← $\theta$ （ $\theta$ は小さい正数）

Policy Improvement
　　policy_stable ← True
　　For each $s \in S$
　　　　a ← $\pi(a|s)$
　　　　 $\pi(s)$ ← argmax $_a \sum_{s'} P(s'|s,a) \bigl(R(s,a,s') + \gamma V(s')\bigr)$
　　　　if $a \neq \pi(s)$
　　　　　　policy_stable ← False
　　if policy_stable, stop and return $V, \pi$ , else go to 2

2019-12-28

RのパッケージAERをインストールする際に詰まった問題の対処法

はじめに

普段はPython を使っているのですが，統計モデリング（個人的には離散選択モデリング）をする際にRのパッケージが充実しているのでRを使おうと思いました．

Sarrias, Mauricio, and Ricardo Daziano. 2017. “Multinomial Logit Models with Continuous and Discrete Individual Heterogeneity in R: The Gmnl Package.” Journal of Statistical Software, Articles 79 (2): 1–46. pdf

この論文に，離散選択モデル（多項ロジットモデルから混合ロジットモデル，潜在クラスロジットモデル，Mixed-Mixed Logit Model まで）を使用できるR のパッケージ "gmnl" の使い方が書いてあり，自分でも実装しようとしたところ，

使われているデータセットを呼び出すときのパッケージ "AER" をインストールする際に詰まってしまったのでブログにメモを残しておきます．

AER は，Applied Econometrics with R の略で応用計量経済学に特化したパッケージです．

警告メッセージ

エディタはJupyter Notebook を使っています．

バージョンは以下のとおりです．

> R.version

platform       x86_64-pc-linux-gnu         
arch           x86_64                      
os             linux-gnu                   
system         x86_64, linux-gnu           
status                                     
major          3                           
minor          4.4                         
year           2018                        
month          03                          
day            15                          
svn rev        74408                       
language       R                           
version.string R version 3.4.4 (2018-03-15)
nickname       Someone to Lean On

AER をインストールしてみます．すると警告メッセージがでました．

> install.packages("AER", dependencies = TRUE)

Warning message:
“dependencies ‘car’, ‘effects’, ‘mvtnorm’ are not available”
also installing the dependencies ‘ks’, ‘dynlm’, ‘rugarch’, ‘sampleSelection’, ‘systemfit’


Warning message in install.packages("AER", dependencies = TRUE):
“installation of package ‘ks’ had non-zero exit status”
Warning message in install.packages("AER", dependencies = TRUE):
“installation of package ‘dynlm’ had non-zero exit status”
Warning message in install.packages("AER", dependencies = TRUE):
“installation of package ‘systemfit’ had non-zero exit status”
Warning message in install.packages("AER", dependencies = TRUE):
“installation of package ‘AER’ had non-zero exit status”
Warning message in install.packages("AER", dependencies = TRUE):
“installation of package ‘rugarch’ had non-zero exit status”
Warning message in install.packages("AER", dependencies = TRUE):
“installation of package ‘sampleSelection’ had non-zero exit status”

どうやらcarというパッケージが使用できないようです．carというパッケージは，回帰分析などの手法に関する関数やデータを集めたものです．そこでcarをインストールしてみます．すると警告がでました．

> install.packages("car", dependencies = TRUE)

Warning message:
“package ‘car’ is not available (for R version 3.4.4)”

carは自分のRのバージョン 3.4.4 では使えないようです．

CRAN で確認すると，バージョン3.5以上が求められていました．

CRAN - Package car

CRAN (Comprehensive R Archive Network)は，R本体や各種パッケージをダウンロードするためのWebサイトのことです．（引用元：CRAN国内ミラーの使い方）

再びインストールにトライ

> install.packages("car", dependencies = TRUE)

also installing the dependencies ‘rematch’, ‘forcats’, ‘hms’, ‘readr’, ‘tidyselect’, ‘cellranger’, ‘progress’, ‘zip’, ‘estimability’, ‘highr’, ‘markdown’, ‘xfun’, ‘miniUI’, ‘webshot’, ‘DBI’, ‘sp’, ‘haven’, ‘data.table’, ‘readxl’, ‘openxlsx’, ‘minqa’, ‘nloptr’, ‘RcppEigen’, ‘effects’, ‘bdsmatrix’, ‘zoo’, ‘knitr’, ‘manipulateWidget’, ‘numDeriv’, ‘mitools’, ‘carData’, ‘abind’, ‘pbkrtest’, ‘quantreg’, ‘maptools’, ‘rio’, ‘lme4’, ‘alr4’, ‘coxme’, ‘leaps’, ‘lmtest’, ‘MatrixModels’, ‘rgl’, ‘sandwich’, ‘SparseM’, ‘survey’

インストールが完了し，carが使用可能になりました．（dependencies = TRUE にすると依存関係にあるパッケージも同時にインストールしてくれるのですが，時間がかかるのでFALSEにしておいても良いかもしれません）

そしてAERのインストールをします．（時間がかかりそうなのでAERのみインストール）

> install.packages("AER")

also installing the dependency ‘Formula’

インストールが完了したので呼び出してみます．

> library(AER)

Loading required package: lmtest

Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


Loading required package: sandwich

Loading required package: survival

必要なパッケージも呼び出して，使用できるようになりました．

データセットを呼びだしてみます．

> data("TravelMode", package = "AER")
> head(TravelMode)

f:id:ryonsd:20191228162222p:plain

これでデータが使えるようになりました．

まとめ

Rはほぼ使ったことがなかったのでパッケージのインストールに手こずりましたが，無事にパッケージを使用できるようになりました．上記の論文の写経を進めたいと思います．

参考文献

2019-12-17

確率分布（工事中）

統計学

確率分布とは，確率変数と確率の対応をとったものです．

ある事象がある値をとるときの確率は〜という言い方をします．
例えば，サイコロを振った時に出る目が1となる確率は $\frac{1}{6}$ という言い方を日常でも使うと思います．
このとき，サイコロを振ったときに出る目を確率変数，1を実現値と呼びます．

確率変数を $X$ ，実現値を $x$ ，確率を $p$ とすると，
確率変数 $X$ の実現値が $x$ となる確率は $p$ というのは次式で表現されます．

$Pr(X=x)=p$

サイコロの例だと， $Pr(X=1)=\frac{1}{6}$ です．

確率分布とは，確率変数と確率の対応をとったものです．

確率変数の実現値が離散か連続かで分類されます．

以下，離散確率分布と連続確率分布について書いています．

import math
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

離散確率分布

離散一様分布

ある事象がある実現値となる確率が起こり得る全ての実現値で同じ場合

def discrete_uniform(N):
    prob = 1. / N
    return prob

# 例えばサイコロの出る目
N = 6
x = np.arange(1, N+1)
probs = [discrete_uniform(N)] * N

plt.bar(x, probs)
plt.ylim(0, 1.)
plt.xlabel('X')
plt.ylabel('P(X)')

f:id:ryonsd:20200523210023p:plain

ベルヌーイ分布

確率変数の実現値が二つの場合の確率分布です．よく，ある事象が起こるか起こらないかを表すときに使われます．パラメータはある事象が起こる確率です．

def bernoulli(p):
    return p[f:id:ryonsd:20200523211401p:plain]

p = 0.5
x = [0,1]
probs = [1 - bernoulli(p), bernoulli(p)]

plt.bar(x, probs)
plt.ylim(0., 1.)
plt.xticks([0, 1])
plt.xlabel('X')
plt.ylabel('P(X)')

f:id:ryonsd:20200523211401p:plain

二項分布

ベルヌーイ分布（パラメータ $p$ ）に従う事象をある回数 $N$ 繰り返したときに，ある事象が起こる回数 $k$ が確率変数

def binominal(k, N, p):
    bin_coef = math.factorial(N) / math.factorial(k) / math.factorial(N-k)
    prob = bin_coef * (p ** k) * (1-p)**(N-k)
    return prob

N = 10
p = 0.5
x = np.arange(1, N+1)
probs = []
for k in x:
    probs.append(binominal(k, N, p))

plt.bar(x, probs)
plt.ylim(0., 1.)
plt.xlabel('X')
plt.ylabel('P(X)')

f:id:ryonsd:20200523211448p:plain

TD法

Q学習

SARSA

（Q学習とSARSAの違い）

Actor-Critic

方策勾配法

方策勾配定理の導出

価値反復法

方策反復法

はじめに

警告メッセージ

最新バージョンのRをインストール

再びインストールにトライ

まとめ

参考文献

離散確率分布

離散一様分布

ベルヌーイ分布

二項分布

カテゴリ分布

多項分布

ポアソン分布

負の二項分布

ポアソン二項分布

幾何分布

超幾何分布

ジップ分布

ゼータ分布

連続確率分布

連続一様分布

指数分布

正規分布（ガウス分布）

ベータ分布

ガンマ分布

ウィッシャート分布

カイ2乗分布

逆ガウス分布

対数正規分布

極値分布（generalized extreme value distribution, GEV)

ワイブル分布

ロジスティック分布

F分布

コーシー分布

アーラン分布

三角分布

ラプラス分布

レイリー分布

パレート分布

一般化双曲型分布

双曲線正割分布