思っちゃったんだからしょうがない

統計学を学ぶ傍でひっそりとやっていきます。 #statistics #R

R で Unscented Kalman Filter (UKF)

こんなのに二日もかけて効率悪すぎた・・・

最近、RでUnscented Kalman Filter (UKF) のシミュレーションプログラムを書いていて、全然うまくいかなったんですけど、やっとこさ上手くいって嬉しかったのでブログに書きます。いきなりカルマンフィルターでしかもUnscentedってなんだよって感じですが、非線形で正規性を持つ状態空間モデルにカルマンフィルターを適用しようじゃないかということで考案されたものです。他に有名なものはExtended Kalman Filter (EKF) や Ensenble Kalman Filter (EnKF) があります。EKFはカルマンフィルターさえ理解すれば簡単にわかると思います（ざっくり言うと、非線形関数のテイラー展開を利用して線形近似してカルマンフィルターを適用するだけです）。EnKFはこれから勉強します。

一方UKFですが、これEKFとはかなり違うアイデアに基づいていて、ぶっちゃけカルマンフィルタという名はついているもののその面影はあまりないような気もします。

さて、ここで考える状態空間モデルは次のような式で表されるとします。

 

{ \displaystyle 　y_t = Z_t(\alpha_t) + \varepsilon_t, \hspace{22mm} \varepsilon_t\sim {\rm N}(0, H_t) \\ \alpha_{t+1} = T_t(\alpha_t) + R_t(\alpha_t) \eta_t, \hspace{10mm} \eta_t \sim {\rm N}(0, Q_t) }

 

ここで Z_t(\alpha_t), T_t(\alpha_t), R_t(\alpha_t) はそれぞれ非線形な関数としましょう。実際モデルを考える上では R_t(\alpha_t)は \eta_tをp次元ベクトルとした場合にp次の単位ベクトルとすることがほとんどだと思います。（そういうものしかみたことがありません。。。もしそうでなかったらごめんなさい）このようなモデルを考えた時のフィルターを考えたい。つまり、「一期前あるいはその時点までの観測値が得られた下での状態の条件付き分布の平均と分散共分散行列の推定」に興味があるんだよ！というのがモチベーションになります。フィルターを考える上でこれを忘れたり、心に止めておかないと何やってんだっけとなることが多いので個人的には大切なポイントだと思ってます。

Unscented Transformation (U変換)

U変換については、その本質的な意味合いをまだ消化しきれていないので、詳細は他の文献にお任せします（これから理解深めなきゃ・・・） 。参考文献からそっくりそのままもってくると、「正規分布に従う確率変数が非線形な変換によって写された先の分布が正規分布でないような場合を仮定し、元の分布、つまり正規分布からサンプリングした点を非線形変換して得られたサンプルを用いて平均と分散を近似する」という試みです。

U変換の手順は次のようなものになります。

{y \in {\rm R}^{p}, x \in {\rm R}^{m} }として、 { y = f(x) }となるような非線形な関数 {f} を考えます。このとき {x} は平均0ベクトル、共分散行列{P_{xx}}を持つp変量正規分布に従うとします。このとき、シグマポイントと呼ばれるサンプルを2m+1個要素として持つ集合 { x_0, x_1, \ldots, x_{2m} } とその重み { w_0, w_1, \ldots, w_{2m} } を次の条件を満足するように取ってきます。

  { \sum_{i=0}^{2m} w_i = 1, \hspace{5mm} \sum_{i=0}^{2m} w_i x_i = \bar{x}, \hspace{5mm} \sum_{i=0}^{2m}w_i( x_i-\bar{x} )( x_i - \bar{x} )^{T} = P_{xx} }

 このとき、 {f(x)} の分布の平均と分散を次のように近似することができます。 {y_i = f(x_i)} として、

{ \bar{y} = \sum_{i=0}^{2m} w_i y_i,  \hspace{5mm} P_{yy} = \sum_{i=0}^{2m} w_i (y_i - \bar{y})( y_i - \bar{y} )^{T}, \hspace{5mm} i = 0,\ldots, 2m }

このときシグマポイントは

{ x_0 = \bar{x}, \hspace{5mm} x_i = \bar{x} + \sqrt{m+k} \sqrt{ P_{xx,i} }, \hspace{5mm} x_{i+m} = \bar{x} - \sqrt{m+k} \sqrt{ P_{xx,i} } }

{ w_i = w_{i+m} = \frac{1-w_0}{2m} , \hspace{10mm} i=1,\ldots,m }

ここで \sqrt{m+k} の {k} は {m+k=3} を満たすように取ると良いでしょうとのことです。また { \sqrt{P_{xx,i}} } は { P_{xx} } の平方根行列のi列目のベクトルとします。（ここではコレスキー分解による平方根行列とします。）

 

最後にU変換をする関数を書いてみたのでそれを載っけて、近似がなかなか上手くいってることを示して終わりにしようと思います。

gistd25a7fa764618e82a814d3c7c35502f6

 

まずは２変量正規分布から乱数を発生させて、次の非線形な関数によって正規分布に従わないような確率ベクトルを生成します。

{\begin{pmatrix} y_1 \\ y_2 \end{pmatrix}  =  f(x) = \begin{pmatrix}  x_1 \cos{x_2} \\ x_1 \sin{x_2}  \end{pmatrix}   }

 このときのデータ点を黒で、平均と分散共分散行列に相当する楕円を青でそれぞれ描いて、赤い線でU変換による平均と分散共分散行列の近似を重ね書きしてみます。

f:id:forest-of-sheep-and-steel:20161012210924p:plain

この関数の場合、かなりの精度で近似できてます。この図のスクリプトは以下。

gistfdb330290067dcf0cb63e3a8e391cb2d

上手くいかない例なども気になるので今後の課題です。次回はU変換によるフィルターの構成とその近似精度がどれくらいなのか図示するところまで書きたいなぁと思っています。

[参考文献] 「カルマンフィルタの基礎」　足立修一・丸太一郎　共著

ターミナルでPATHの設定を間違えた時

MySQLを触れるようになろうと思ってインストールやらなんやら試していたら、PATHの設定を間違えてコマンドが全然通らなくなって焦ってて、やっと直せたので備忘録的にメモ。

gistf7f6d547661b601fb9bb88298ca71ce4

 

最初に設定されている基本的なPATHである↑を再設定すればOKとのこと。

いやぁ、焦った。

データ解析のための統計モデリング入門　２章 - ①

一回目から（中身のない一回目）時間空きすぎで自分でも結局やらないような気はしてましたが、やります。できるだけ、少しづつ。

今日はタイトルの通り久保拓弥さんの「データ解析のための統計モデリング入門」の２章の内容をまとめていきたいと思います（が、やってみて半分までしか進みませんでした）。１章は、気になる人は読んでみてください笑

 

この章のポイントは下記だと思ってます。

  • データの発生メカニズムが確率分布に従っているという仮定を理解すること

慣れない人にとっては確率分布って何、という所がもやもやする所だと思います。僕の考えですが、ざっくり言えば「起こりうる事象に対して、もれなく生起確率を割り当てるもの」だという理解でいいと思っています。もちろん、事象の確率の総和が１になるとか、確率が必ず０以上であるとか、厳密には色々と制約がついてきますが、そこら辺に興味がある方は数理統計の入門書など見ていただければと思います。

 

さて、２章では架空の植物から得られた種子数データを使って、このデータが、あるポアソン分布を使ってうまい具合に表現できるよねということを確かめていきます。この、うまい具合に表現できる、というのは、再現できる確率が高い、という意味で僕は理解をしています。

使われているデータはこんな感じでした。

gist98bbd458ba96f99215d5583a4cfc36ad

Rのコードを貼り付けたものですが、

生態学データ解析 - 本/データ解析のための統計モデリング入門、ここで本に扱われているデータをダウンロードできるので、ここから持ってきました。

ls() という関数は今R上で保存されているオブジェクトを全て列挙する関数です。立ち上げた時点では何も表示されないはずですが、load()関数でデータを読み込んだのでdataというオブジェクトが保存されていることが確認できます。よく今何があるんだっけ、とか、作りすぎて消したいけどどれを消していいんだっけ、見たいな時に僕は使います。

読み込んだらとりあえず要約値を見たり図示したりしてデータがどんなもんかを把握するのが、データ解析の第一歩だ、ということがよく言われているので、そうしてみましょう。

summary()関数は四分位点と平均値を返してくれて、table()関数は数値毎に何個あるのかを表形式で返してくれます。また、var()関数は分散を返してくれます。ぱっと見て、平均が3.56で平均に近い値が多く見られているといった感じですね。分散は約2.99でした。つまり、この中から１つデータをとってきた時に、平均すると、平均3.56と取り出したデータとの差の二乗が2.99です、ということです。なんで差の二乗なんだって所ですが、今の所は分散の計算をする時に「各データと平均との差の二乗」を使っているから、ということにしておきます。でもこれだとややこしくて、じゃあ一体実際の数値としてどのくらいなの？という話になるので、これの平方根を取って単位を直してあげます。これが標準偏差です。

標準偏差を使って言い直すと、このデータから１つ値をとってきた時に、平均すると、平均3.56との差が標準偏差1.73くらいあるでしょう。ということになります。大分すっきりしましたね。

ヒストグラムで図示もしてみます。

f:id:forest-of-sheep-and-steel:20160921210343j:plain

 数値的にみるのと違って、視覚的に見ると把握しやすいですよね←
本の図と若干違っていますが、それは引数breaksを指定していないので0と1がごっちゃになってるためです。指定しないとこうなるよーってことでこのまま載せました。
 
ここで、この章のポイントだと勝手に挙げた「データの発生メカニズムが確率分布に従っていると仮定している」ことを思い出して、今扱っている種子数のデータが従っている確率分布って何なんだろう、という問題について考えてみましょう。
  • 「種の数」なので、負の値はとらない（種の数が-3ということはありえない）
  • 同様に、1.4などの少数点以下に0以外の値を持たない。
ということに着目すると、この条件に合う確率分布を持ってきて当てはめてみようかな、という発想になるわけです。（こんなこと言われても、これって経験が重要ですよね。）
シミュレーション用のデータなのでこっから割とサクサクいくわけですが、この条件に合う確率分布といえば、そうだポアソン分布があるじゃないか！ということで、ポアソン分布を当てはめてみましょう。
とその前に、ポアソン分布って何なのという所ですが、定義の数式は今回は省いて、「パラメータλによって形が変わる確率分布」ということにして話を進めます。形が変わるとは一体・・・、という人のために次の図を描いてみました。（本にも同じような図があります。）

f:id:forest-of-sheep-and-steel:20160921222003j:plain

 

ポアソン分布のパラメータλを1,3,5,8とした時に0~10をとる確率をプロットした図です。「形が変わる」というのが視覚的にわかると思います。λの値が変わると線が描く形も変わってますね。縦に入っている点線はλの値と対応していて、ポアソン分布の場合はλの値が平均と分散と相当します。
余談ですが、最初は一番確率が高い点に点線を描いてたんですが、平均＝最頻値とは限らないということに気づかされました。学びました←。
λが取り得る数だけポアソン分布がある、つまり無数のポアソン分布の中から種子数のデータに一番合うポアソン分布を選ぼう、という話になってきたわけですが、今回は最後に見た感じ一番合いそうなポアソン分布を「視覚的」に探して終わりにしたいと思います。

f:id:forest-of-sheep-and-steel:20160921221909j:plain

gist65465013b1fdab935be94045ddba473c

ヒストグラムに重ね書きしてます。par(new=T)と入力すると、次に出力された図をそのまま重ね書きしてくれます。軸とか軸名とかも全部。なので後から出力される図はそれらを全部消す形で書いてます。こういう重ね書きの仕方は本当は良くなくて、できるならpoints()とかlines()で書いたほうがいいんですけど、縦軸の値を合わせたくなかったので力技で書いてます。ちなみに、各ポアソン分布については縦軸を合わせて書いてます。

これをみると、図の中で、データの分布と同じ形をしてそうなのはλ=3.56かな〜というあたりがつきますね。（出来レース感がすごい笑）

ということで今回はここまでです。次は「見た目」、じゃなくて「数値的」に最適なポアソン分布を探そうという所からです。

 

ブログって書くの疲れますね。内容や文章など拙いもので申し訳有りません。何かご指摘があれば是非いただければと思っています！
ここまで読んでくださった方がいたら、本当に有難うございました。これからも細々と書いていこうと思います。では。m(_ _)m

統計学を勉強していく様を見せるブログスタート

ブログを始める理由

こんにちは。統計学を専攻しているしがない学生です。勉強したことは発信して、炎上させて理解を深めることが重要だということに気づいたので、ならブログをやろうと思い立ち始めました◀️

直近の目標

次の二つのテーマについてちょこちょこ書いていこうと思っています。

１、「データ解析のための統計モデリング入門」(久保拓弥, 岩波)

２、状態空間モデル,   参考文献は色々

あとはその都度、備忘録的に勉強したことを書いていければなと。

 

モチベーション

理由にも書いた通り、勉強したことは発信してボコボコにされてこそ身につくものだと思っているので、そういう場をネット上にも設けた方がいいと思った次第です。

 

使用するソフト

基本的にはRを使っていきます。できるだけコードも綺麗に書いて載せたい。けどプログラミングのイロハはよくわからないので見辛いものになるかもしれません。

おそらくRの使い方自体の話も、特に覚えておきたいことがあればそれをテーマに記事を書くこともあるでしょう。

 

そもそもな話

私自身の話をすると、一応院で勉強してるくせに、数学が得意ではない、要領が悪い、集中力がない、と三拍子そろったダメ学生と自覚してます。これは何か変えなきゃいかんということでたまたま友達が始めたブログの話を聞いて、これだ！と思ったんですよね。統計学のブログって色々あるので正しさとか厳密さとか読みやすさというのは実力のある方々にお任せして（おい）、のほほんと間違いながらも自分がやったことを残していく場として活用したいっていうのが本音です。

もし、どこかでどなたかのお役に立つことがあれば嬉しいことこの上ないのですが、正直そこは今の目標ではありません！

 

次回より

次の記事から本格的に中身を作ってアップしていこうと思います。どんな人が見るのかなーとか、需要あんのか、とかすごい不安ですけど、きままに頑張ります。