数学モデルのパラメータを推定するアニメーションを作る
この記事はこんな人にオススメです
  • モデルのパラメータ推定をしている人
  • 理系大学生の人（アニメーションは資料作成等にお役立てください）
  • 数学や工学が好きな人

こんにちは．けんゆーです（@kenyu0501_)
Twitterの方でも発信しましたが，以下．

時系列データの解析で良く行われる「時間周波数解析」と，「その解析窓の内側で最小二乗法を使ったモデルパラメータ推定」のプログラミングを公開しようと思います．
時系列データの解析をリアルタイムで行えるようにするためアニメーションを作っています．
（アニメーションというか最終的に解析したデータの表示を考えたいところですが，，，）

時間周波数解析は，生体信号を解析する上でも基本的なものなので知っていて損はないと思います．
周波数が高くなれば緊張とか，低くなれば安静とか，そういった定性的な指標を獲得するためにも大事なツールの一つです．
（参考：「時間周波数分析とは！？」）

さらに，最小二乗法は計算コストが低い上に，データの傾向を把握するモデルのパラメタを決定するアルゴリズムでして，こいつも同時に使えたらデータ解析が面白くなると思います．（ほとんどの工学系では，一次直線近似とかで習うと思います）

解析のアニメーション動画について

改良を重ねて，アニメーションももっと見やすい形に随時変更中ですが，出し惜しみはするものではないと思うのでその都度UPしていきます．

以下では，このアニメーションの説明をしておきます．
（ブログに載せているため，GIFにしてますが，プログラムを回すとmp4出力です．）

一番上の図は生体信号（脳波）の生波形です．
フィルタリング処理など，何も施していない時間的なデータです．

この時系列データに対して，解析するための窓（赤い線）を設定します．
窓幅は１秒間（512 sampling点数）です．
もちろん，自分のデータの解析したい時間に調整することが可能です．
この解析する窓は，ランニングウインドウとか，スライディングウインドウと呼ばれます．
窓内での解析処理を終えると，1秒分横に移動して，さらに次の解析を行います．

このウインドウごとに高速フーリエ変換を行います．
このプログラムでは，ランニングウインドウごとに支配的な周波数を二つ選びます．
（対象が脳波なので，0.5 ~ 30.0 Hz の間にしてます．）

ちなみにモデルはこれです．

アニメーションが上手く動くかどうか試しただけなので，自分で好きな数式を入れたら良いと思います．

選んだ二つのピーク（ω1とω2）を，3×3で表示した9つのグラフのうち，青のプロットで示してます．
モデル右辺のωに入れて，残りは推定します．
あっ，ωは周波数ピーク値に2πをかけてます．

残りのパラメータ(A, B, C, P1_1, P2_1, P2_1, P2_2)は，最小二乗法によってパラメータ同定を行なってます．
赤で示したプロットです．
ここで，生波形は使わずに，ピークを含んだ以下の周波数帯域の情報を用いたIFFTの波形によってパラメタを決めてます．

10個に帯域を分けてます．
今回のパラメタ推定の詳しい狙いなどは，以下の奮闘記に書いています．
「最小二乗法によるパラメータ同定奮闘記2」

Pythonプログラムについて

こちらがプログラムになります．
公開してますが，他のユーザが使用するのを意識してないので雑です．
素人が組んだプログラムなので素人には読めると思います（笑）
（できないよりは回れば良いと思っている派なので勘弁ください，，，）
クラスとかそういうオブジェクト指向たっぷりの組み方もいつか勉強しておきます．

やりたいことが無くなって時間が余ったら見やすいコードに修正します．

アニメーションの意義

理系大学生ならば，発表資料等で，「解析方法」や「実験結果」を説明すると思います．
やはりその時には，動きのある資料の方が良かったりする場合が多いと思います．
価値の高い素敵な発表も，難しすぎて聴衆に理解されないままだと悔しいので，出来る限り納得させる工夫ができたら良いですね．

ブログに公開するのを了承してくださる方で，そういったご相談があれば，一緒に考えますよ〜！
Twitterの方にDM下さい．

では〜！