<第 20 回 ITS シンポジウム 2022>

地域全体のアクティビティを再現するナウキャストモデルの構築
花房比佐友*1

小宮粋史 2*1

一瀬恭平 2*1

堀口良太 2*1

株式会社アイ・トランスポート・ラボ*1

論文概要

本研究は、都市・地域全体を対象として、直近時間帯における「滞在人口」等実測データに基づいて個人のアクティビティを再現するためのナウキャストモデルを構築することを目的とする。本研究においては、対象エリア内におけるトリップチェーンと滞在行動で構成されるアクティビティの総需要を未知変数として、ある時間帯・地域における滞在人口、道路・公共交通ネットワーク上の断面交通量などの観測値を制約条件として関係をエントロピーモデルで定式化し、アクティビティ需要の推計を行う手法を考案した。本稿では、ナウキャストモデルのコンセプトおよび手法、リアルタイム推計にあたっての課題と対応について議論する。

Development of an Activity-Based Nowcast Model for RegionWide
Hisatomo Hanabusa*1

Tadashi Komiya*1

Kyohei Ichinose*1

Ryota Horiguchi*1

i-Transport Lab. Co., Ltd.*1

Abstract
This study aims to develop a method of an activity-based nowcast model for region wide based on real-time data suchas estimated population distributions and traffic volumes. This method uses an entropy model to estimate individual’s
daily activity demands consisting of trip chains and stays as unknown variables, with the estimated population
distribution and traffic counts in real-time as constraints. In this paper, we describe the concept and the methodology ofthe model. On the other hand, we discuss the issue to be tackled for real-time estimation.

Keyword: Activity-based simulation, Trip chain, Nowcast, Entropy maximization method

1. はじめにMaaS（Mobility as a Service）の推進は、近年多様化する地域住民や旅行者などの移動ニーズに対応するサービスの形態として、社会実装に向けた様々な取り組みがなされている。MaaS は電車・バス等の公共交通をはじめとした多様な移動サービスを組み合わせて検索、予約や決済などを一括で行うサービスであるといえるが、どの地域でも同じ手続きで利用できるシステムであることが望ましい一方、交通インフラや抱えている社会課題は地域によって大きく異なることを踏まえた対応が必要である。当然のことながら、個人の行動パターンもまた地域性があり、通勤や買物など日々の行動を行うための移動手段が異なることから、提供されるサービスもその地域に最適化されたものが提案されなければならない。
2019 年 12 月初旬に第 1 例目の感染者が報告されたことを端に世界的流行（パンデミック）となった
新型コロナウイルス感染症（COVID-19）1は、国内のみならず世界中で日々の生活に大きな影響を与えてきた。日本国内においては、パンデミック初期における外出自粛や、その後テレワーク導入など生活様式の変化によって、外出・移動量はパンデミック以前にくらべて著しく低下した2。今までの消費行動を前提とした社会に大きな影響を与えることになり、経済刺激策と感染拡大防止策の両立を試みる中で、モノやサービスに対する需要と供給のバランスが不安定な状態となっている。このことをきっかけに、今後も人々のライフスタイルへの嗜好性変化はあり続けると考えられ、職住や余暇活動の多様化に適用した企業活動や都市構造のあり方を継続して議論していく必要がある。加えて、様々なリソースの規模や粒度が小さくなっていくことが予想され、大きな混乱を引き起こすことのないよう直近の需給ギャップの変化に対応できる仕組みも求められる。このよ
うな社会的変化に対応し、交通環境を最適化することをと可能とする MaaS の実現には、より時間的・空間的に高い解像度で個人の行動を把握・理解し、今後起こりうる結果を予測しながら対策を打つ手段や方針を検討するための環境整備が必要であると考えられる。
本研究における目的は、観測された人流の状況から地域全体の各個人の行動（アクティビティ）需要がどのように変化したかを推定し、可視化する手法を確立することである。加えてその推定を逐次リアルタイムで行うことを試みる。これを本稿ではマクロナウキャスト手法と呼ぶ。ここで提案する方式は、前述した交通環境を最適化するための個人の行動変化の理解を支援するものであるが、同時に、近未来を予測するための機能を持たせて、需要の変化に対応する策を評価するためのシミュレーションを実現することも目指す。第 2 章では、ナウキャスト手法の概要として本研究で定義する個人のアクティビティや本手法のコンセプトを説明する。第 3 章では、実際にアクティビティを推計するモデルを示し、解法の導出を行う。第 4 章では、リアルタイム処理を行う上での課題とその対処について議論し、最後に5 章で今後の方針について述べる。
2. マクロナウキャスト手法
マクロナウキャスト手法は、対象地域で直近時間帯の実測値が示す状況を再現するよう個人のアクティビティ需要を求め、その状況をシミュレーション上で再現・可視化する仕組みである。以下、アクテ
ィビティの定義とマクロナウキャストの方針について述べる。
2-1.アクティビティの定義
本研究において定義する個人のアクティビティとは、ある個人 1 日の行動を拠点間の移動（トリップ）と滞在の組み合わせで構成するものである。トリップは、出発点、出発時刻、移動手段、目的地、到着時刻で構成されたもので定義される。一方、滞在については、同じ位置における滞在開始時刻（＝トリップの到着時刻）と滞在終了時刻（＝次のトリップの出発時刻）で表現され、滞在時間は滞在終了時刻と滞在開始時刻の差分で計算できる。このことから、
個人𝑎の行動を定義するアクティビティ𝐴𝑎 は以下の通り記述できる。
𝐴𝑎 = {(𝑡1𝑎 , 𝑝1𝑎 , 𝑐1𝑎 ), … , (𝑡𝑖𝑎 , 𝑝𝑖𝑎 , 𝑐𝑖𝑎 ), … , (𝑡𝑛𝑎 , 𝑝𝑛𝑎 , 𝑐𝑛𝑎 )} (1)
ここで、𝑡𝑖𝑎 ：アクティビティ上の𝑖番目の時刻
𝑝𝑖𝑎 ：アクティビティ上の𝑖番目の位置
𝑐𝑖𝑎 ：アクティビティ上の𝑖番目の交通モード
アクティビティを構成する各点の時刻においては、秒単位での記述やある単位時間で区切った時間帯で構成してもよい。また位置においては緯度経度などの座標あるいはメッシュ等で区切られた空間上の定義で与えることができる。交通モードにおいては、徒歩、バス、鉄道、自動車など、個人が選択し得る手段を任意に設定する。1 日のアクティビティを構成する各点の時間解像度においては、数秒単位から数時間単位、空間解像度においては、数 m 単位から数 km 単位など様々であるが、自宅から勤務地、勤務地から商業施設など、滞在場所が把握できるレベ

ルで定義されている必要がある。また駅の利用者数など、通過交通量による調整が必要な場合は、経路が把握できる解像度が必要となる。図 1 にアクティビティデータ時空間図上に展開したイメージを示す。

図 1 アクティビティデータのイメージ
2-2.マクロナウキャストの方針
マクロナウキャスト手法は、現在の人流状況に応じて、対象地域のアクティビティ需要をリアルタイム推計し、可視化することを目的とする。全体のシステムとしては、まず直近の観測値が入力されることをトリガーにして、観測値の状況を再現できるようアクティビティ需要を推定する。その後、推定したアクティビティ需要をシミュレーションモデルに入力し、個人の行動を模擬して地域全体の状況を可視化する。図 2 に処理手順を示す。に推定されたアクティビティ需要を基に再推定され、逐次過去のアクティビティ需要のパターンが更新されていく。
ここで、前提としているのは、需要推計に用いられるアクティビティパターンが既知の情報としてあらかじめ用意されることで、ここではパーソンプローブデータ等の実データに基づいて作成される疑似データを適用することを想定している。システムの運用を考えると、平日、休日における平均的な状況やイベント時の特異日など、過去のデータから想定し得るいくつかのアクティビティパターンを用意しておくことが望ましい。本手法においては、マクロナウキャストを行うにあたって各個人のアクティビティにおける出発時間帯や滞在場所についての変更を考慮していないが、この方式を採用することによって、過去から現在までの交通状況を再現することができると同時に、数時間後の交通状況も予測（フォアキャスト）することができるようになる。
この仕組みを実現するためには、アクティビティ需要を推計するモデルと、推計したアクティビティ需要を基に各個人の行動を模擬し、地域全体の混雑状況等を再現するシミュレーションモデルが必要になる。ここで想定するシミュレーションモデルは、都市部全域など広域で個人の 1 日の行動をシミュレーションするマルチエージェントモデルで、該当するものとしてはたとえば MATSim3が挙げられる。
3. アクティビティ需要の推計モデル
本章では、前章で述べたアクティビティ需要の推定モデルについて説明する。図 3 にアクティビティ需要推計のイメージを示す。

図 2 マクロナウキャスト処理
ここである地域において 1 時間毎に滞在人口が観測されている場合を想定する。マクロナウキャスト
では、その地域に関わる（滞在、通過など）アクティビティパターンを基に、観測状況に合う需要を 1時間毎に推定する。このとき、今期（直近 1 時間）に推定するアクティビティ需要は、1 期前（2 時間前）

図 3 アクティビティ需要の推計イメージ本手法では、対象エリア内のアクティビティ総需要を未知変数として、ある時間帯・地域における滞在人口、交通ネットワーク上の断面交通量の観測値に一致させることを制約条件とし、エントロピー最大化法によりアクティビティ需要を逆推定する。
エントロピー最大化法による交通需要の推定手法については、たとえば Willumsen4、小根山ら5、三輪ら6が提案したモデルが挙げられる。本研究においては、小根山らのモデルを参考に、1 日の行動が記述されたアクティビティ需要を推定する目的で構築した確率密度関数𝑃を以下に示す。

観測値で定義している時間帯ℎは、マクロナウキャストを行う現在から過去の時間帯であり、推定のために利用する直近の過去に取得された観測値を利用する目的で導入する。次に𝑃について対数関数を求め、スターリングの公式を適用することによって以下の式を得ることができる。
log 𝑃 ≈ − ∑{𝑞𝜔 log 𝑞𝜔 − 1}𝑚+ ∑ 𝑞𝜔 𝑙𝑜𝑔 𝑞̂𝜔𝜔+ 𝛾1 [− ∑{𝑣𝑎 (ℎ) 𝑙𝑜𝑔 𝑣𝑎 (ℎ) − 1}
𝑎,ℎ𝑞𝜔

𝑃=[𝑞!
𝑞̂𝜔∏()∏𝜔 𝑞𝜔 !
∑𝜔 𝑞̂𝜔

]𝜔+ ∑ 𝑣𝑎 (ℎ) 𝑙𝑜𝑔 𝑣̂𝑎 (ℎ)]

𝑣!
𝑣̂𝑎 (ℎ)∙[∏()∏𝑎,ℎ 𝑣𝑎 (ℎ)!
∑𝑎,ℎ 𝑣̂𝑎 (ℎ)

𝑣𝑎 (ℎ) 𝑟1𝑎,ℎ]𝑎,ℎ∙ [∏𝜎!

𝑧,ℎ 𝜎𝑧 (ℎ)!

∏𝑧,ℎ (𝜎̂𝑧 (ℎ)∑𝑧,ℎ 𝜎̂𝑧 (ℎ))𝜎𝑎 (ℎ) 𝑟2]+ 𝛾2 [− ∑{𝜎𝑧 (ℎ) 𝑙𝑜𝑔 𝜎𝑧 (ℎ) − 1}
𝑧,ℎ…(2)+ ∑ 𝜎𝑧 (ℎ) 𝑙𝑜𝑔 𝜎̂𝑧 (ℎ)] + 𝑐𝑜𝑛𝑠.
𝜔𝑧,ℎ𝑣 = ∑ 𝑣𝑎 (ℎ)

𝑞 = ∑ 𝑞𝜔…(3)𝑎,ℎ…(4)𝜎 = ∑ 𝜎𝑧 (ℎ)
𝑧,ℎ…(5)

ここで、𝑞：アクティビティ推定需要の総量𝑞𝜔 ：アクティビティ𝜔の推定需要𝑞̂𝜔 ：アクティビティ𝜔の先験需要
𝑣：推定観測交通量の総量
𝑣𝑎 (ℎ)：計測箇所𝑎、時間帯ℎの推定観測交通量𝑣̂𝑎 (ℎ)：計測箇所𝑎、時間帯ℎの観測交通量
𝜎：推定観測滞在人口の総量𝜎𝑧 (ℎ)：ゾーン𝑧内、時間帯ℎの推定観測滞在人口𝜎̂𝑧 (ℎ)：ゾーン𝑧内、時間帯ℎの観測滞在人口
𝑟1、𝑟2：重み付けパラメータ𝑞𝜔 や𝑞̂𝜔 で定義されるアクティビティ𝜔は、あらかじめ各地点の時刻、位置、交通モードが用意されたもので、OD 交通量の OD と同様一つのパターンとして扱い、出発時刻等は推定対象としない。また𝑞̂𝜔
は先験需要であるが、マクロナウキャストを行う場
合は、1 期前に推定したアクティビティ需要に置き換えて当期のアクティビティ需要を推定する。一方、…(6)

ここで、上記の式を次の条件下で最大化することを考える。
𝑣𝑎 (ℎ) = ∑ 𝑝𝑎𝜔 (ℎ)𝑞𝜔…(7)𝜔𝜎𝑧 (ℎ) = ∑ 𝑝𝑧𝜔 (ℎ)𝑞𝜔…(8)𝜔𝑝𝑎𝜔 (ℎ)：リンク a を時間帯ℎに通過するアクティビティ𝜔の選択確率𝑝𝑧𝜔 (ℎ)：ゾーン z で時間帯ℎに滞留する
アクティビティ𝜔の選択確率
この問題は、ラグランジュの未定乗数法によって解を求めることができる。ラグランジュ関数Lは以下のように定義することができる。

L = log 𝑃 + ∑ 𝜆𝑎 (ℎ) {𝑣𝑎 (ℎ) − ∑ 𝑝𝑎𝜔 (ℎ)𝑞𝜔 }𝑎,ℎ𝜔+ ∑ 𝜇𝑧 (ℎ) {𝜎𝑧 (ℎ) − ∑ 𝑝𝑧𝜔 (ℎ)𝑞𝜔 }
𝑧,ℎ

𝜔…(9)

なお、Lで導入した𝜆𝑎 (ℎ)および𝜇𝑧 (ℎ)はラグランジュ乗数である。ここで、一階導関数が満たすべき最適条件（Karush-Kuhn-Tacker 条件）は、以下の通りとなる。

−1𝜎̂𝑧 (ℎ)𝑌𝑧 (ℎ) 𝛾2 = ∑ 𝑝𝑧𝜔 (ℎ) ∙ 𝑞̂𝜔 ∏ 𝑋𝑎 (ℎ)𝑝𝑎𝜔(ℎ)
𝜔𝑎,ℎ∙ ∏ 𝑌𝑧 (ℎ)𝑝𝑧𝜔 (ℎ)
𝜕𝐿= − 𝑙𝑜𝑔 𝑞𝜔 + 𝑙𝑜𝑔𝑞̂𝜔 − ∑ 𝜆𝑎 (ℎ)𝑝𝑎𝜔 (ℎ)
𝜕𝑞𝜔𝑎,ℎ− ∑ 𝜇𝑧 (ℎ)𝑝𝑧𝜔 (ℎ) = 0…(10)𝑧,ℎ𝜕𝐿= −𝛾1 𝑙𝑜𝑔 𝑣𝑎 (ℎ) + 𝛾1 log 𝑣̂𝑎 (ℎ) + 𝜆𝑎 (ℎ) = 0
𝜕𝑣𝑎 (ℎ)…(11)𝜕𝐿= −𝛾2 𝑙𝑜𝑔 𝜎𝑧 (ℎ) +𝛾2 𝑙𝑜𝑔 𝜎̂𝑧 (ℎ) + 𝜇𝑧 (ℎ) = 0
𝜕𝜎𝑧 (ℎ)…(12)
さらに−𝜆𝑎 (ℎ)と−𝜇𝑧 (ℎ)から以下の式をおくこと
で、アクティビティ𝜔の推定需要𝑞̅𝜔 、観測箇所𝑎の推定観測交通量𝑣̅𝑎 (ℎ)、ゾーン𝑧の推定観測滞在人口𝜎̅𝑧 (ℎ)は𝑋𝑎 (ℎ)と𝑌𝑧 (ℎ)の関数にすることができる。
𝑋𝑎 (ℎ) = 𝑒𝑥𝑝[−𝜆𝑎 (ℎ)] …(13)𝑌𝑧 (ℎ) = 𝑒𝑥𝑝[−𝜇𝑧 (ℎ)] …(14)よって、𝑞̅𝜔 = 𝑞̂𝜔 ∙ ∏𝑎,ℎ 𝑋𝑎 (ℎ)𝑝𝑎𝜔(ℎ) ∙ ∏𝑧,ℎ 𝑌𝑧 (ℎ)𝑝𝑧𝜔(ℎ)−1𝑣̅𝑎 (ℎ) = 𝑣̂𝑎 (ℎ)𝑋𝑎 (ℎ) 𝛾1−1𝜎̅𝑧 (ℎ) = 𝜎̂𝑧 (ℎ)𝑌𝑧 (ℎ) 𝛾2…(19)𝑧,ℎ…(15)…(16)以上の非線形連立方程式 を解き、𝑋𝑎 (ℎ)および𝑌𝑧 (ℎ)を式（15）に代入することでアクティビティ需要𝑞̅𝜔 を求めることができる。
4. リアルタイム推計の課題と対応
リアルタイム推計においては、決められた更新頻度で計算が完了している必要があり、本研究で提案する本手法においても、計算負荷を軽減する方策を検討する必要がある。本手法においては、前述の通り、観測値を基にアクティビティ需要を求める計算を行うことになるが、対処としては、観測値の集計時間単位や空間解像度を調整する、未知変数であるアクティビティ需要のパターン数を少なくする、などの方法が考えれる。
本研究においては、アクティビティパターンをクラスタリング手法によって集約させ、集約したアクティビティから代表的な需要パターンを定義し、未知変数を少なくすることで計算負荷を低減させる。
前述の通りアクティビティは時間と位置、移動手段から構成され、1 日の行動を 3 次元の時空間上に展開することができる。図 4 にアクティビティパターンの集約時の類似性評価イメージを示す。方針としては、アクティビティを構成する各点の平面空間上の位置が近い、かつ移動時の交通モードが同じである状況がどの程度あるかを評価することでアクティビティパターンを集約化する。

…(17)

なお、𝑋𝑎 (ℎ)および𝑌𝑎 (ℎ)は、式（15）と式（16）を式（7）に、式（15）と式（17）を式（8）に代入した次の式を解くことによって得られる。
−1𝑣̂𝑎 (ℎ)𝑋𝑎 (ℎ) 𝛾1 = ∑ 𝑝𝑎𝜔 (ℎ) ∙ 𝑞̂𝜔 ∏ 𝑋𝑎 (ℎ)𝑝𝑎𝜔(ℎ)
𝜔∙ ∏ 𝑌𝑧 (ℎ)𝑝𝑧𝜔 (ℎ)

𝑎,ℎ…(18)

図 4 アクティビティの類似性評価イメージ𝑧,ℎ

本研究では、類似した行動を評価する手法として
Dynamic Time Warping7（以下、DTW）を適用する。

DTW は、音声データなど時系列で計測される 2 つの波形データの類似性を評価する手法として適用されているもので、位相のずれや時間軸方向のずれを許容しながら波形間の距離を見積もることができる。
ここで DTW による類似性評価方法について説明する。アクティビティ𝐴𝑎 は n 個の移動点を、𝐴𝑏 は m 個の移動点をもつと仮定した場合、それぞれの軌跡上のある 2 地点𝑠𝑖𝑎 、𝑠𝑗𝑏 を直線で結ぶユークリッド距離は、ℝ2 上の xy 座標𝑥(𝑠𝑖𝑎 )、𝑦(𝑠𝑖𝑎 )、𝑥(𝑠𝑗𝑏 )、𝑦(𝑠𝑗𝑏 )に置き換えて以下のように計算できる。
2𝑑(𝑠𝑖𝑎 , 𝑠𝑗𝑏 ) = √(𝑥(𝑠𝑖𝑎 ) − 𝑥(𝑠𝑗𝑏 )) + (𝑦(𝑠𝑖𝑎 ) − 𝑦(𝑠𝑗𝑏 ))2…(20)
類似性評価にあたって留意すべき点は、各個人の交通手段が異なることで、アクティビティの時空間軌跡自体は類似性が高いと評価されていても異なる交通手段で滞在地点へ移動する場合は集約の対象とはならない。また、滞在する地域が同じまたは近くても滞在時間帯が大きく異なる場合は類似性が低いと判断されるべきであろう。そのため、2 地点間の距離を求める際には、同じ交通モードで移動してない、あるいは各点における時刻が大きく乖離している場合は集約の対象としない処理とした。よって制約条件を加えた 2 地点間の距離𝑑 ′ (𝑠𝑖𝑎 , 𝑠𝑗𝑏 )は以下の通りとなる。
𝑑 (𝑠𝑖𝑎 , 𝑠𝑗𝑏 ) = {𝑑(𝑠𝑖𝑎 , 𝑠𝑗𝑏 )′∞𝑎𝑏
𝑖𝑓 𝑐𝑖−1= 𝑐𝑗−1𝑎𝑛𝑑 |𝑡𝑖𝑎 − 𝑡𝑗𝑏 | ≤ 𝜏𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒…(21)

ここで導入した𝜏は、比較する 2 地点間の時間差をどの程度許容するかを判定するパラメータである。
𝐴𝑎 、𝐴𝑏 間の DTW 距離𝐷(𝐴𝑎 , 𝐴𝑏 )は、動的線形計画法を基に以下のように計算ができる。
参考文献
1東京都での新型コロナウイルス感染症（COVID-19）の流行（2020 年 1～5 月）(IASR Vol.41p146-147:2020 年
8 月号),https://www.niid.go.jp/niid/ja/diseases/ka/corona-virus/2019-ncov/2502-idsc/iasr-in/9818486d01.html
2令和３年版交通政策白書,第２部コロナ禍を乗り越え,進化する交通,https://www.mlit.go.jp/sogoseisaku/transport/content/001408994.pdf.
3 MATSim, https://www.matsim.org/
4Willumsen,L G, Estimating time-dependent trip
matrices from traffic counts, Proceedings of the 9th𝐷(𝐴𝑎 , 𝐴𝑏 ) = 𝑓(𝑛, 𝑚)𝑓(𝑖, 𝑗) =…(22)

𝑑 (𝑠𝑖𝑎 , 𝑠𝑗𝑏 ) + min {′𝑓(𝑖 − 1, 𝑗 − 1),𝑓(𝑖 − 1, 𝑗), }𝑓(𝑖, 𝑗 − 1)…(23)

1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚
𝑓(0,0) = 0, 𝑓(𝑖, 0) = 𝑓(0, 𝑗) = ∞𝑓(𝑛, 𝑚)は各地点間のコストを格納する𝑛 × 𝑚の表を示し、最小コストの Warping Path を探索する際に利用される。なお、計算の過程で𝑓(𝑖 − 1, 𝑗 − 1) =𝑓(𝑖 − 1, 𝑗) = 𝑓(𝑖, 𝑗 − 1) = ∞となる場合は、比較したアクティビティはお互いに集約できないと判断される。本研究においては、以上の処理から計算されるDTW 距離を類似性の評価指標として、階層型クラスタリング処理でアクティビティパターンを集約する。
5. 今後の展開
本稿では、広域ナウキャストモデルとして直近の滞在人口や通過交通量を基に各個人のアクティビティ需要を推定する手法について提案し、リアルタイム運用するにあたっての課題と対処方法について議論を行った。今後は疑似人流データを基にした適用性検証を通じて推定精度の評価やリアルタイム処理のための対処方法についての検討とその効果についての検証を行っていく予定である。
本稿は、国立研究開発法人情報通信研究機構の「ウイルス等感染症対策に資する情報通信技術の研究開発 課題 C アフターコロナ社会を形成する ICT 副
題：多様な都市活動を支援する予測情報共有型時空間リソース有効活用技術の研究開発」（222C02）で取り組んだ研究成果である。ここに本プロジェクトの関係者各位に深く謝意を表す。
International Symposium on Transportation andTraffic Theory, pp. 397-411, 1984.
5小根山裕之, 桑原雅夫, 路側観測交通量からの時間変化
する OD 交通量の推定,交通工学,Vol.32,No.2,pp.516,1997.
6三輪富生,山本俊行,竹下知範,森川高行,プローブカーの速度情報を用いた動的 OD 交通量の推定可能性に関する研究,土木学会論文集 D,Vol.64, No.2,pp.252-265,2008.
7H. Sakoe and S. Chiba, Dynamic Programming Algorithm
Optimization for Spoken Word Recognition, IEEE Transactions
on Acoustics, Speech, and Signal Processing, Volume: 26,
Issue: 1, February 1978.