ギャラクシースーパーはてなブログ

ギャラクシースーパーはてなブログ

ギャラクシースーパーノヴァ子の日記だお

【SolidWorks】寸法の一括変更

同じサイズの複数穴の径を一気に変更したいという場合に,今まで一つ一つ寸法変更しており,人生の貴重な時間をかなり無駄にしていました. 物事を独習で身に付けようとすると,不合理なフィッティングが生じてしまって,「どうしてその機能使わないの・・・?」という状況がしばしば起こります. 知ってる人からすれば当然のものばかりなのですが,たくさん機能があるソフトウェアだと,そのアイコンは認知してたけどなんかしらんが使ったことなかった...ということがたまにありますよね. 逐一身に付けていくしかない.

等しい値拘束を使いましょう

名無し電気設備(@nanashinanoni__)さんに教えていただきました.感謝. たとえば二つ以上の穴の径を同一にしたい場合は,ShiftまたはCtrlを押しながらクリックして穴を同時選択した状態で,「等しい値拘束」を入れましょう. 一つの寸法を変更すれば残りの寸法も一気に変更されます.

穴は穴ウィザードで作りましょう

ボルト穴は穴ウィザードを使おう.なんかしらんけど通し穴は今まで部品形状スケッチの中で描いていました.これはマジで無駄だった.「なんかしらんけど(不合理な行動)をしていた」というのがありすぎて泣きたくなる.

作って覚える SOLIDWORKSの一番わかりやすい本(Amazon)

【ロボティクス】運動学・ヤコビ行列・擬似逆行列の覚え書き

位置・姿勢に関する運動学

各関節にアクチュエータが搭載されたロボットアームを考える. 関節数を nとして,各関節変位を n次元ベクトル \boldsymbol{q}で表す. アームの手先位置と姿勢をベクトル \boldsymbol{x}で表す. \boldsymbol{x} m次元ベクトルとしておく.ふつう m=6であるが,運動が平面上に限られて,かつ姿勢角度を問わない場合, m=2となるといった状況もある.逆に双腕アームを有するロボットで,両手先の位置姿勢を指定したい場合は 6\times 2 = 12自由度になったりする.

ロボットアームを使用する上で,我々はアーム手先の位置・姿勢を直接指定したり,もしくはそれらが何らかの時間的な軌道を追従するように動かすことで,作業を行わせたい. 我々が直接操作できるのは各関節のアクチュエータの変位や速度といった値であるから,各関節を動かしたときに手先がどのように動くか,または手先に目的とする動作を行わせるためには各関節の変位や速度をどのように指定すればよいのかといった,関節-手先の関係を把握しておく必要がある.

手先位置・姿勢は各関節変位の関数として次のように表すことができる. \begin{align} \boldsymbol{x}(t)=\boldsymbol{f}\left(\boldsymbol{q}(t)\right) \end{align} これが順運動学であり,関数 \boldsymbol{f}は一般に非線形関数である.たとえば簡単なケースとして,2つの回転関節で構成された平面2自由度アームの手先位置を関節角 \theta _ 1 \theta _ 2で表そうとすると,これらを変数とする \sin \cosが出現することからイメージがつく. 順運動学はロボットの幾何学的な設計パラメータから求めることができる.隣接する関節同士の座標変換は簡単に行えるので,たとえばDH記法によって各座標系間のねじれや並進のずれを表現し,関節間の同次変換行列を求めることで任意の座標系から見た手先位置・姿勢を関節変位の関数として表すことが可能である.

逆運動学は \boldsymbol{f}の逆関数を手にすることであって,すなわち \begin{align} \boldsymbol{q}(t)=\boldsymbol{f}^{-1}\left(\boldsymbol{x}(t)\right) \end{align} という関係がわかればよい.しかしながら,先ほどの2自由度アームの例からわかるように,逆三角関数を用いて必要に応じて場合分けをしながら \boldsymbol{q}(t)を求めなければならず,これが6自由度くらいになると大変である. そこで解析的に逆運動学を解くのではなく,ニュートン法などの反復計算によって解を求めるということがふつう行われる.そこで役に立つのが微分運動学である.

微分運動学

手先自由度と関節数が等しい場合

手先と関節の(角)速度間の関係は線形で書ける.位置についての順運動学の式を時間微分し,表れた係数部分を行列としてまとめることで次の微分順運動学の式が得られる. \begin{align} \dot{\boldsymbol{x}}(t)=\boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t) \end{align} 行列 \boldsymbol{J}\left(\boldsymbol{q}(t)\right)はヤコビ行列である. 手先6自由度に対して関節数が6,といったように m=nとなっている状況であって,さらに特異姿勢でなければ,ヤコビ行列は正方かつ正則となるから,逆行列が存在して,微分逆運動学を求めることができる. \begin{align} \dot{\boldsymbol{q}}(t)=\boldsymbol{J}^{-1}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{x}}(t) \end{align} 説明順序が前後するが,特異姿勢とはヤコビ行列が非正則となる状態を意味する.この場合,ヤコビ行列の逆行列は存在しないから,目的とする手先速度を実現するための関節速度が得られず,ロボットは動作することができない.また,逆行列の各要素には行列式の逆数がかかることを思い起こすと,特異姿勢に近い状態では行列式の値がほぼ零であるため,逆行列の各要素は過大な値となる.したがって, \dot{\boldsymbol{x}}を微小に変化させるだけであっても \dot{\boldsymbol{q}}が極めて大きい値になり得る.関節速度が過大になるのは安全性やロボット本体にかかる負荷といった観点から考えてもよろしくない.

手先自由度と関節数が異なる場合

手先自由度 mと関節の数 nが異なる場合はどうすればよいのか. このときヤコビ行列は非正方となり,逆行列を考えることができない.冗長系,すなわち関節の数が多く m \lt nとなっている場合はヤコビ行列は横長の長方形となるし, m>nであれば縦長となる. 必要な手先自由度を実現するためには,関節数は自由度の数以上でなければならないから,ふつうは m \leq nとなるようにロボットアームを作る.運動学的にはこれで十分である.2次元平面上において手先位置を自在に指定したいのに,もし関節が一つしか無いのであれば,手先位置はリンク長で決まる円周上に限定されてしまう(手先位置は2次元ベクトルで表すことができるけれども,運動は1次元である). また,関節数 nはアクチュエータの個数とイコールであるとは限らないことに注意する. m \leq nとなっているロボットアームにおいても,全ての関節にアクチュエータが搭載されているとは限らない.関節数 n未満の個数のアクチュエータを用いる系を劣駆動系と呼ぶ.平面2自由度アームの二つの回転関節のうち,一つにしかアクチュエータが搭載されておらず,もう一つはフリージョイントであるとする.このような場合においても,手先位置は二つの関節角で決まる.すなわち冗長系と劣駆動系は両立する.一つのアクチュエータによって生成される運動によってフリージョイント側の関節角をどのように制御するか,というのはまた別の問題である.

関節数の方が少ない場合:最小二乗法

このときヤコビ行列は縦長となり,連立方程式の数(手先自由度)が未知数(関節速度)の数よりも多い.この場合,すべての方程式を満たす関節速度のセットは存在せず(不能), \dot{\boldsymbol{x}}(t)=\boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t)の関係がもはや成り立たない.そこで, \dot{\boldsymbol{x}}(t) \boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t)の差をなるべく小さくするような関節速度 \hat{\dot{\boldsymbol{q}}}(t)を望ましい答えだとみなすという方針を採る.すなわち最小二乗法を用いればよい.  \hat{\dot{\boldsymbol{q}}}(t)=\boldsymbol{J}^{+}\dot{\boldsymbol{x}}(t)と表すことができたとき,行列 \boldsymbol{J}^{+}を擬似逆行列と呼ぶ. 最小二乗法では差の二乗和を最小化する.ベクトルの大きさの二乗は自身との内積であるから,目的関数を \begin{align} Q=\left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right)^{T}\left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right) \end{align} とする.合成関数の微分によって \begin{align} \frac{\partial Q}{\partial \dot{\boldsymbol{q}}}&= \frac{\partial}{\partial \dot{\boldsymbol{q}}} \left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right) \frac{\partial}{\partial \left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right)} \left[ \left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right)^{T}\left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right) \right] \\ &=\boldsymbol{J}\left(\boldsymbol{q}\right)^{T} 2\left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right)\\ &=2\boldsymbol{J}\left(\boldsymbol{q}\right)^{T}\left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right) \end{align} これがゼロになるときの \dot{\boldsymbol{q}}を求めればよいから, \begin{align} 2\boldsymbol{J}\left(\boldsymbol{q}\right)^{T}\left(\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right) &= 0 \\ \boldsymbol{J}\left(\boldsymbol{q}\right)^{T}\boldsymbol{J}\left(\boldsymbol{q}\right)\dot{\boldsymbol{q}} &=\boldsymbol{J}\left(\boldsymbol{q}\right)^{T} \dot{\boldsymbol{x}} \end{align} ここでグラム行列 \boldsymbol{J}\left(\boldsymbol{q}\right)^{T}\boldsymbol{J}\left(\boldsymbol{q}\right) n\times nの正方行列であり, \boldsymbol{J}が列フルランクであれば, {\rm rank}(\boldsymbol{J}^{T}\boldsymbol{J})=nであって逆行列を持つ. したがって \begin{align} \dot{\boldsymbol{q}} &=\left(\boldsymbol{J}\left(\boldsymbol{q}\right)^{T}\boldsymbol{J}\left(\boldsymbol{q}\right)\right)^{-1}\boldsymbol{J}\left(\boldsymbol{q}\right)^{T} \dot{\boldsymbol{x}} \end{align} としてそれらしい \dot{\boldsymbol{q}}を得ることができた.このときの \left(\boldsymbol{J}\left(\boldsymbol{q}\right)^{T}\boldsymbol{J}\left(\boldsymbol{q}\right)\right)^{-1}\boldsymbol{J}\left(\boldsymbol{q}\right)^{T} が, \boldsymbol{J}が縦長で列フルランクであるときの擬似逆行列である. これは順運動学を満たしている(順運動学の式の両辺に左から擬似逆行列を作用させてみればよい).

関節数の方が多い場合:ラグランジュ未定乗数法

この場合はヤコビ行列が横長となり,連立方程式の数(手先自由度)よりも未知数(関節速度)の数が多い.この連立方程式をがんばって解こうとしても,どうしても未知数が定まらないまま残ってしまう(不定). つまり等式 \dot{\boldsymbol{x}}(t)=\boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t)を満たす \dot{\boldsymbol{q}}が無数に存在する. このように解が不定の状況では,無数に存在する解に何らかの制約を課すことで一つに絞り込むことができる. そこでノルム |\dot{\boldsymbol{q}}|を最小化することを考える.ただし,このとき等式 \dot{\boldsymbol{x}}(t)=\boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t)が満たされていなければならない. すなわち等式制約付きの最小化問題を解くことになる.ラグランジュ未定乗数法によって目的関数は \begin{align} Q=\dot{\boldsymbol{q}}^{T}\dot{\boldsymbol{q}}+\boldsymbol{\lambda}^{T} \left( \dot{\boldsymbol{x}}(t)-\boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t) \right) \end{align} と書ける. \begin{align} \frac{\partial Q}{\partial \dot{\boldsymbol{q}}}= 2\dot{\boldsymbol{q}} - \boldsymbol{J}^{T}\boldsymbol{\lambda} &= 0 \\ \dot{\boldsymbol{q}} &= \frac{1}{2} \boldsymbol{J}^{T}\boldsymbol{\lambda}\\ \end{align} となるから, \dot{\boldsymbol{x}}(t)=\boldsymbol{J}\left(\boldsymbol{q}(t)\right)\dot{\boldsymbol{q}}(t)= \frac{1}{2} \boldsymbol{J} \boldsymbol{J}^{T}\boldsymbol{\lambda}. 行列 \boldsymbol{J}\left(\boldsymbol{q}\right)\boldsymbol{J}\left(\boldsymbol{q}\right)^{T} m\times mの正方行列であり, \boldsymbol{J}が行フルランクであれば, {\rm rank}(\boldsymbol{J}\boldsymbol{J})^{T}=mであって逆行列を持つ.したがって,  \boldsymbol{\lambda}=2\left( \boldsymbol{J} \left(\boldsymbol{q}\right) \boldsymbol{J}\left(\boldsymbol{q}\right)^{T} \right)^{-1}\dot{\boldsymbol{x}}となる.これにより, \begin{align} \dot{\boldsymbol{q}} = \frac{1}{2} \boldsymbol{J}^{T}\boldsymbol{\lambda} = \frac{1}{2} \boldsymbol{J}^{T}2\left( \boldsymbol{J} \boldsymbol{J}^{T} \right)^{-1}\dot{\boldsymbol{x}} =\boldsymbol{J}^{T}\left( \boldsymbol{J} \boldsymbol{J}^{T} \right)^{-1}\dot{\boldsymbol{x}} \end{align} を得る.行列 \boldsymbol{J}^{T}\left( \boldsymbol{J} \boldsymbol{J}^{T} \right)^{-1}が, \boldsymbol{J}が横長で行フルランクであるときの擬似逆行列である.

\begin{align} \dot{\boldsymbol{q}} =\boldsymbol{J}^{T}\left( \boldsymbol{J} \boldsymbol{J}^{T} \right)^{-1}\dot{\boldsymbol{x}} \end{align} の両辺に左から \boldsymbol{J}を作用させると,順運動学の式となることがわかる.

特異姿勢となってランク落ちした場合

階数因数分解による方法

どちらもランクが r m\times r行列 \boldsymbol{B} r\times n行列 \boldsymbol{C}を用いて \boldsymbol{J}=\boldsymbol{B}\boldsymbol{C}と階数因数分解する.  \boldsymbol{J}はランク落ちしており,これまでに挙げた手法では擬似逆行列を計算することができない.グラム行列のランクは元の行列のランクに等しく,つまり {\rm rank}(\boldsymbol{J})={\rm rank}(\boldsymbol{J}^{T}\boldsymbol{J})={\rm rank}(\boldsymbol{J}\boldsymbol{J}^{T})だから, \boldsymbol{J}がランク落ちならば \boldsymbol{J}^{T}\boldsymbol{J} \boldsymbol{J}\boldsymbol{J}^{T}もランク落ちとなり,これらの逆行列が存在しないからである.しかしながら, \boldsymbol{B}は列フルランク, \boldsymbol{C}は行フルランクであるから,これらは擬似逆行列を持つ. したがって \begin{align} \boldsymbol{J}^{+}=\left(\boldsymbol{BC}\right)^{+}=\boldsymbol{C}^{+}\boldsymbol{B}^{+}=\boldsymbol{C}^{T}\left( \boldsymbol{C} \boldsymbol{C}^{T} \right)^{-1}\left(\boldsymbol{B}^{T}\boldsymbol{B}\right)^{-1}\boldsymbol{B}^{T} \end{align} と \boldsymbol{J}の擬似逆行列を計算することができる.  \left(\boldsymbol{BC}\right)^{+}=\boldsymbol{C}^{+}\boldsymbol{B}^{+} \boldsymbol{B} \boldsymbol{C}のランクが等しいとき成立する.

SR逆行列(Singularity-Robust Inverse)

手先速度を所望の値にすることよりも,特異点付近での安定性を優先する場合,最小化する目的関数を次のように変更すればよい.

\begin{align} Q=\left(\boldsymbol{J}\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right)^{T}\left(\boldsymbol{J}\dot{\boldsymbol{q}}-\dot{\boldsymbol{x}}\right) +\lambda \dot{\boldsymbol{q}}^{T}\dot{\boldsymbol{q}} \end{align}

そうすると結局, \begin{align} \dot{\boldsymbol{q}}=\left( \boldsymbol{J}^{T}\boldsymbol{J} + \lambda \boldsymbol{I}\right)^{-1} \boldsymbol{J}^{T} \dot{\boldsymbol{x}} \end{align} となる. \boldsymbol{J}^{T}\boldsymbol{J}がランク落ちであっても,  \lambda \boldsymbol{I}のおかげで \left( \boldsymbol{J}^{T}\boldsymbol{J} + \lambda \boldsymbol{I}\right)^{-1} の逆行列は求まる. 減衰因子 \lambdaは可操作度を用いるなどして,特異姿勢への近さに応じて調整する.

とりあえず一旦ここまで

擬似逆行列については,特異値分解との関係もまとめたい.可操作度を考えるにあたっても必要だし.

参考書籍

  1. ロボット制御基礎論(Amazon) 冗長アームといえばこれ.
  2. 実践ロボット制御(Amazon) ロボット制御基礎論をベースにまとめられた新しい本.ロボットアームの入門的内容を概観するのであればこの本から始めてもよいと思う. ただ当たり前のことだが一冊で勉強が完結するということはないので,ロボット制御基礎論などの他の本を一緒に持っておくのがよい.
  3. “巧みさ”とロボットの力学 プレミアムブックス版(Amazon) 冗長アームやハンドについて,いかに器用な動作を実現するか?というテーマで2008年に書かれた本.深層学習以後,有本先生たちがどのようなことを考えられているのかは,2017年のロボット制御学ハンドブック(Amazon)最終章にまとめられている.
  4. ヒューマノイドロボット(改訂2版)(Amazon) ヒューマノイドを題材にして逆運動学の数値解法やMATLAB実装例が紹介されており,イメージがつきやすい.

【現代制御】可制御性

はじめに

現代制御,面白いですよね. と言っても私は大人になってから社会人の教養として雰囲気を嗜んだだけなので,教科書を読みながら(フーン...そうなんだ...)程度の感想とともにスルーしてしまっているところが結構あります. 大人になってから独学したので...と言い訳しましたが,それ自体は理解の精度に対して本質的ではなく,仮に大学で履修していたとしても似たようなものだっただろうという気もします.

現代制御では線形時不変システムの可制御性を,可制御行列と呼ばれる正方行列の正則性を調べることによって判定することができます. 教科書によっては「可制御行列がフルランクであれば可制御」,もしくは「可制御行列の行列式が非零であれば可制御」といった書かれ方をされている本がありますが,どれも同じことです(えっ,という方は線形代数の教科書を読み返しましょう). 可制御行列を用いて可制御性が判定できることの証明は,易しい入門書では説明が省略されていることもあります.備忘録としてまとめておきます.

参考書籍

例題で学ぶ現代制御の基礎(Amazon)

この本は必要な情報のまとまりも良く,個人的に読みやすいと感じています.

線形時不変システムの状態空間表現

さて,1入力1出力である線形時不変システムを考える. 状態方程式と出力方程式は \begin{align} \dot{\boldsymbol{x}}(t)&=\boldsymbol{A}\boldsymbol{x}(t)+\boldsymbol{b}u(t) \\ y(t)&=\boldsymbol{c}^T\boldsymbol{x}(t) \end{align} である.ここで状態 \boldsymbol{x}(t) \boldsymbol{b}および \boldsymbol{c} n次元のベクトルであって, \boldsymbol{A} n\times n行列. 状態方程式の解は \begin{align} \boldsymbol{x}(t)=e^{\boldsymbol{A}t} \boldsymbol{x}_{0}+ \int_0 ^{t} e^{\boldsymbol{A}(t-\tau)} \boldsymbol{b} u(\tau) d\tau \end{align} です.

可制御性とその判定方法

システムが可制御であるとは,制御入力を与えることによって有限の時間で任意の初期状態から目標とする状態へ変化させられることを意味する. 可制御であるかどうかは, n\times nの正方行列 
\begin{bmatrix}
\boldsymbol{b} \ \boldsymbol{Ab} \ \cdots \ \boldsymbol{A}^{n-1}\boldsymbol{b}
\end{bmatrix}
が正則であるかどうかによって判定可能. この行列を可制御行列と呼ぶ.

なぜ可制御性の判定となるのか

必要性

システムが可制御であれば,可制御行列が正則となることを示す. 可制御であるシステムでは,適当な制御入力 u(t)を与えることで,ある時刻 t_fにて \boldsymbol{x}(t_f)=\boldsymbol{x}_fとなるから,

\begin{align} \boldsymbol{x}(t _ f) &=\boldsymbol{x} _ f =e^{\boldsymbol{A}t _ f} \boldsymbol{x} _ 0 + \int _ 0 ^{t _ f} e^{\boldsymbol{A}(t _ f-\tau)} \boldsymbol{b} u(\tau) d\tau \\ e^{\boldsymbol{A}t _ f} \boldsymbol{x} _ 0 - \boldsymbol{x} _ f &=-\int _ 0 ^{t _ f} e^{\boldsymbol{A}(t _ f-\tau)} \boldsymbol{b} u(\tau) d\tau \end{align}

ここで,行列指数関数の性質として, \left(e^{\boldsymbol{A}t _ f} \right)^{-1}= e^{-\boldsymbol{A}t _ f}であることに注意する.これを使って上式を書き換えると, \begin{align} \boldsymbol{x} _ 0 - e^{-\boldsymbol{A}t _ f}\boldsymbol{x} _ f &=-e^{-\boldsymbol{A}t _ f}\int _ 0 ^{t _ f} e^{\boldsymbol{A}(t _ f-\tau)} \boldsymbol{b} u(\tau) d\tau \\ &=-\int _ 0 ^{t _ f} e^{-\boldsymbol{A}t _ f+\boldsymbol{A}(t _ f-\tau)} \boldsymbol{b} u(\tau) d\tau \\ &=-\int _ 0 ^{t _ f} e^{-\boldsymbol{A}\tau} \boldsymbol{b} u(\tau) d\tau \end{align} となる.

また,行列指数関数は適当な時間関数を係数として n個の項の和に展開できる. \begin{align} e^{\boldsymbol{A}t} = r _ 0 (t) \boldsymbol{I} + r _ 1 (t) \boldsymbol{A} + \cdots + r _ {n-1}(t)\boldsymbol{A}^{n-1} \\ e^{-\boldsymbol{A}t} = r _ 0 (-t) \boldsymbol{I} + r _ 1 (-t) \boldsymbol{A} + \cdots + r _ {n-1}(-t)\boldsymbol{A}^{n-1} \end{align} これを使って \begin{align} \boldsymbol{x} _ 0 - e^{-\boldsymbol{A}t _ f}\boldsymbol{x} _ f &=-\int _ 0 ^{t _ f} \left(r _ 0 (-\tau) \boldsymbol{I} + r _ 1 (-\tau) \boldsymbol{A} + \cdots + r _ {n-1}(-\tau)\boldsymbol{A}^{n-1} \right) \boldsymbol{b} u(-\tau) d\tau \\ &= \left( -\int _ 0 ^{t _ f} r _ 0 (-\tau) \boldsymbol{b} u(\tau) d\tau \right) + \cdots + \left( -\int _ 0 ^{t _ f} r _ {n-1} (-\tau) \boldsymbol{A}^{n-1} \boldsymbol{b} u(\tau) d\tau \right) \\ &= s _ 0 \boldsymbol{b}+ \cdots + s _ {n-1} \boldsymbol{A}^{n-1} \boldsymbol{b} \\ &= \begin{bmatrix} \boldsymbol{b} & \boldsymbol{A}\boldsymbol{b} & \cdots & \boldsymbol{A}^{n-1}\boldsymbol{b} \end{bmatrix} \begin{bmatrix} s _ 0 \\ s _ 1 \\ \vdots \\ s _ {n-1} \end{bmatrix} = \boldsymbol{U}_c \boldsymbol{s} \end{align} と書ける.ただし,途中で \begin{align} s _ i = -\int _ 0 ^{t _ f} r _ i (-\tau) u(\tau) d\tau \end{align} とおいた.

\begin{align} \boldsymbol{x} _ 0 - e^{-\boldsymbol{A}t _ f}\boldsymbol{x} _ f &= \boldsymbol{U}_c \boldsymbol{s} \end{align} ここで登場した \boldsymbol{A} \boldsymbol{b}で構成されている行列 \boldsymbol{U} _ cが可制御行列である. 上式の左辺は n次元ベクトルであり, \boldsymbol{s}も同じく n次元ベクトルであるから,この可制御行列はフルランクでなければならない. (可制御行列が非正則,すなわちランク落ちしていると,右辺の次元は nよりも小さくなり等式が成立しない.)

十分性

可制御行列が正則であれば,システムが可制御であることを示す.

(つづく)

【モデル予測制御】ステップ応答のモデル推定値を有限個の和で表現する

目的

モデル予測制御の学習のために参考資料を読んでいたところ,サラッと出て来たステップ応答の有限和の表現(式(10))で躓いてしまった.備忘のために式の導出のノートを残しておく.

※というわけで,モデル予測制御と銘打っていますがモデル予測制御のモの字も出て来ません.私も悲しいのです.あしからず.

参考資料

概要

  • 足して引く(正味ゼロ)操作によって和の書き換えを行う.
  • 入力が加えられてから,時間 s以上経過した後の応答は一定値に収束すると仮定する.事前に定めたこの仮定によって, \sum_{k=s}^\inftyの無限和の項を消す.

本題

時刻 t+jでのステップ応答として,次に示す式(7)を得る.

 y_M(t+j)=\sum_{k=0}^{j-1}a_{j-k}\Delta u(t+k)+\sum_{k=1}^{\infty}a_{j+k}\Delta u(t-k)

ここで第1項が現在時刻 t以降の入力による応答であり,第2項は時刻 t-1以前の過去の入力による応答である. モデル予測制御では制御入力区間を定めて,時刻t+1以降の入力のうち,時刻t+N_M-1までの操作量u(t+1), \cdots, u(t+N_M-1)を求める.t+N_M以降は変化させずにu(t+N_M-1)のままとする. すなわち,  \Delta u(t+i)=0 \hspace{15pt} {\rm for}\hspace{5pt} i \geq N_M 以上から式(7)は次のように書き直せる.(式9)  y_M(t+j)=\sum_{k=0}^{N_M-1}a_{j-k}\Delta u(t+k)+\sum_{k=1}^{\infty}a_{j+k}\Delta u(t-k) モデルによる将来の推定値y_M(t+j)を計算するためには無限個の過去の操作量が必要となってしまう.現実のアルゴリズムとして実装するために,第2項の無限個の和を有限個の和に書き直したい. 第2項は次のように書き直せる.

\begin{align} \sum _ {k=1}^{\infty}a _ {j+k}\Delta u(t-k) &= \sum _ {k=1}^{s-1}a _ {j+k}\Delta u(t-k) + \sum _ {k=s}^{\infty}a _ {j+k}\Delta u(t-k) \\ &= \sum _ {k=1}^{s-1}(a _ {j+k}-a _ k+a _ k)\Delta u(t-k) + \sum _ {k=s}^{\infty}a _ {j+k}\Delta u(t-k) \\ &= \sum _ {k=1}^{s-1}a _ k\Delta u(t-k) + \sum _ {k=1}^{s-1}(a _ {j+k}-a _ k)\Delta u(t-k) + \sum _ {k=s}^{\infty}a _ {j+k}\Delta u(t-k) \\ &= \sum _ {k=1}^{s-1}a _ k\Delta u(t-k) + \sum _ {k=1}^{s-1}(a _ {j+k}-a _ k)\Delta u(t-k) + \sum _ {k=s}^{\infty}a _ {j+k}\Delta u(t-k) + \sum _ {k=s}^{\infty}a _ {k}\Delta u(t-k) - \sum _ {k=s}^{\infty}a _ {k}\Delta u(t-k) \\ &= \left (\sum _ {k=1}^{s-1}a _ k \Delta u(t-k) + \sum _ {k=s}^{\infty}a _ {k}\Delta u(t-k) \right) + \sum _ {k=1}^{s-1}(a _ {j+k}-a _ k)\Delta u(t-k) + \left( \sum _ {k=s}^{\infty}a _ {j+k}\Delta u(t-k) - \sum _ {k=s}^{\infty}a _ {k}\Delta u(t-k) \right) \\ &= \sum _ {k=1}^{\infty}a _ k \Delta u(t-k) + \sum _ {k=1}^{s-1}(a _ {j+k}-a _ k)\Delta u(t-k) + \sum _ {k=s}^{\infty}(a _ {j+k}-a _ {k})\Delta u(t-k) \\ &= y _ M(t) + \sum _ {k=1}^{s-1}(a _ {j+k}-a _ k)\Delta u(t-k) + \sum _ {k=s}^{\infty}(a _ {j+k}-a _ {k})\Delta u(t-k) \\ \end{align}

ここで,上の式の第3項は時刻tからsステップ以上過去の入力による成分である. sステップ経過すると,応答は一定値になると仮定した.すなわち,  a_i = a_s \hspace{15pt} {\rm for} \hspace{5pt} i \geq s であった.第3項の応答係数はそれぞれ

\begin{align} &a _ {j+s}, a _ {j+s+1}, a _ {j+s+2}, \cdots \\ &a _ {s}, a _ {s+1}, a _ {s+2}, \cdots \end{align}

であり,添え字は全てs以上であるから,応答係数は全てa_sになる.したがって,差a_{j+k}-a_{k}=0となり,第3項は消える. 以上を踏まえて,計算結果を元の式(9)に戻すと,  y_M(t+j)=\sum_{k=0}^{N_M-1}a_{j-k}\Delta u(t+k)+ y_M(t) + \sum_{k=1}^{s-1}(a_{j+k}-a_k)\Delta u(t-k) となって,これは式(10)に一致した.

多次元正規分布の条件付き分布

はじめに

ガウス過程回帰を導出する上で必要になる多次元正規分布の条件付き分布についてまとめておく. 教科書「ガウス過程と機械学習」を参考に,式変形の各ステップをなるべく省略せずに記した.

参考資料

持橋先生,大場先生の「ガウス過程と機械学習」先行公開 (γ2版)の第2章. 19年3月発売予定らしいが,サポートページにて一部公開されている.先行公開原稿が素晴らしいのに加えて,サポートページの内容の充実っぷりがすごい. ただ,(私の勘違いかもしれないが)19年2/4現在において,公開版では後述のように微妙な誤りがあるので注意.正式版では修正されることを期待.

19/03/12追記:書籍版では一部修正されていました.また,正誤表がサポートページにて公開されています.

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

多次元正規分布

 D次元のベクトル \boldsymbol{x}=\left(x _ 1, \cdots, x _ D \right)が平均 \boldsymbol{\mu},共分散行列 \boldsymbol{\Sigma}の正規分布に従うとき,以下のように表す.

\begin{align} \boldsymbol{x}&\sim \mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}, \boldsymbol{\Sigma} \right) \\ \mathcal{N}(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})&=\frac{1}{\left(\sqrt{2\pi}^D \sqrt{|\boldsymbol{\Sigma}|}\right)}\exp\left(-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^T \boldsymbol{\Sigma}^{-1} \left(\boldsymbol{x}-\boldsymbol{\mu}\right)\right) \end{align}

多次元正規分布の条件付き分布

 \boldsymbol{x}を二つのベクトル \boldsymbol{x}\ _ 1 \boldsymbol{x}\ _ 2に分ける.  \boldsymbol{x}=\left(x\ _ 1, \cdots, x\ _ D \right)から,最初の L次元 \boldsymbol{x}\ _ 1=\left(x\ _ 1, \cdots, x\ _ L \right)を抜き出し,残りを \boldsymbol{x}\ _ 2=\left(x\ _ {L+1}, \cdots, x _ D \right)とする.

このとき \boldsymbol{x} _ 1を固定したときの条件付き分布 p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right)は次のように書ける.

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right)=\mathcal{N}\left(\boldsymbol{\mu} _ 2+\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1}\left(\boldsymbol{x} _ 1-\boldsymbol{\mu} _ 1\right), \boldsymbol{\Sigma} _ {22}-\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1}\boldsymbol{\Sigma} _ {12}\right) \ \end{align}

導出

同時分布と条件付き分布の関係(乗法定理)は

\begin{align} p\left(\boldsymbol{x} _ 1,\boldsymbol{x} _ 2 \right)=p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right)p\left(\boldsymbol{x} _ 1 \right)\end{align}

だった.条件付き分布は

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right)=\frac{p\left(\boldsymbol{x} _ 1,\boldsymbol{x} _ 2 \right)}{p\left(\boldsymbol{x} _ 1\right)} \ \end{align}

であり,今 \boldsymbol{x} _ 1は条件として固定されているので, \boldsymbol{x} _ 2の関数 p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) p\left(\boldsymbol{x} _ 2,\boldsymbol{x} _ 1 \right)に比例している(分母の p\left(\boldsymbol{x} _ 1 \right)には依存しない). すなわち p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) \propto p\left(\boldsymbol{x} _ 1,\boldsymbol{x} _ 2 \right)である.

さて, \boldsymbol{x}\sim \mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}, \boldsymbol{\Sigma} \right) であった.同時分布 p\left(\boldsymbol{x} _ 2,\boldsymbol{x} _ 1 \right) \boldsymbol{x} _ 1 \boldsymbol{x} _ 2が「同時に」得られる確率なのだから, p\left(\boldsymbol{x} _ 2,\boldsymbol{x} _ 1 \right)=\mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}, \boldsymbol{\Sigma} \right) である.

 \boldsymbol{x}を二つのベクトルに分割したので, \boldsymbol{\mu} \boldsymbol{\Sigma}も二つに分割して,次のように表すことができる.

\begin{align} \boldsymbol{x}= \left( \begin{array}{c} \boldsymbol{x} _ {1} \\ \boldsymbol{x} _ {2} \end{array} \right) \sim \mathcal{N}\left( \left( \begin{array}{c} \boldsymbol{\mu} _ {1} \\ \boldsymbol{\mu} _ {2} \end{array} \right) , \left( \begin{array}{cc} \boldsymbol{\Sigma} _ {11} & \boldsymbol{\Sigma} _ {12}\\ \boldsymbol{\Sigma} _ {21} & \boldsymbol{\Sigma} _ {22} \end{array} \right) \right) \end{align}

精度行列,すなわち共分散行列の逆行列 \Lambdaを次のように定義する.

\begin{align} \Lambda= \begin{pmatrix} \boldsymbol{\Lambda} _ {11} & \boldsymbol{\Lambda} _ {12} \\ \boldsymbol{\Lambda} _ {21} & \boldsymbol{\Lambda} _ {22} \end{pmatrix} = \begin{pmatrix} \boldsymbol{\Sigma} _ {11} & \boldsymbol{\Sigma} _ {12} \\ \boldsymbol{\Sigma} _ {21} & \boldsymbol{\Sigma} _ {22} \end{pmatrix}^{-1} \end{align}

この精度行列を用いることによって,元の正規分布を次のように表せる.

\begin{align} \mathcal{N}\left(\boldsymbol{x}|\boldsymbol{\mu}, \boldsymbol{\Sigma} \right) = p\left(\boldsymbol{x} _ 2,\boldsymbol{x} _ 1 \right) &\propto \exp\left(-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^T \boldsymbol{\Sigma}^{-1} \left(\boldsymbol{x}-\boldsymbol{\mu}\right)\right) \\ &= \exp\left(-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^T \boldsymbol{\Lambda} \left(\boldsymbol{x}-\boldsymbol{\mu}\right)\right) \\ &= \exp\left(-\frac{1}{2} \begin{pmatrix} \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \\ \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \end{pmatrix}^T \begin{pmatrix} \boldsymbol{\Lambda} _ {11} & \boldsymbol{\Lambda} _ {12}\\ \boldsymbol{\Lambda} _ {21} & \boldsymbol{\Lambda} _ {22} \end{pmatrix} \begin{pmatrix} \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \\ \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \end{pmatrix} \right) \end{align}

次に \expの括弧の中を展開したいのだが,要素がブロックに分割されたベクトルの転置については注意しておこう.次に示すように,中身を転置した上で,さらに各ブロックについて転置を取る必要がある.

\begin{align} \begin{pmatrix} \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \\ \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \end{pmatrix}^T = \begin{pmatrix} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T & \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \end{pmatrix} \end{align}

したがって,

\begin{align} \begin{pmatrix} \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \\ \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \end{pmatrix}^T \begin{pmatrix} \boldsymbol{\Lambda} _ {11} & \boldsymbol{\Lambda} _ {12}\\ \boldsymbol{\Lambda} _ {21} & \boldsymbol{\Lambda} _ {22} \end{pmatrix}\begin{pmatrix} \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \\ \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \end{pmatrix}&= \begin{pmatrix} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T & \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \end{pmatrix}\begin{pmatrix} \boldsymbol{\Lambda} _ {11} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) +\boldsymbol{\Lambda} _ {12} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \\ \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) +\boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \end{pmatrix} \\ &= \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {11} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) + \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {12} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)+ \\ & \hspace{14pt} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) + \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \\ &= \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {11} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) + 2\left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)+ \\ & \hspace{14pt} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \\ \end{align}

最後の式変形では, \boldsymbol{\Lambda}\ _ {12}=\boldsymbol{\Lambda}\ _ {21}であり,中央の二つの項が同一であることから導かれる. ここで一旦, \expの中に戻して眺めてみる. \expの括弧内の和は,それぞれの \expの積に分解できるから,

\begin{align} &\exp\left( \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {11} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) + 2\left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) + \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \right)\\ &= \exp\left( \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {11} \left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right) \right) \exp\left( 2\left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) + \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \right)\\ \end{align}

となる. -\frac{1}{2}を省略していることに注意.

ここで \exp\left(
\left( \boldsymbol{x}\ _ {1} - \boldsymbol{\mu}\ _ {1} \right)^T
\boldsymbol{\Lambda}\ _ {11} \left( \boldsymbol{x}\ _ {1} - \boldsymbol{\mu}\ _ {1} \right)
\right)には \boldsymbol{x} _ 2 が含まれておらず, p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right)はこの項に依存しない.したがって引き続き比例関係のみに注目すれば,

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) &\propto p\left(\boldsymbol{x} _ 1,\boldsymbol{x} _ 2 \right)\\ &\propto \exp\left( -\frac{1}{2} \left( 2\left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) + \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \right) \right) \end{align}

となる. \expの括弧の中身をさらに展開しよう.

\begin{align} &2\left( \boldsymbol{x} _ {1} - \boldsymbol{\mu} _ {1} \right)^T \boldsymbol{\Lambda} _ {21} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) + \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\mu} _ {2} \right) \\ &= \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} - \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} - \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} + \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} + \\ &\hspace{15pt} 2\boldsymbol{x} _ {1}^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} -2\boldsymbol{x} _ {1}^T \boldsymbol{\Lambda} _ {21} \boldsymbol{\mu} _ {2} -2\boldsymbol{\mu} _ {1}^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} +2\boldsymbol{\mu} _ {1}^T \boldsymbol{\Lambda} _ {21} \boldsymbol{\mu} _ {2} \end{align}

ここで, \boldsymbol{x}\ _ {2}が含まれていない項は,先ほどと同様に独立な \expの項として分離することができ,ただの係数となって比例関係から無視することができる.そうすると生き残る項は

\begin{align} &\boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} - \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} - \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} + 2\boldsymbol{x} _ {1}^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} -2\boldsymbol{\mu} _ {1}^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} \\ &= \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} - \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} - \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} + 2\left(\boldsymbol{x} _ {1}^T-\boldsymbol{\mu} _ {1}^T\right) \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} \\ &= \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} - \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} - \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} + 2\left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right)^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} \end{align}

ここで,各項はそれぞれ内積なので,転置を取っても値が変わらない.また,精度行列は対称行列であることから,

\begin{align} \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} &= \left( \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} \right)^T \\ &= \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} \\ 2\left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right)^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} &= \left( 2\left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right)^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} \right)^T \\ &= 2\boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) \end{align}

となる.そうすると上の式は,次のように書ける.

\begin{align} &\boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} - \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} - \boldsymbol{\mu} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} + 2\left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right)^T \boldsymbol{\Lambda} _ {21} \boldsymbol{x} _ {2} \\ &= \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} - 2\boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} + 2\boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) \\ &= \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} -2\boldsymbol{x} _ {2}^T \left( \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) \right) \end{align}

見やすくするためにベクトル \boldsymbol{\Lambda}\ _ {22} \boldsymbol{\mu}\ _ {2}
+\boldsymbol{\Lambda}\ _ {21} \left(\boldsymbol{x}\ _ {1}-\boldsymbol{\mu}\ _ {1}\right) \boldsymbol{a}とおく. そうして上の式を平方完成すると,

\begin{align} &\boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} -2\boldsymbol{x} _ {2}^T \boldsymbol{a} = \left( \boldsymbol{x} _ {2} - \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{a} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{a} \right) - \boldsymbol{a}^T \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{a} \end{align}

ここで, \boldsymbol{a}^T \boldsymbol{\Lambda}\ _ {22}^{-1} \boldsymbol{a} \boldsymbol{x}\ _ {2}を含んでいないから,これまでと同様に独立な \expの項として分離することができ,ただの係数となって比例関係から無視することができる.

結局,比例関係は次のようになる.

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) &\propto p\left(\boldsymbol{x} _ 1,\boldsymbol{x} _ 2 \right)\\ &\propto \exp\left( -\frac{1}{2} \left( \boldsymbol{x} _ {2} - \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{a} \right)^T \boldsymbol{\Lambda} _ {22} \left( \boldsymbol{x} _ {2} - \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{a} \right) \right) \end{align}

したがって, p\left(\boldsymbol{x}\ _ 2|\boldsymbol{x}\ _ 1 \right)は次の正規分布に従う.共分散行列が \boldsymbol{\Lambda} _ {22}の逆行列であることに注意.

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) &\sim \mathcal{N} \left( \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{a} , \boldsymbol{\Lambda} _ {22}^{-1} \right) \\ &= \mathcal{N} \left( \boldsymbol{\Lambda} _ {22}^{-1} (\boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) , \boldsymbol{\Lambda} _ {22}^{-1} \right) \\ &= \mathcal{N} \left( \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {22}^{-1}\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) , \boldsymbol{\Lambda} _ {22}^{-1} \right) \\ \end{align}

精度行列 \boldsymbol{\Lambda}を共分散行列 \boldsymbol{\Sigma}に戻したいのだが,そのためにブロック行列の逆行列を求める公式を使う.  \boldsymbol{M}=\left(
\boldsymbol{\Sigma}\ _ {22}-\boldsymbol{\Sigma}\ _ {21}\boldsymbol{\Sigma}\ _ {11}^{-1}\boldsymbol{\Sigma}\ _ {21}
\right)^{-1}
とおくと,

\begin{align} \Lambda= \begin{pmatrix} \boldsymbol{\Lambda} _ {11} & \boldsymbol{\Lambda} _ {12}\\ \boldsymbol{\Lambda} _ {21} & \boldsymbol{\Lambda} _ {22} \end{pmatrix} = \begin{pmatrix} \boldsymbol{\Sigma} _ {11} & \boldsymbol{\Sigma} _ {12}\\ \boldsymbol{\Sigma} _ {21} & \boldsymbol{\Sigma} _ {22} \end{pmatrix}^{-1} \end{align}

において,

\begin{align} \boldsymbol{\Lambda} _ {22} &= \boldsymbol{M} = \left( \boldsymbol{\Sigma} _ {22}-\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1}\boldsymbol{\Sigma} _ {21} \right)^{-1} \\ \boldsymbol{\Lambda} _ {21} &= -\boldsymbol{M}\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1} \end{align}

となる.二つを組み合わせれば

\begin{align} \boldsymbol{\Lambda} _ {22}^{-1} \boldsymbol{\Lambda} _ {21} &= -\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1} \end{align}

であるから,

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) &\sim \mathcal{N} \left( \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {22}^{-1}\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) , \boldsymbol{\Lambda} _ {22}^{-1} \right) \\ &= \mathcal{N} \left( \boldsymbol{\mu} _ {2} -\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) , \boldsymbol{\Sigma} _ {22}-\boldsymbol{\Sigma} _ {21}\boldsymbol{\Sigma} _ {11}^{-1}\boldsymbol{\Sigma} _ {21} \right) \\ \end{align}

に辿り付く.これがゴールである.完.

誤植?

19年2/4現在において,参考資料の先行公開版では微妙な誤りが見られた.少し混乱してしまったので,一応まとめておく.(私の勘違いや計算ミスだったらご指摘ください)

19/03/12追記:書籍版では一部修正されていました.

平方完成の直前,式(2.57)の6行目(γ2版)式(2.56)の5行目(書籍版)

 \expの中身にだけ注目する.

\begin{align} \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} -2\boldsymbol{x} _ {2}^T \left( \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) \right) \end{align}

のはずであるが,これが

\begin{align} \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} -2 \left( \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) \right)\boldsymbol{x} _ {2} \end{align}

となっていた.

19/03/12追記:書籍版では以下のようになっているが,誤り.全ての項が二次形式になるはずであるが,後半が(二次形式)ベクトルの形になっている.

19/03/12追記2:これについてはサポートページの正誤表に記載されています.

\begin{align} \boldsymbol{x} _ {2}^T \boldsymbol{\Lambda} _ {22} \boldsymbol{x} _ {2} -2\boldsymbol{x} _ {2}^T \left( \boldsymbol{\Lambda} _ {22} \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) \right)\boldsymbol{x} _ {2} \end{align}

式(2.58)(γ2版)式(2.57)(書籍版)以降,精度行列で表された正規分布

19/03/12追記:書籍版では修正されていました.

正しくは

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) &\sim \mathcal{N} \left( \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {22}^{-1}\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) , \boldsymbol{\Lambda} _ {22}^{-1} \right) \ \end{align}

だと思うが,

\begin{align} p\left(\boldsymbol{x} _ 2|\boldsymbol{x} _ 1 \right) &\sim \mathcal{N} \left( \boldsymbol{\mu} _ {2} +\boldsymbol{\Lambda} _ {22}^{-1}\boldsymbol{\Lambda} _ {21} \left(\boldsymbol{x} _ {1}-\boldsymbol{\mu} _ {1}\right) , \boldsymbol{\Lambda} _ {22} \right) \ \end{align}

となっていた.共分散行列は \boldsymbol{\Lambda}\ _ {22}の逆行列 \boldsymbol{\Lambda}\ _ {22}^{-1}になるはずである.

混合モデルの導入

 \displaystyle{}

はじめに

須山ベイズ本を参考に,混合モデルの導入部分を勉強したノート.離散確率分布の復習も含む.隠れマルコフモデルについて学びたかったのだが,そのために混合モデルについて押さえようと思ったらそれなりのボリュームになってしまったので一旦記事としてまとめておく.

あくまでも勉強ノートなので,教科書ありきの文章である.これを解説として読んでもらうのではなく,須山本を読む上での補足資料として参考になれば幸いである.単なる写経ではなく,式変形は細かくやっていたりするので,そのあたり誰かの役に立つかもしれない.

参考文献

予備知識

離散確率分布:ベルヌーイ分布とカテゴリ分布

混合モデルでは,データを構成するこの点たちは分布1から発生しており,他の点たちは分布2から発生しており・・・と,複数の異なる分布によってデータが構成されていると考える.データがどの分布から発生したのか,ということを考えるのはいわゆるクラスタリング問題であり,このとき有用になるのがベルヌーイ分布とカテゴリ分布である.これらを予備知識として復習しておこう. ベルヌーイ分布は2つのクラスタのみを考えるケースであり,ベルヌーイ分布を K次元に拡張したものがカテゴリ分布である.すなわちカテゴリ分布は K個のクラスタのうち,どのクラスタに所属するのかを示す確率分布になっている. 須山本ではカテゴリ分布と表記されているが,一般にはカテゴリカル分布と呼ばれるようである.マルチヌーイ分布という言い方もあるらしい.

クラスタ数2の場合:ベルヌーイ分布

まずベルヌーイ分布であるが, 0 1のどちらかを取る二値変数 xの生成分布である.これは 0から 1までの連続値を取るパラメータ \mu\in(0,1)によって形が決まり,次のように書く.

\begin{align} {\rm Bern}(x|\mu)=\mu^x(1-\mu)^{1-x} \end{align}

 x=0のときは (1-\mu)^{1}になるし, x=1のときは \mu^1になるから,これはつまり \muの値は x=1が得られる確率そのものだ,ということである.冗長だが,折角なので具体的な数字を入れてどうなるか見ておこう.  \mu=0.5の場合を考えてみると, \begin{align} {\rm Bern}(x|0.5)=0.5^x\times0.5^{1-x} \end{align} となって, x=0となる確率は  {\rm Bern}(0|0.5)=0.5^0\times0.5^{1}=0.5 となる.逆に x=1となる確率も \begin{align} {\rm Bern}(1|0.5)=0.5^1\times0.5^{0}=0.5 \end{align} となるから,この場合は 0 1は半々の確率で得られるということになる.

クラスタ数Kの場合:カテゴリ分布

ベルヌーイ分布ではクラスタが2つの場合について考えたが, K個の場合はどう表せるだろうか.

ここで変数 \boldsymbol{s}を考える(これは後々見るように,混合分布の文脈では直接観測に引っかからない潜在変数となる).

このベクトルの要素は s_k\in\{0,1\},すなわち 0 1かどちらかの値を取る離散的な変数であり,かつ \Sigma_{k=1}^Ks_k=1を満たす.つまり, k番目が 1であれば,他の要素は全て 0ということになる.このようにしてクラスタへの所属を表すことを 1 of  K表現と呼ぶ.

ベルヌーイ分布では1つのパラメータ \mu\in(0,1)によって分布が決まっていた.カテゴリ分布では,これを K個に拡張して,混合比率 \boldsymbol{\pi}=\left(\pi_1,\cdots,\pi_K\right)によって分布の形を決める. \pi_k\in(0,1)であり, \Sigma_{k=1}^K\pi_k=1を満たす. s_kが離散変数であるのに対して, \pi_kは連続値を取ることに注意しておこう.

以上で準備が整った.カテゴリ分布は,次のように書ける.

\begin{aligned} {\rm Cat}(\boldsymbol{s}|\boldsymbol{\pi})=\Pi _ {k=1}^{K} \pi _ {k}^{s _ {k}} \end{aligned}

クラスタ数 K=2の場合について考えてみよう.このとき潜在変数は  \boldsymbol{s}=\{s_1, s_2\}=\{s_1, 1-s_1\} となるし,混合比率は  \boldsymbol{\pi}=\{\pi_1, \pi_2\}=\{\pi_1, 1-\pi_1\} である.したがって

\begin{align} {\rm Cat}(\boldsymbol{s}| \boldsymbol{\pi}) &= \Pi _ {k=1}^K \pi_k^{s _ k} \\ &= \Pi _ {k=1}^2 \pi_k^{s _ k} \\ &= \pi _ 1^{s _ 1}\cdot\pi_2^{s _ 2} \\ &= \pi _ 1^{s _ 1}\cdot(1-\pi _ 1)^{1-s _ 1} \end{align}

となって,ベルヌーイ分布に一致することが確かめられた. 2つのクラスタのうち,1番目のクラスタが選ばれる確率を計算したいなら, \boldsymbol{s}=\{s_{1}, s_{2}\}=\{1, 0\}だから,

\begin{align} {\rm Cat}(\boldsymbol{s}|\boldsymbol{\pi}) &= \pi _ 1^{s _ 1}\cdot(1-\pi _ 1)^{1-s _ 1} \\ &=\pi _ 1^{1}\cdot(1-\pi _ 1)^{0} \\ &=\pi _ 1 \end{align}

となって \pi_1が残る. s_1 s_2が半々の確率で選ばれる場合, \pi_1=0.5である.ベルヌーイ分布と同様に,結局 \boldsymbol{\pi}によって \boldsymbol{s}の各要素が出現する確率が決まることが確認できた.

連続確率分布:ベータ分布とディリクレ分布

カテゴリ分布を決めるのは混合比率 \pi_k\in(0,1)であった.次はこの \piを得る確率について考えたい. \pi 0から 1までの値を取る連続変数であり,このような変数を生成する確率分布がベータ分布とディリクレ分布である.今回は知識として眺めるだけに留めておく.「そういうもん」と思って欲しい.

ベータ分布

ベルヌーイ分布で出て来た \mu\in(0,1)の分布を与えてくれるのがベータ分布である.パラメータ a bを使って次のように書く.  {\rm Beta}(\mu|a,b)=C_B(a,b)\mu^{a-1}(1-\mu)^{b-1}  C_B(a,b) a bの値に応じて決まる規格化係数で,ガンマ関数によって構成される.

ディリクレ分布

カテゴリ分布で出て来た混合比率 \pi_k\in(0,1)の分布を与えてくれるのがディリクレ分布である.パラメータ \boldsymbol{\alpha}を使って次のように書く.

\begin{align} {\rm Dir}(\boldsymbol{\pi}|\boldsymbol{\alpha}) =C _ D(\boldsymbol{\alpha})\Pi _ {k=1}^K \pi_k^{\alpha _ k-1} \end{align}

規格化係数 C_D(\boldsymbol{\alpha})はやはりガンマ関数によって構成される.

混合モデル

混合モデルでは,データが複数の異なる分布によって生成していると考える. データを生成する確率分布の個数(クラスタ数)を Kとして, N個のデータ \boldsymbol{X}=\{\boldsymbol{x_1},\cdots, \boldsymbol{x_N} \}の生成過程を考えよう. ちなみに, \boldsymbol{X}の1点である \boldsymbol{x_n}は多次元空間上の1つの点であり,これ自体がベクトルである.たとえば,3次元空間上にばらまかれている点の集団がデータ \boldsymbol{X}であれば,その一つ一つの点は3次元ベクトルであって \boldsymbol{x _ n}=(x _ {n,x},x _ {n,y},x _ {n,z})^Tのように書ける. 以下では混合モデルを考える上で道具となる量について見ていく.

クラスタの混合比率 \pi

データのうち,それぞれの分布(クラスタ)がどの程度含まれているか, 0から 1までの値で表すことにする.これをベクトルとしてまとめたものが混合比率(mixture weight) \piである.  \boldsymbol{\pi}=\left(\pi _ 1,\cdots,\pi _ K\right)^T ここで \pi_K\in(0,1)である.また全体で 1となるように \Sigma_{k=1}^K\pi_k=1とする.さて,予備知識のところで見たように, \piは"混合比率"なのだから, 0から 1までの連続値を取ることに注意しておこう.これに対して,この後出てくる潜在変数 sはどのクラスタに属するかを 0 1で表す離散的な量である.基本的なことではあるが確認のため記しておくと,このような離散的な集合は波括弧を用いて \{\cdot\}で表す.何でこんなことをくどくど書くのかというと,単に私の脳が混乱したからにすぎない.  \boldsymbol{\pi}も確率変数であって,事前分布 p(\boldsymbol{\pi})に従って生成されるとする.

モデルパラメータ \theta

たとえば各分布が正規分布の場合,平均値と分散がモデルパラメータとなる.分布ごとにパラメータは異なるから, k番目の分布のパラメータセットをベクトル \boldsymbol{\theta}_kで表す.正規分布の場合であれば, \boldsymbol{\theta} _ k=\mathcal{N}\left(\mu _ k, \sigma _ k\right)となる. \boldsymbol{\theta} _ kは事前分布 p(\boldsymbol{\theta} _ k)に従って生成されるとする.

潜在変数 s_n

観測データ \boldsymbol{X}=\{\boldsymbol{x_1},\cdots, \boldsymbol{x_N} \}のある点 \boldsymbol{x_n} K個の分布のうちどれによって生成するかという対応関係を決める変数が s_nであり,潜在変数と呼ばれる. s_n=\{s_{n,1}, \cdots, s_{n,K}\} であり,予備知識で見たとおり 1 of  K表現によってどのクラスタに所属しているかを表す.これが N個分あるのだから,その集合を \boldsymbol{S}=\{\boldsymbol{s} _ 1, \cdots, \boldsymbol{s} _ N\}と表す. また,データがどのクラスタに所属しているか(どのクラスタから生成したか)?というのは,そのデータがどのモデルパラメータに対応しているのかを考えることと同じである.

グラフィカルモデル

ここで考えているデータの生成過程は下に示すグラフィカルモデルで表すことができる.

undefined.jpg

 \boldsymbol{x} _ n \boldsymbol{\Theta} \boldsymbol{s} _ nに依存する.注意しておくべきなのは, \boldsymbol{\Theta} \boldsymbol{s} _ nの間には依存関係が無いということ.各クラスタの分布がどのような形になっているかということ( \boldsymbol{\Theta})と,どの分布が選ばれるかということ( \boldsymbol{s} _ n)は無関係である.

余談であるが,図のように N個の量を一つのノードで表し,それらを四角で囲ってから個数を Nと記す省略表記方法はプレート表現(plate notation)と言われる.また, \boldsymbol{x}_nが観測値であることはノードを塗りつぶすことによって示す.このあたりは須山本p.24や,PRMLの8章で解説がなされている.

混合分布の構成

 \boldsymbol{x} _ nがどのクラスタ(モデルパラメータ)に属するかは潜在変数 s _ nで決まるのであった. \boldsymbol{x} _ n s _ nとモデルパラメータ全体の集合 \boldsymbol{\Theta}=\{\boldsymbol{\theta} _ 1,\cdots\boldsymbol{\theta} _ k\}に依存して決まるので, \boldsymbol{x} _ nの生成確率は p(\boldsymbol{x} _ n|s_n,\boldsymbol{\Theta})と書ける. 予備知識として見たように, K個のクラスタについてのカテゴリ分布は次のように書けるのだった.

\begin{align} {\rm Cat}(\boldsymbol{s}|\boldsymbol{\pi})=\Pi _ {k=1}^K \pi _ k^{s _ k} \end{align}

これを参考にして, \boldsymbol{x}_nの生成確率は次のように書ける.

\begin{align} p(\boldsymbol{x} _ n|s _ n,\boldsymbol{\Theta})=\Pi _ {k=1}^K p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ k)^{s _ {n,k}} \end{align}

となる.くどいようだが,やはりひとまず具体的に見ていこう. K=2の場合を考えると, \boldsymbol{s}_n=\{s_{n,1}, s_{n,2}\}=\{s_{n,1}, 1-s_{n,1}\}なんだから,

\begin{align} p(\boldsymbol{x} _ n|s _ n,\boldsymbol{\Theta}) &=\Pi _ {k=1}^2 p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ k)^{s _ {n,k}} \\ &=p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 1)^{s _ {n,1}}\cdot p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 2)^{1-s _ {n,1}} \end{align}

ということである. \boldsymbol{x}_nが1番目のクラスタから生成されたとすると, \boldsymbol{s}_n=\{s_{n,1}, s_{n,2}\}=\{1, 0\}だから,

\begin{align} p(\boldsymbol{x} _ n|s _ n,\boldsymbol{\Theta}) &=\Pi _ {k=1}^2 p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ k)^{s _ {n,k}} \\ &=p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 1)^{s _ {n,1}}\cdot p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 2)^{1-s _ {n,1}} \\ &=p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 1)^{1}\cdot p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 2)^{0} \\ &=p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 1) \end{align}

となって p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ 1)が残る.  p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ k)はあるモデルパラメータを選んだときに \boldsymbol{x}\ _ nが得られる確率である.モデルが正規分布だとすればモデルパラメータは \boldsymbol{\theta} _ k=\{\boldsymbol{\mu} _ k, \boldsymbol{\Sigma
} _ k\}であり, k=1,\cdots,Kに対して,

\begin{align} p(\boldsymbol{x} _ n|\boldsymbol{\theta} _ k) = \mathcal{N}(\boldsymbol{x} _ n|\boldsymbol{\mu} _ k, \boldsymbol{\Sigma} _ k ) \end{align}

となる.

さて,潜在変数 s 1であるか 0であるかの確率は混合比率 \piで決まり,

\begin{align} {\rm Cat}(\boldsymbol{s}|\boldsymbol{\pi})=\Pi _ {k=1}^K \pi _ k^{s _ k} \end{align}

だった.そして \pi _ 1 0.5か, 0.1かといったように, \boldsymbol{\pi}の各要素がどのような値を取るかというのは事前分布 p(\boldsymbol{\pi})によって決まると考えるのだった.この事前分布としてディリクレ分布を使う.すなわち  p(\boldsymbol{\pi})={\rm Dir}(\boldsymbol{\pi}|\boldsymbol{\alpha}) である.

同時分布

以上で下準備が終わった.それでは N個のデータ \boldsymbol{X}=\{\boldsymbol{x\ _ 1},\cdots, \boldsymbol{x\ _ N} \},潜在変数 \boldsymbol{S}=\{\boldsymbol{s} _ 1, \cdots, \boldsymbol{s} _ N\},モデルパラメータ \boldsymbol{\Theta}=\{\boldsymbol{\theta} _ 1,\cdots\boldsymbol{\theta} _ k\},混合比率 \boldsymbol{\pi}=\left(\pi\ _ 1,\cdots,\pi\ _ K\right)^Tの同時分布 p(\boldsymbol{X},\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi})を考える. 乗法定理によって条件付き分布に書き換えていくと,

\begin{align} p(\boldsymbol{X},\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi}) &= p(\boldsymbol{X},\boldsymbol{S},\boldsymbol{\Theta}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ &= p(\boldsymbol{X},\boldsymbol{\Theta}|\boldsymbol{S},\boldsymbol{\pi})p(\boldsymbol{S}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ &= p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}, \boldsymbol{\pi})p(\boldsymbol{\Theta}|\boldsymbol{S},\boldsymbol{\pi})p(\boldsymbol{S}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ \end{align}

ここでグラフィカルモデルを思い出す.

undefined.jpg

図から明らかなように, \boldsymbol{\Theta} \boldsymbol{s} _ nの間には依存関係が無い. \boldsymbol{\Theta} \boldsymbol{s} _ nにも \boldsymbol{\pi}にも無関係なのだから,上の同時分布の式で出て来た p(\boldsymbol{\Theta}|\boldsymbol{S},\boldsymbol{\pi})は,単なる p(\boldsymbol{\Theta})に等しい.つまり,

\begin{align} p(\boldsymbol{\Theta}|\boldsymbol{S},\boldsymbol{\pi}) = p(\boldsymbol{\Theta}) \end{align}

である.もう一つ,図から明らかなように, \boldsymbol{x} _ n \boldsymbol{\Theta} \boldsymbol{s} _ nに依存するけれど \boldsymbol{\pi}には直接依存しない.したがって,

\begin{align} p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}, \boldsymbol{\pi}) = p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}) \end{align}

となる.というわけで同時分布の式は,次のように書き換えられる.

\begin{align} p(\boldsymbol{X},\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi}) &= p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}, \boldsymbol{\pi})p(\boldsymbol{\Theta}|\boldsymbol{S},\boldsymbol{\pi})p(\boldsymbol{S}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ &= p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta})p(\boldsymbol{\Theta})p(\boldsymbol{S}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ \end{align}

ここで, p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta})

\begin{align} p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}) &= p(\boldsymbol{x} _ 1,\cdots, \boldsymbol{x} _ N|\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N,\boldsymbol{\Theta}) \\ &=p(\boldsymbol{x} _ 1|\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N,\boldsymbol{\Theta}) \cdots p(\boldsymbol{x} _ N|\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N,\boldsymbol{\Theta}) \\ \end{align}

であって, \boldsymbol{x} _ n \boldsymbol{S}のうち \boldsymbol{s} _ nにしか依存しないということを考慮すると,

\begin{align} p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}) &= p(\boldsymbol{x} _ 1,\cdots, \boldsymbol{x} _ N|\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N,\boldsymbol{\Theta}) \\ &=p(\boldsymbol{x} _ 1|\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N,\boldsymbol{\Theta}) \cdots p(\boldsymbol{x} _ N|\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N,\boldsymbol{\Theta}) \\ &=p(\boldsymbol{x} _ 1|\boldsymbol{s} _ 1,\boldsymbol{\Theta}) \cdots p(\boldsymbol{x} _ N|\boldsymbol{s} _ N,\boldsymbol{\Theta}) \\ &= \Pi _ {n=1}^N p(\boldsymbol{x} _ n|\boldsymbol{s} _ n,\boldsymbol{\Theta}) \end{align}

 p(\boldsymbol{S}|\boldsymbol{\pi})についても同様.

\begin{align} p(\boldsymbol{S}|\boldsymbol{\pi}) &= p(\boldsymbol{s} _ 1,\cdots, \boldsymbol{s} _ N|\boldsymbol{\pi}) \\ &=p(\boldsymbol{s} _ 1|\boldsymbol{\pi}) \cdots p(\boldsymbol{s} _ N|\boldsymbol{\pi})\\ &= \Pi _ {n=1}^N p(\boldsymbol{s} _ n|\boldsymbol{\pi}) \end{align}

次に p(\boldsymbol{\Theta})であるが,

\begin{align} p(\boldsymbol{\Theta}) &= p(\boldsymbol{\theta} _ 1,\cdots, \boldsymbol{\theta} _ K) \\ &=p(\boldsymbol{\theta} _ 1) \cdots p(\boldsymbol{\theta} _ k)\\ &= \Pi _ {k=1}^K p(\boldsymbol{\theta} _ k) \end{align}

したがって同時分布の式は

\begin{align} p(\boldsymbol{X},\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi}) &= p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta}, \boldsymbol{\pi})p(\boldsymbol{\Theta}|\boldsymbol{S},\boldsymbol{\pi})p(\boldsymbol{S}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ &= p(\boldsymbol{X}|\boldsymbol{S},\boldsymbol{\Theta})p(\boldsymbol{\Theta})p(\boldsymbol{S}|\boldsymbol{\pi})p(\boldsymbol{\pi}) \\ &= \left( \Pi _ {n=1}^N p(\boldsymbol{x} _ n|\boldsymbol{s} _ n,\boldsymbol{\Theta}) p(\boldsymbol{s} _ n|\boldsymbol{\pi}) \right) \left( \Pi _ {k=1}^K p(\boldsymbol{\theta} _ k) \right) p(\boldsymbol{\pi}) \end{align}

こうして須山本の式(4.5)に辿り着きました.めでたし.

事後分布

データが所属するクラスタの推定を行いたい場合,事後分布 p(\boldsymbol{S}|\boldsymbol{X})を求めれば良い.

\begin{align} p(\boldsymbol{S}|\boldsymbol{X}) &=\int \int p(\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi}|\boldsymbol{X}){\rm d}\boldsymbol{\Theta}{\rm d}\boldsymbol{\pi}\\ p(\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi}|\boldsymbol{X}) &= \frac{p(\boldsymbol{X},\boldsymbol{S},\boldsymbol{\Theta},\boldsymbol{\pi})} {p(\boldsymbol{X})} \end{align}

であるが,この計算は一般に困難とされている.その解決法として,なんらかの近似アルゴリズムによって推論を行おう,ということでMCMCや変分推論に話が繋がっていく.

【強化学習】Bellman方程式の導出

はじめに

Qrunchに投稿した記事に少し追記したものです.というのも,現状Qrunchだと外部からの検索で上手く引っかからないよう(つまりGoogle検索で出てこない).QrunchはQrunchコミュニティ内での気軽な技術系情報共有を目的にしているところがコンセプトっぽいので,あえて閉鎖的なコンテンツにしているということなのかもしれません.個人の勉強メモとして引き続きQrunchも活用していきたいと思っていますが,多少なりとも誰かの役に立ちそうな内容のものは積極的に外部に公開しても良いのかなと考えて,はてなも活用していこうと思っています.(Qrunchには「クロス投稿」機能があるので,それを使うのが良さそう)

モチベーション

強化学習のお勉強のため,Pythonで実装しようぜ系の本をしばらく読んでいたのですが,肝心のBellman方程式はわりとサラッと導出されていて,自分で手を動かして式を追ってみないとなんだかよくわからない.そうするにあたって瀧先生の深層学習本の強化学習の章における説明が個人的に分かりやすかったので,これを参考にBellman方程式導出の過程をまとめておく.この記事の目的はあくまでもBellman方程式を導出するまでのひとつひとつの式変形のステップを省略せずにまとめておくことなので,強化学習そのものについては次節で挙げている書籍を参考にされてください.

Bellman方程式ってなに

Bellman方程式とは,価値関数を状態 sと次の状態 s'についての漸化式として表したものである.再帰的に価値関数を表すことの必要性は書籍「Pythonで学ぶ強化学習」に易しい解説があり,分かりやすいと思う.

Bellman方程式の導出の流れ

量の書き換えがいくつか出て来て最初は混乱したのですが,やっていることは単純で,以下の二点にまとめられる.

  1. 報酬の再帰性を使って,価値関数を漸化式の形に書く.

  2. マルコフ性によって依存関係を整える.

参考書籍

メイン

サブ

計算ルール(確率,期待値の基本)

Bellman方程式の導出で必要になる計算ルールを予備知識としてまとめておく. これらに加えて,強化学習の問題設定がマルコフ決定過程であることを利用して話が進んでいきます. 確率の基礎にまつわる話は個人的に「プログラミングのための確率統計」が直感的で分かりやすいので,脳が混乱したときには読み返している.

周辺化

基本

 \displaystyle{
P(A)=\sum_{B}P(A,B)
}

ですよね.同時分布をある確率変数について和を取れば,括弧からそいつが消えます.得られた P(A)を周辺分布という.

Bellman方程式の導出では,たとえば遷移確率を周辺化によって次のように式変形する. 強化学習では基本的に条件付き分布ばかりを扱いますが,上で見た基本式のおしりに条件がくっついてるだけ.

 \displaystyle{
P^{\pi}(s'|s)=\sum_{a}P(s',a|s)
}

同時確率と条件付き確率の関係(乗法定理)

基本

 P(A,B)=P(A|B)P(B)

これを使って先ほどの周辺化の例を,  \displaystyle{
P^{\pi}(s'|s)=\sum _ aP(s',a|s)=\sum _ aP(s'|a,s)P(a|s)
}

というように変形していったりする.

条件付き期待値

基本1

条件付き期待値の定義はこう.


\displaystyle{
E[Y|a]= \sum _ {b} bP(Y=b|a)
}

基本2


\displaystyle{
\sum_aE[Y|a] P(a)=E[Y]
}

こうなることを一応確認しておく.

\begin{align} \sum_aE[Y|a]P(a) &= \sum_a\sum_bbP(b|a)P(a) \\ &= \sum_a\sum_bbP(a,b) \\ &= \sum_bb\left(\sum_aP(a,b)\right) \\ &= \sum_bbP(b) \\ &= E[Y=b] \end{align}


\displaystyle{
V ^ \pi (s)=E _ {p,\pi}[R(t)|s=s(t)]=\sum _ a E[R(t)|s,a]P(a|s)
}

本題

登場人物(強化学習に出てくる基本概念)

遷移確率

状態 sで行動 aを取ったときに状態 s'に遷移する確率. 
P(s'|s,a)=P(s'=s(t+1)|s=s(t),a=a(t))

取り得る全ての行動について足し合わせると(計算ルールのところで確認した周辺化),状態 sから s'への遷移確率になる.

  P^{\pi}(s'|s)=\sum_{a}P(s'|s,a) P(a|s)=\sum_a\pi(a|s)P(s'|s,a)

ここで,状態 sにおいて行動 aを選択する確率は方策と呼ばれ,特別に \piと書く.

 \pi(a|s)=P(a|s)

利得(累積報酬)

 \displaystyle{
R(t)=r(t+1)+\gamma r(t+2)+\gamma^2 r(t+3)+\cdots=\sum_{k=0}^\infty \gamma^k r(t+k+1)
}

 \gammaは割引率.累積報酬とか割引報酬和,割引累積報酬なんて呼ぶのが意味が分かりやすくて好きだが,単に報酬や利得と呼ばれることも多いっぽい.これは第1項と残りを分けることで再帰的に書くことができる.

 \displaystyle{
\begin{align}
R(t)&=r(t+1)+\gamma \left(r(t+2)+\gamma r(t+3)+\cdots\right)\\
&=r(t+1)+\sum_{k=0}^\infty \gamma^k r(t+k+2) \\
&=r(t+1)+\gamma R(t+1)
\end{align}
}

ちなみに r=r(t+1)は,行動 aによって,状態 s=s(t)から状態 s'=s(t+1)になったときに得る即時報酬.貰うタイミングは時刻 t+1

即時報酬の期待値

大文字の Rで書いているが,上の累積報酬ではなく時刻 t+1の即時報酬についての話である.

  R(s,a,s')=E_P [r(t+1)|s',s,a]=\sum_{rr}P(r|s',s,a)

上で書いているように r=r(t+1)なので,右辺の意味は以下のとおり.  \sum_{r(t+1)}r(t+1)P(r(t+1)|s',s,a)

状態 sでの報酬の期待値は,上記の量を行動 aと次の状態 s'について期待値を取る.要するに,行動や次の状態にかかわらず,「ただ状態 sにいるだけで」これくらいの報酬が得られるだろうと期待できるということ.ここで重みとして用いる確率は, aを取る確率と, s'となる確率なのだから,すなわち方策と遷移確率ですね.  R^\pi (s)=\sum_{a,s'}\pi(a|s)P(s'|s,a)R(s,a,s')

価値関数

状態価値関数 V

ある状態 sであったとき,得られるであろう累積報酬 R(t)の期待値を状態価値関数と呼ぶ.これによって各状態の良し悪し(価値)が評価できる.  V^\pi(s)=E_{P, \pi}[R(t)|s] ここでの添字 P \piの意味は,さきほど出て来た即時報酬の期待値と同様に,遷移確率 Pと方策 \piで期待値を取っているよということ.

行動価値関数 Q

状態 sにおいて行動 aを取ったときの累積報酬 R(t)の期待値.  Q^\pi(s,a)=E_P[R(t)|s,a]

 V Qの関係

計算ルールの条件付き期待値のところで見たように,  V^\pi(s)=E_{P, \pi}[R(t)|s]=\sum_aE_P[R(t)|s,a]P(a|s)=\sum_a\pi(a|s)E_P[R(t)|s,a]

と書き換えられる.最右辺の条件付き期待値はまさに Qだから,  V^\pi(s)=\sum_a\pi(a|s)Q^\pi(s,a)

ようやくBellman方程式へ

状態価値関数 VについてのBellman方程式

さきほどの「 V Qの関係」で出て来たように,

 V^\pi(s)=E_{P, \pi}[R(t)|s] =\sum_a\pi(a|s)E_P[R(t)|s,a]

予備知識の計算ルールである条件付き期待値の性質を使って,条件に s'を追加する.

 E_P[R(t)|s,a]=\sum_{s'} E_P[R(t)|s',s,a]P(s'|s,a)

するとVは

 \displaystyle{
\begin{align}
V^\pi(s)
&=E_{P, \pi} [R(t)|s] \\
&=\sum_a\pi(a|s)E_P[R(t)|s,a] \\
&=\sum_{a,s'}\pi(a|s)P(s'|s,a)E_P[R(t)|s',s,a]
\end{align}
}

となりますね.

累積報酬を再帰的に書くと  R(t)=r(t+1)+\gamma R(t+1)

だったから,これを上の式に入れてやる.そうすると,

 \displaystyle{
\begin{align}
V^\pi(s)&=\sum_{a,s'}\pi(a|s)P(s'|s,a)E_P[R(t)|s',s,a] \\
&=\sum_{a,s'}\pi(a|s)P(s'|s,a)E_P[r(t+1)+\gamma R(t+1)|s',s,a] \\
&=\sum_{a,s'}\pi(a|s)P(s'|s,a)\left(E_P[r(t+1)|s',s,a]+E_P[\gamma R(t+1)|s',s,a]\right)
\end{align}
}

となって,即時報酬 r(t+1)の期待値が出て来た.これは上で見たように

  R(s,a,s')=E_P [r(t+1)|s',s,a]

と表していたので,これを使ってさらに書き換えると

 \displaystyle{
\begin{align}
V^\pi(s)&=\sum_{a,s'}\pi(a|s)P(s'|s,a)\left(E_P[r(t+1)|s',s,a]+E_P[\gamma R(t+1)|s',s,a]\right) \\
&=\sum_{a,s'}\pi(a|s)P(s'|s,a)\left(R(s,a,s')+E_P[\gamma R(t+1)|s',s,a]\right)
\end{align}
}

ここで出て来た R^\pi (s)

 R^\pi (s)=\sum_{a,s'}\pi(a|s)P(s'|s,a)R(s,a,s')

なのでした.したがって

 \displaystyle{
V^\pi(s)=R^\pi (s)+\sum_{a,s'}\pi(a|s)P(s'|s,a)E_P[\gamma R(t+1)|s',s,a]
}

この第2項目には遷移確率が含まれているから,

 \displaystyle{
\sum_{s'}\left(\sum_{a}\pi(a|s)P(s'|s,a)\right)E_P[\gamma R(t+1)|s',s,a]=\gamma \sum_{s'}P^\pi(s'|s)E_P[R(t+1)|s',s,a]
}

さてここで, R(t+1)は状態 s'=s(t+1)から s''=s(t+2)となった時点以降の累積報酬である.これは状態 sと行動 aには依存しない(マルコフ性). そういうわけで上の式の中の条件付き期待値は,

 E_P[R(t+1)|s',s,a]=E_P[R(t+1)|s']=V^\pi (s')

となり,これは結局,次の状態 s'での状態価値関数 V(s')であった. したがって最終的に,

 V^\pi(s)=R^\pi (s)+\gamma \sum_{s'}P^\pi(s'|s)V^\pi (s')

となって,これが状態価値関数についてのBellman方程式です.

行動価値関数 QについてのBellman方程式

行動価値関数は

 Q^\pi(s,a)=E_P[R(t)|s,a]

なのでした.これを Vに対してやったのと同様に,計算ルールである条件付き期待値の性質を使って s'について和を取ります. それから R(t)=r(t+1)+\gamma R(t+1)を代入して二つの項に分けます. そうすると第1項には R(s,a,s')が現れます. 第2項はマルコフ性によって条件付き期待値から s aが除かれて,その結果 V^\pi(s')が出て来る.

 \displaystyle{
\begin{align}
Q^\pi(s,a)&=E_P[R(t)|s,a] \\
&=\sum_{s'}P(s'|s,a)E_P[R(t)|s,a,s'] \\
&=\sum_{s'}P(s'|s,a)E_P[r(t+1)+\gamma R(t+1)|s,a,s'] \\
&=\sum_{s'}P(s'|s,a)E_P[r(t+1)|s,a,s']+\sum_{s'}P(s'|s,a)E_P[\gamma R(t+1)|s,a,s'] \\
&=\sum_{s'}P(s'|s,a)R(s,a,s')+\gamma \sum_{s'}P(s'|s,a)E_P[R(t+1)|s,a,s'] \\
&=\sum_{s'}P(s'|s,a)R(s,a,s')+\gamma \sum_{s'}P(s'|s,a)E_P[R(t+1)|s'] \\
&=\sum_{s'}P(s'|s,a)R(s,a,s')+\gamma \sum_{s'}P(s'|s,a)V^\pi(s')
\end{align}
}

 V Qの関係」で見たように

 V^\pi(s)=\sum_a\pi(a|s)Q^\pi(s,a)

だったから,

 V^\pi(s')=\sum_{a'}\pi(a'|s')Q^\pi(s',a')

である.これを使って書き換えれば,

 \displaystyle{
\begin{align}
Q^\pi(s,a)&=\sum_{s'}P(s'|s,a)R(s,a,s')+\gamma \sum_{s'}P(s'|s,a)\sum_{a'}\pi(a'|s')Q^\pi(s',a') \\
&=\sum_{s'}P(s'|s,a)R(s,a,s')+\gamma \sum_{s'}\sum_{a'}\pi(a'|s')P(s'|s,a)Q^\pi(s',a')
\end{align}
}

となる.これが行動価値関数 QについてのBellman方程式.

Sutton本との比較

Sutton本と瀧深層学習本はそれぞれ量の表記が異なるので,同じ方程式であるのか一目見るだけでは私は脳のメモリが足りずよく分からない.というわけで互いに同一の結果を得ていることを確認しておく.Sutton本だとp.76の式(3.10)で VについてのBellman方程式が示されている.

 \displaystyle{
V^\pi(s)=\sum_{a}\pi(s,a)\sum_{s'}P_{ss'}^{a}\left[R_{ss'}^{a}+\gamma V^\pi (s')\right]
}

まず \pi(s,a)であるが,p.56にて説明がなされており,これはこれまで見てきた \pi(a|s)と同じものである.Sutton本では, \pi s aの関数であることを明示した表記を取っていると私は解釈した. 次に P_{ss'}^{a} R_{ss'}^{a}はp.72の式(3.6)と(3.7)で説明がされているように,

 P_{ss'}^{a}=P(s'|s,a)

 \displaystyle{
R_{ss'}^{a}=E[r_{t+1}|s,a,s']
}

である.ここで r_{t+1}と書かれているのは, r_{t+1}=r(t+1)ということである. 以上を踏まえて,瀧深層学習本での表記に合わせるようにSutton本のBellman方程式を書き換えてみる.

 \displaystyle{
\begin{align}
V^\pi(s)&=\sum_{a}\pi(s,a)\sum_{s'}P_{ss'}^{a}\left[R_{ss'}^{a}+\gamma V^\pi (s')\right] \\
&=\sum_{a}\pi(a|s)\sum_{s'}P(s'|s,a)\left(E_P[r(t+1)|s,a,s']+\gamma V^\pi (s')\right)
\end{align}
}

この時点で,あぁ,同じやんけ,と安心するのだが,折角なので最後まで書き換える.ここで以下の3つを思い出す.

 E_P [r(t+1)|s',s,a]=R(s,a,s')

 R^\pi (s)=\sum_{a,s'}\pi(a|s)P(s'|s,a)R(s,a,s')

  P^{\pi}(s'|s)=\sum_a\pi(a|s)P(s'|s,a)

これらを使って書き換えると,

 \displaystyle{
\begin{align}
V^\pi(s)=R^\pi (s)+\gamma \sum_{s'}P^\pi(s'|s)V^\pi (s')
\end{align}
}

となって,瀧深層学習本のBellman方程式に一致した.安心した.