2017-06-12

R で Nelder and Mead 法を実装してみる（勉強のため）

久しぶりの更新です。事情によりNelder and Mead法を実装する必要が出たのでやってみました。参考文献はNelder and Mead(1965)と

ネルダーミード法はRのoptim関数で使われていたり、調べたらRのパッケージとしても公開されていたり、別に自分で実装しなくてもいい人はたくさんいると思いますが、私はそういうわけにも行かなかったのでやりました。そしてついでにブログに書きます。

この方法の一番の特徴は導関数を求めなくてもいいという点ですね。実行するために必要なのは

関数形： $y=f(x)$
初期値： $x_{init}$
パラメータ $\alpha, \gamma, \rho, \sigma$
判定量： $criterion$

になります。今回のシミュレーションでは論文のResultの1個目としました。（下記）

$f(x) = 100(x_2 - x_1^2)^2 + (1 - x_1)^2$
$x_{init} = (-1.2, 1)$
$\alpha = 1, \gamma=2, \rho = 0.5, \sigma=0.5$
$criterion = \sqrt( \sum_{i=1}^{n+1} ( f(x_i) - \bar{ f(x_i) } )^2 )$

関数定義とサンプルコードも載せておきます。でも初期値依存するので場合によっては正解に収束していないケースも・・・。

Nelder and Mead 1965 for R

Sample Nelder and Mead 1965 1 for R

2016-12-21

Rで日本地図を書く（たぶん最終回）

一体Rで何をしたいんだというのは置いといて、自分の作業内容でRで日本地図を柔軟に書きたいなぁと思ったので、これまで書いてきたコードをちょこっと改良しました。多分今回作ったのでだいたいやりたいことはできるはずなので最終回のはず。

一応本の方もちゃんと読み進めています。一応。

今回ので何ができるようになったのかというと、

指定した都道府県のみをかけるようにした
描画範囲を自動と任意とで使い分けができるようにした

の２点です。

まずは、以下のコードで世界地図がのデータから日本のデータを引っ張ってきます。この時に全部描く用と県毎に描く用のデータを分けて作っています。データについては前回か前々回のブログを参照ください。

gist5ec1f0eec3d282aa9db3bcb9892d486e

これでmap.dataとmap.Allというリストができました。これらを使ってプロットする関数を次のように2Dと3Dに分けてそれぞれ定義します。

gist58eb3c17bb62660d01a8bbdea6496c74

gist02d8484b4d91bf698fdb7c692d6c0131

すごい手抜きなブログですが、内容は以上です。

こんな感じでプロットできます。

１つ目は3dプロットの全国地図

もう１つのは2dプロットの熊本、神奈川、埼玉、青森、新潟、岡山、京都、岐阜の地図です。

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

年内にもう一回更新したいなと思ってます。普通に統計の話で・・・。

以上でーす。

更新：ちなみに、shapeファイルを扱うRのパッケージもあるので、こんなめんどくさいことは実はやらなくてもいいのです笑

2016-11-01

Rで日本地図を書く（３次元で）②

よくよく見たら、この前書いた日本地図、県境がかなり適当。細々した島も無い。特に沖縄が適当。ってことで、もっといい地図はないのものか探してきました。んで、見つけたのがこのデータ！

www.naturalearthdata.com

今回は上から５グループ目の

Admin 1 – States, Provinces

という見出しの Download as scale ranks というデータを使いました。

さっそくコードを。。

gist025531edb17547aa2803766ad18887a4

このデータは県境もあるし、結構細かいし、なかなか良いと思う。二次元にプロットするとこんな感じ。

gist57b20e6b5ba46b69e03f6fed1b45a176

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

３次元で写したい場合はこちら。

gist283c11451c49f1949bbd652d1710bfa6

結果はこんな感じ。いい感じ。。

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

そもそも何がやりたかったかというと、天気予報とかでよくある降水量のバーがついた地図を作りたかったのです。降水量のデータはちょっと載せられないのでサンプルデータで書いた結果だけ。やってみればそんなに苦労せずできますね！

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

そういえば最近、PVが20を超えました笑

2016-10-30

全く関係ない話２（対数関数の性質）

高校で初めて対数関数を習ったとき、それを使うメリットがわからなかった。今は理解しているつもりだけど、こんな説明があったら理解の助けになるのかなと思って作図。またまたx軸, y軸を平行にしてどのように値が変化されているのか確認する。普通の二次元プロットでもわかる人はわかるんだろうけど・・・。

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

gisteef85984474dbb1c667ca76f1eb0bed9

上の横軸がxの値、下の横軸がyの値ですね。ちなみに、この時対数の底の値はネイピア数eとしています。適当にxの値を２つ持ってきたときに、それを対数変換によってy軸に移しても、それぞれのx軸での大小関係は保たれたままy軸に移ってる事がわかると思います。

これが結構大切で、大小関係を評価したいとき、対数を取っても問題ないというのはこの性質からきているんですね。

もう一点は、対数関数の場合、xは0より大きいという制約がありますね。それを踏まえて、この対数関数は、すごく小さい数（1より小さく0より大きい）は負のより大きい数字に、すごく大きい数は正のより小さい数字にしてくれるという性質があります。

0.001 -> -6.907, 1000 -> 6.907,

のように、0.001と1000ではその差が999.9もあったのに、２つの対数をとると、大小関係を保ったまま、-6.907, 6.907とその差が13.815となっています。極端に大きい数字とか、極端に小さい値を扱う時、対数をとると普段使っているようなスケールの数字に変換してくれるんですね。しかも大小関係は保存してくれるという嬉しさです。

とりあえず以上。。。

2016-10-26

Rで日本地図を書く（３次元で）

かなりニッチな需要だと思いつつ、Rで３次元の中の平面に日本地図を書こう、ということをやっていきます。目的はと言いますと、天気予報なんかである、日本地図の対応する地点に観測された降雨量の棒グラフを書きたい、みたいなことがあって、今回はじゃあ３次元の中に平面の日本地図を書くというところまでやっていきます。

必要なパッケージは３つ、Nippon, spsurvey, rgl の２つです。とりあえず読み込みましょう

gistfa75446b9c7707ba6829d351f02361e9

Nipponは他の方のブログでもよく紹介されているように、日本地図を書いてくれる関数がありますが、この記事は、いやいや俺は座標の行列から日本地図を書きたいんだよ！っていう人向けでもありますね。

Nipponのパッケージを使えば地図がかける、つまり座標はこのパッケージの中にあるわけです。まずはそれを抽出していきましょう。そのために、Nipponパッケージがあるフォルダーにアクセスする必要があります。windowsのことはよくわからないのでごめんなさい。macの方は

/Library/Frameworks/R.framework/Versions/3.3/Resources/library/Nippon/shapes

ここにshpファイルがあって、その中にプロットに必要な情報が色々と入っています。

では、座標データだけ抽出していきましょう。次のコードを実行すればOKです。

giste07d64552229c7dbbaa125891e151fd4

これで座標を取り出せました。例えばこれを使って北海道をプロットするとこんな感じになります。

gist83ccd26a79fe22632aa844d6ece88e21

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

あとは力技でちょっとカッコ悪いんですが、3次元平面にプロットしていくだけです。(三軸目の値は全て0にしています。)

gistde0302bc23bbbbb84dc95b42c60f99b3

どうしても、色々いじりたい場合はパッケージの関数一発で書くよりも、座標データだけ欲しい状況が出てくるわけです。こんな感じのニッチな需要を持っている人の役に立てばいいなーと思います。

結果はこんな感じでした。でもこの地図、沖縄がかなり省略されてるみたいですね。

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

2016-10-12

R で Unscented Kalman Filter (UKF)

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

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

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

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

$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 }$