2状態Gaussian HMMで圧力からON/OFFを推定する

はじめに

設備のセンサデータ（電流・圧力・温度など）を持っていて，こんなことをしたい場面はないでしょうか．

• 設備の調子を数値で追いたい． 期間どうしで，圧力や電流の水準がどれだけ変わったかを比べたい
• 稼働中と停止中を，別々の監視項目として評価したい． 稼働中の値は負荷時の性能を，停止中の値はベースラインを表します．意味が違うので，混同せずにそれぞれ評価したい．
• 稼働率や起動回数そのものを知りたい

どれをやるにも，その前に各時刻に設備が動いていたか止まっていたかの情報が必要です．ところが現実のログには，運転信号が残っていないことが珍しくありません．あるのはセンサ値の列だけ．だったらセンサ値そのものから稼働/停止を切り分けるしかない――というのがこの記事の出発点です．

なぜ切り分けをサボれないのか

全体の平均は，停止中の平均 $\mu_{\mathrm{OFF}}$ と稼働中の平均 $\mu_{\mathrm{ON}}$ を稼働率 $r$ で混ぜた加重平均です．

\bar{x} = (1 - r)\,\mu_{\mathrm{OFF}} + r\,\mu_{\mathrm{ON}}

$\mu$ がまったく動かなくても，$r$ が動けば $\bar{x}$ は動きます．停止中がほぼ0，稼働中が8 bar台のような信号では，稼働率が数ポイント動くだけで全体平均は大きく変わります．つまり稼働と停止を混ぜたままの平均は，値の指標ではなくほぼ稼働率の関数です．値そのものがまったく動いていなくても，大きく変化したように見えてしまいます（後半の実験で実例を示します）．

この記事でやること

• 稼働/停止のラベルを教師なしで推定します．使うのは2状態Gaussian HMMです
• パラメータの学習（EMアルゴリズム）と状態の推論（forward-backward）を，「なぜその式になるのか」から順に導出します
• 地下鉄車両の空気圧縮機の公開データ（MetroPT-3，約150万点）で，切り分けの効果を確認します

なぜHMMなのか

素直な発想は「ある値以上ならON」という閾値判定です．波形を眺めて線を1本引けば，たいていの区間は正しく分かれます．問題は，その線がいつまで正しいかです．

• 水準そのものが動きます． 経年変化や運転条件の違いでON時の値がずれても，固定した閾値は追従しません．期間ごとに引き直すなら，それは人手の作業になります
• 1点のノイズで判定が反転します． 各時刻を独立に見るので，値が一瞬だけ線をまたぐと状態も一瞬だけ切り替わります
• 判定の確からしさが残りません． 閾値を9 barとすると，10.5 barの点も9.01 barの点も同じ「ON」になります．後者は境界すれすれで怪しいのですが，0/1に潰れた後では区別がつきません

しかも，これらは例外を出さずに壊れがちです．判定がずれても計算は正常に終わるため，稼働率だけが静かにずれていきます．

隠れマルコフモデル（HMM）は，この3点にそれぞれ答えます．

• 各状態に対応する観測分布を正規分布としてデータから学習します．OFFとONの水準が変化すれば，推定される分布もそれに応じて変化します
• 遷移確率によって時間的なつながりをモデル化できるため，前後と整合しない一瞬の状態反転を起こりにくくできます
• 状態を0/1ではなく確率で持てます．「ONらしさ0.98」と「ONらしさ0.51」を区別できます

では，HMMとはどんなモデルで，どんなパラメータを持つのでしょうか．

モデルとパラメータ

手元にあるのは観測値の列$X_1, \dots, X_T$だけです．知りたいのは各時刻の状態$S_t \in \{\mathrm{OFF}, \mathrm{ON}\}$ですが，これは記録されていません（隠れ状態）．しかも「ONのとき圧力はいくつか」も分かっていません．状態も，状態を決める基準も，どちらも未知という状況です．

そこで次のモデルを仮定します．パラメータは以下の4種類を使用します．

記号	意味
$\pi_k = P(S_1 = k)$	初期状態：最初の時刻に状態 $k$ である確率
$A_{ij} = P(S_t = j \mid S_{t-1} = i)$	遷移確率：状態 $i$ の次に $j$ になる確率
$b_t(k) = p(X_t \mid S_t = k) = \mathcal{N}(X_t;\, \mu_k,\, \sigma_k^2)$	観測分布：状態 $k$ のときに $X_t$ が観測される確率密度．平均 $\mu_k$ と標準偏差 $\sigma_k$ がパラメータ

このモデルの下では，観測と状態が同時に起きる確率が1本の式で書けます．

p(X_{1:T},\, s_{1:T} \mid \theta) =
\underbrace{\pi_{s_1}}_{\text{最初の状態}}\;
\underbrace{\prod_{t=2}^{T} A_{s_{t-1} s_t}}_{\text{状態の遷移確率}}\;
\underbrace{\prod_{t=1}^{T} b_t(s_t)}_{\text{各時刻の観測尤度}}

以降の式はすべてここから出発します．$\theta = (\pi, A, \mu, \sigma)$ をまとめてパラメータと呼びます．やることは2段階です．

では，このパラメータをどうやって決めるのでしょうか．

学習 ― EMアルゴリズム

最尤推定と2つの課題

決め方は最尤推定です．手元の観測がもっとも起こりやすくなる $\theta$ を選びます．ただし状態列が見えないので，尤度はあり得る状態列すべての和になります．

p(X_{1:T} \mid \theta) = \sum_{s_{1:T}} p(X_{1:T},\, s_{1:T} \mid \theta)

しかしながらここで2つの課題に直面します．

課題1：和の項数が$2^T$個ある． 状態列はOFF/ONの2択が$T$個並んだものなので，$T = 150$万では到底列挙できません．しかもこれは「最良の状態列を選ぶ」数ではなく，$\theta$を1つ評価するたびに足し尽くす必要のある数です．

課題2：閉形式の更新式が得られない． 対数の中に和が残るため（$\log(a+b)$は分解できません），パラメータごとに項が分離しません．$\mu_k$で微分すると「状態の事後確率を重みにした加重平均」の形にはなりますが，その重み自体が$\theta$に依存します．$\mu$を求めるには状態確率が要り，状態確率を求めるには$\mu$が要る――解ではなく不動点方程式です．

状態が見えないのに，どうやって学習すればよいのでしょうか．ここで発想を変えます．仮に状態列が全部見えていれば，パラメータは「遷移を数える」「平均を取る」だけで求まります．ならば，見えない状態を確率で埋めて数えればよい．これがEMアルゴリズムの考え方です．では，なぜそれでよいと言えるのでしょうか．

なぜEMでよいのか ― ELBO

状態列の上の任意の分布$q(s)$に対して，対数尤度は恒等的に2つの項へ分解できます．

\log p(X \mid \theta) = \mathrm{ELBO}(q, \theta) + \mathrm{KL}\bigl(q \,\|\, P(s \mid X, \theta)\bigr)

ここでELBO（Evidence Lower BOund，証拠下界）は次の量です．

\mathrm{ELBO}(q, \theta) = \mathbb{E}_q\bigl[\log p(X, s \mid \theta)\bigr] + H(q),
\qquad
H(q) = -\,\mathbb{E}_q\bigl[\log q(s)\bigr]

KLダイバージェンスは常に0以上なので，ELBOは対数尤度の下界になります．そしてKLが0になるのは，$q$が事後分布$P(s \mid X, \theta)$に一致するときだけです．

直接触れない対数尤度の代わりに，この下界を$q$と$\theta$で交互に持ち上げます．

• Eステップ — $\theta^{\text{old}}$を固定して$q$について最大化します．最適解は$q(s) = P(s \mid X, \theta^{\text{old}})$で，下界が対数尤度に接します
• Mステップ — その$q$を固定して$\theta$について最大化します

すると1往復で

\log p(X \mid \theta^{\text{old}})
= \mathrm{ELBO}(q^*, \theta^{\text{old}})
\le \mathrm{ELBO}(q^*, \theta^{\text{new}})
\le \log p(X \mid \theta^{\text{new}})

となり，対数尤度は決して下がりません（左の等号はEステップで接したこと，真ん中はMステップで上げたこと，右はELBOが常に下界であることです）．これが「EMでよい」ことの保証です．ただし保証されるのは対数尤度の単調非減少までで，大域最適解への収束は保証されません．初期値によって異なる局所的な停留点に収束する可能性があります．

では，Eステップでは具体的に何を計算するのでしょうか．

Eステップ ― 事後確率γとξを求める

Eステップの答えは$q(s) = P(s \mid X, \theta^{\text{old}})$，つまり状態列の事後分布です．ただし$2^T$通りの状態列すべてに確率を割り当てる必要はありません．後のMステップで実際に使うのは，周辺化した2つの事後確率だけです．

\gamma_t(k) = P(S_t = k \mid X_{1:T},\, \theta^{\text{old}}),
\qquad
\xi_t(i,j) = P(S_t = i,\, S_{t+1} = j \mid X_{1:T},\, \theta^{\text{old}})

$\gamma_t(k)$は「時刻$t$に状態$k$だった確率」，$\xi_t(i,j)$は「$t$で$i$，次の時刻に$j$だった確率」です．0/1ではなく0〜1の実数で，これが「見えない状態を確率で埋める」の中身です．

では，この$\gamma$と$\xi$をどう計算するのでしょうか．定義どおりでは，またもや$2^T$本の和が要ります．ここでforward-backwardアルゴリズムが登場します．

forward-backward ― γとξをO(T)で求める

条件付き確率の定義から，分子の同時確率に集中します（分母は分子を$k$で足せば出ます）．その分子を，時刻$t$を境に乗法定理で切ります．マルコフ性（$S_t$さえ分かれば未来の観測は過去の観測に依存しない）で条件を整理すると，きれいに2つの因子へ分かれます．

p(X_{1:T},\, S_t = k)
= \underbrace{p(X_{1:t},\, S_t = k)}_{\alpha_t(k)}\;
  \underbrace{p(X_{t+1:T} \mid S_t = k)}_{\beta_t(k)}

• $\alpha_t(k)$（forward変数）：これまでの観測がすべて起き，いま$k$にいる確率．過去からの証拠を運びます
• $\beta_t(k)$（backward変数）：いま$k$にいると仮定したとき，この先の観測が起きる確率．未来からの証拠を運びます

$\alpha$と$\beta$の定義は天下りではなく，求めたい量を乗法定理で切ったら自然に現れたものです．どちらも漸化式で1パスずつ計算できます．

forwardは先頭から右へ進みます．

\alpha_t(k) = \Bigl(\underbrace{\textstyle\sum_i \alpha_{t-1}(i)\, A_{ik}}_{\text{予測：遷移の全パターンを足す}}\Bigr)\;
\underbrace{b_t(k)}_{\text{更新：観測を取り込む}}

「1つ前の状態で場合分けして（全確率の公式），遷移し，そこで自分の観測を取り込む」という読み方です．初期値は$\alpha_1(k) = \pi_k\, b_1(k)$です．

backwardは末尾から左へ進みます．

\beta_T(k) = 1,
\qquad
\beta_t(k) = \sum_j
\underbrace{A_{kj}}_{k \to j \text{ へ遷移し}}\;
\underbrace{b_{t+1}(j)}_{\text{そこで観測し}}\;
\underbrace{\beta_{t+1}(j)}_{\text{その先が続く}}

初期値が全状態で1なのは，時刻$T$より先には観測が存在せず，未来側から掛ける尤度因子がないからです．したがって，積の単位元である1を初期値として$\beta_T(k) = 1$とします．

forwardとbackwardは，どの観測を誰が持つかをちょうど分担しています．時刻$t$で掛け合わせると，全観測を1回ずつ使ったことになります．

観測列       x₁   x₂   ┄┄   x_t │ x_{t+1}   ┄┄   x_T
                                │
forward     α₁ ─▶ α₂ ─▶ ┄ ─▶ α_t                     過去を蓄積しながら右へ
backward                       β_t ◀─ ┄ ◀─ β_T (=1)  未来を蓄積しながら左へ
合流                    γ_t ∝ α_t · β_t              時刻 t で全観測が揃う

掛けて正規化すれば$\gamma$が，そこに遷移と観測を1個挟めば$\xi$が出ます．

\gamma_t(k) = \frac{\alpha_t(k)\,\beta_t(k)}{\sum_{k'} \alpha_t(k')\,\beta_t(k')},
\qquad
\xi_t(i,j) \propto \alpha_t(i)\, A_{ij}\, b_{t+1}(j)\, \beta_{t+1}(j)

これで$2^T$通りの状態列を列挙する必要がなくなりました．計算量は一般の$K$状態HMMでは$O(TK^2)$，今回のような2状態では$O(4T) = O(T)$です．課題1はここで解消します．

$\gamma$と$\xi$が手に入りました．では，これを使ってパラメータをどう更新するのでしょうか．

Mステップ ― Q関数を解析的に解く

ELBOの第1項をQ関数と呼びます．

Q(\theta \mid \theta^{\text{old}}) = \mathbb{E}_q\bigl[\log p(X, s \mid \theta)\bigr]
\qquad \text{よって} \qquad
\mathrm{ELBO}(q, \theta) = Q(\theta \mid \theta^{\text{old}}) + H(q)

Mステップでは$q$を固定するので$H(q)$は定数です．したがってELBOの最大化は$Q$の最大化に帰着します．$Q$を展開すると，各項の重みがちょうど$\gamma$と$\xi$になります．

Q(\theta \mid \theta^{\text{old}}) =
\sum_k \gamma_1(k)\log \pi_k
+ \sum_{i,j} \Xi_{ij} \log A_{ij}
+ \sum_{t}\sum_k \gamma_t(k) \log b_t(k),
\qquad \Xi_{ij} = \sum_{t}\xi_t(i,j)

ここが要点です．課題2で失った「パラメータごとに項が分かれる」という性質が，期待値の世界では復活しています．$\pi$の項・$A$の項・$(\mu, \sigma)$の項が独立なので，別々に解けます．しかも解析解があります．勾配法のような数値最適化は要りません．

• $\pi$と$A$には「和が1」という制約があるので，ラグランジュ乗数法で解きます
• $\mu$と$\sigma$は制約がないので，素直に微分して0と置くだけです

結果はこうなります．

\pi_k \leftarrow \frac{\gamma_1(k)}{\sum_{k'} \gamma_1(k')},
\qquad
A_{ij} \leftarrow \frac{\Xi_{ij}}{\sum_{j'} \Xi_{ij'}}

\mu_k^{\text{new}} = \frac{\sum_t \gamma_t(k)\, X_t}{\sum_t \gamma_t(k)},
\qquad
(\sigma_k^2)^{\text{new}} = \frac{\sum_t \gamma_t(k)\,\bigl(X_t - \mu_k^{\text{new}}\bigr)^2}{\sum_t \gamma_t(k)}

分散には更新後の$\mu_k^{\text{new}}$を使います．

この形の意味は，状態が見えていた場合の答えと並べると分かります．

パラメータ	状態が見えていたら	Mステップの更新式
$\pi_k$	最初の状態が$k$だったか（0か1）	最初の時刻の状態確率 $\gamma_1(k)$
$A_{ij}$	$i \to j$の遷移回数 ÷ $i$にいた回数	遷移の期待回数 $\Xi_{ij}$ ÷ $i$にいた期待回数
$\mu_k$	状態$k$だった時刻の観測値の平均	$\gamma$を重みにした加重平均
$\sigma_k^2$	同じ時刻の分散	$\gamma$を重みにした加重分散

「回数」が「期待回数」に，「平均」が「加重平均」に変わっただけです．$\gamma$がすべて0/1なら，左右の列は完全に一致します．

$\mu_k$の式は課題2の不動点方程式と同じ形ですが，意味が違います．あちらは重みが未知数に依存したままの解けない式でした．こちらは$\gamma$を$\theta^{\text{old}}$で固定してあるので，右辺はただちに計算できます．循環を反復に変えたわけです．

系列が複数に分かれる場合

実データでは，記録が長時間途切れることがあります．途切れた前後を「1ステップの遷移」として数えると遷移確率$A$が歪むため，ギャップで系列を切り，独立した複数の観測列として扱います（後半の実験では332区間に分かれました）．

このとき更新式は，区間$m = 1, \dots, M$にわたって十分統計量を足し合わせる形になります．

\pi_k^{\text{new}} = \frac{1}{M}\sum_{m=1}^{M} \gamma_1^{(m)}(k),
\qquad
A_{ij}^{\text{new}} = \frac{\sum_m \Xi^{(m)}_{ij}}{\sum_{j'}\sum_m \Xi^{(m)}_{ij'}}

$\mu$と$\sigma$も同様に，全区間の$\gamma$で重み付けた和を取ります．要点は$\pi$が「各区間の先頭がどちらだったかの平均」になることと，区間をまたぐ$\xi$は計算しないことです．そのため，区間数が多いほど$\pi$は0か1かに寄らず中間的な値になります．

学習の全体の流れ

部品が揃ったので，学習の全体はこうなります．

1. 初期値θを置く
2. E ステップ: いまのθでforward-backwardを回し，γ, ξを求める
3. M ステップ: γ, ξを重みにθを更新する（解析解）
4. 対数尤度の改善が閾値を下回るまで2〜3を繰り返す

実用上の注意が2つあります．

• 初期値． 前述のとおりEMは初期値によって異なる停留点に収束しうるため，初期値が結果に効きます．観測値の下位25%の平均を$\mu_{\mathrm{OFF}}$，上位25%の平均を$\mu_{\mathrm{ON}}$に置き，2つの山のそれぞれの側から出発させます
• **ラベルの固定．**EMが返す状態番号に意味はありません（ラベルスイッチング）．学習後に$\mu$の昇順で並べ替え，小さい方をOFFと呼びます．並べ替えるのは名前だけで，境界の位置は変わりません

推論 ― 学習済みHMMでON/OFFを求める

学習が終われば$\theta$は固定です．ON/OFFの判定は，学習済みの$\theta$でforward-backwardをもう一度回すだけです．得られる$\gamma_t(k) = P(S_t = k \mid X_{1:T})$が求めたかった状態確率で，0.5を超えた時刻をONと判定します．

状態確率には2種類あり，区別が必要です．

	定義	使う観測	用途
filtering	$P(S_t \mid X_{1:t})$	現在時刻まで	リアルタイム判定
smoothing	$P(S_t \mid X_{1:T})$	未来も含む	過去ログの解析

過去ログの解析では未来の観測も使えるので，smoothingを使います．OFFとONの分布がよく分離している列では両者はほぼ一致し，差が出るのは観測が曖昧なときです．

なお，定義どおりに実装すると$\alpha$は確率の積なので急速に0へ近づき，150万点では倍精度の下限を割って計算が壊れます．各時刻で正規化しながら進めるスケーリングという定石で回避できます（結果は変わりません）．

MetroPT-3で試す

データセット

MetroPT-3（UCI Machine Learning Repository，CC BY 4.0）を使います．地下鉄車両の空気圧縮機（APU）のセンサログです．

項目	値
行数	1,516,948
期間	2020-02-01〜2020-09-01
サンプリング	約10秒間隔
列	アナログ7列（圧力・温度・電流）＋デジタル8列

注意点は時間軸の欠測です．50秒を超えるギャップが331箇所あり，期間の約17.6%が欠けています．途切れた前後を1ステップの遷移として数えると$A$が壊れるので，ギャップで系列を分割し，独立な観測列の集まりとして学習します．全期間では332区間に分かれました．

対象列はコンプレッサ本体の圧力TP2です．全期間を時間軸の中点で2等分し，期間ごとに独立に学習します．

	期間	行数	区間数
前半	2020-02-01 00:00〜2020-05-16 23:43	756,844	159
後半	2020-05-17 01:31〜2020-09-01 03:59	760,104	173

推定されたパラメータ

初期状態確率と観測分布です．どちらの期間も10回前後の反復で収束しました（前半10回，後半6回）．

パラメータ	前半	後半	読み方
$\pi_{\mathrm{OFF}},\ \pi_{\mathrm{ON}}$	0.491, 0.509	0.416, 0.584	区間の先頭がどちらだったか．区間が150以上あるので中間的な値になる
$\mu_{\mathrm{OFF}}$	−0.0134 bar	−0.0121 bar	停止中はほぼ0（ゼロ点オフセットぶん負）
$\mu_{\mathrm{ON}}$	8.6522 bar	8.6520 bar	負荷運転中の圧力
$\sigma_{\mathrm{OFF}}$	0.0027 bar	0.0029 bar	停止中はほぼ一定
$\sigma_{\mathrm{ON}}$	1.9470 bar	1.6803 bar	起動・停止の立ち上がりを含むため広い

遷移確率は行列で見ると分かりやすくなります．行が遷移元，列が遷移先で，各行の和は1です．

A^{\text{前半}} =
\begin{array}{c|cc}
 & \text{OFF} & \text{ON} \\ \hline
\text{OFF} & 0.9934 & 0.0066 \\
\text{ON}  & 0.0457 & 0.9543
\end{array}
\qquad
A^{\text{後半}} =
\begin{array}{c|cc}
 & \text{OFF} & \text{ON} \\ \hline
\text{OFF} & 0.9899 & 0.0101 \\
\text{ON}  & 0.0437 & 0.9563
\end{array}

対角成分がどちらも0.95以上で，状態が続きやすいことが読み取れます．ここから導かれる量が次の表です．

派生量	前半	後半
OFFの期待継続 $1/(1 - A_{\mathrm{OFF},\mathrm{OFF}})$	150.6サンプル（約25分）	98.6サンプル（約16分）
ONの期待継続 $1/(1 - A_{\mathrm{ON},\mathrm{ON}})$	21.9サンプル（約3.7分）	22.9サンプル（約3.8分）
稼働率（推定ONの割合）	12.83%	19.01%

閾値をどこにも書いていないのに，「停止はほぼ0 bar・稼働は8.65 bar・1回動くと約4分・止まると16〜25分」という設備の姿がパラメータとして出てきます．

興味深いのは2期間の差です．$\mu_{\mathrm{ON}}$は8.6522 → 8.6520とほぼ変化しない一方，稼働率は12.83% → 19.01%と大きく増えています．OFFの期待継続が25分から16分に縮んだこととも整合します．冒頭で述べたとおり，状態を分けずに平均を取ると，この「稼働率の変化」が「値の変化」に見えてしまいます．実際，全体平均は1.100 → 1.635 bar（+48.7%）と大きく動きました．

可視化

推定結果を元の波形に重ねたものが次の図です．赤がON，灰色がOFFと推定された区間で，上段が前半，下段が後半を示しています．

赤は圧力が持ち上がっている区間を覆っていて，立ち上がり・立ち下がりもONに含まれています．ONの観測値が8.65 bar付近に限らず0〜10 barに散らばるため，$\sigma_{\mathrm{ON}}$は$\sigma_{\mathrm{OFF}}$より3桁近く大きくなります．

また，下段のほうがONが密です．稼働率12.83% → 19.01%という差が，波形の上でも見て取れます．

おわりに

今回やったことをまとめます．

• センサ値だけから稼働/停止を切り分けるために，2状態Gaussian HMMを使いました．
• パラメータは最尤推定で決めます．状態が見えないため直接は解けず，ELBOを根拠とするEMアルゴリズムで反復的に推定しました．Eステップで事後確率$\gamma, \xi$を求め，Mステップで解析解により更新します．
• $\gamma, \xi$の計算にはforward-backwardを使いました．過去の証拠（$\alpha$）と未来の証拠（$\beta$）を別々に蓄積して掛け合わせることで，$2^T$通りの列挙が$O(TK^2)$（2状態なら$O(T)$）になります．
• 推論は，学習済みパラメータでforward-backwardをもう一度回すだけです．smoothing $P(S_t \mid X_{1:T})$が0.5を超えた時刻をONとしました．
• MetroPT-3の圧力データを2期間に分けて学習すると，稼働中の圧力$\mu_{\mathrm{ON}}$は8.6522 → 8.6520 barとほぼ不変な一方，稼働率は12.83% → 19.01%に増えていました．状態を分けない全体平均はこれを +48.7% の増加として評価してしまいます．

最初にHMMを使う利点を挙げましたが，正直なところ「使ってみたかった」というのが一番の動機でした笑

数式の理解が曖昧なところも残っているので，誤りやお気づきの点があればコメントで教えていただけると助かります．

最後まで読んでいただきありがとうございました．

参考文献

• Rabiner, L. R. (1989). A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proceedings of the IEEE, 77(2).
• Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Chapter 13.
• Davari, N., Veloso, B., Ribeiro, R. P., Gama, J. — MetroPT-3 Dataset, UCI Machine Learning Repository (CC BY 4.0)