Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2015-01-30

二次元データの(線形)補間を行う

| 21:39 | 二次元データの(線形)補間を行うを含むブックマーク

一次元のケースは二年も前に書いた

でよいとして、二次元のケースをなんとかしたい。したいというか、しなければならない状況になった。

ググりまくったところ

というパッケージに入っている関数を使えばよさげだということがわかった。


この2つのパッケージは、手元の補間したいデータの構造に応じて使い分ければよさそうだった。具体的には、以下のように適当にデータを作っておく。

library(dplyr)
library(tidyr)
#適当なサンプルデータ
X0 <- 1:10
Y0 <- 1:10
df <- expand.grid(x=X0, y=Y0) %>% mutate(z=x^2-y^2+y)
ls <- list(x=X0, y=Y0, z=matrix(df$z, nrow=length(X0), ncol=length(Y0)))

これは、dfがデータフレームで、lsがリストのデータ構造になっている、データ形式だけが違う同一のデータだ。中身の意味は見ればだいたいわかるように、x, yという二次元座標値に値zが割り振られている、そういうもんだ。

> df %>% head
  x y  z
1 1 1  1
2 2 1  4
3 3 1  9
4 4 1 16
5 5 1 25
6 6 1 36
> ls %>% head
$x
 [1]  1  2  3  4  5  6  7  8  9 10
$y
 [1]  1  2  3  4  5  6  7  8  9 10
$z
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    1   -1   -5  -11  -19  -29  -41  -55  -71   -89
 [2,]    4    2   -2   -8  -16  -26  -38  -52  -68   -86
(... 途中省略 ...)
 [9,]   81   79   75   69   61   51   39   25    9    -9
[10,]  100   98   94   88   80   70   58   44   28    10

このような2つの異なるデータ形式の同一のデータがあり、それらの二次元補間をしようと思った場合には、

  • データフレーム(data.frame)の場合はakimaパッケージ
  • リスト(list)のの場合はfieldsパッケージ

という使い分けで補間してやるのがよさそうだった。実際は、データ形式を揃えればどちらか片方のパッケージだけでいいんだけども。rlistパッケージあたりを使えば簡単にデータの変換できるのか?


実際にたしかめるために、まず、補間する点(x, y)の組を以下のように与えておく。そして、各点に対する線形な補間値を計算させてみる。また、あわせて答えの一致も確認する。

#補間する点
x <- c(2, 4, 8)
y <- c(3, 7, 2)

akimaパッケージの場合

akimaパッケージは不規則グリッド*1上にあるデータのための補間パッケージだ。このパッケージで二次元線形補間をするためには

を用いる。interpp関数の方が"点ごと"の補間に対応して、interp関数は、下記コードのコメントに書いたようにx×yの直積集合の全ての元に対して補間を実行してくれる。

library("akima")
#点(2,3)(4,7)(8,2)での線形補間
inter_pp <- interpp(df$x, df$y, df$z, x, y)
#直積集合(2,4,8)×(3,7,2)での線形補間
inter_p  <- interp( df$x, df$y, df$z, x, y)

結果はリストで帰ってきて、zに補間値が入っている。

> inter_pp$z
[1]  -2 -26  62
> inter_p$z
     [,1] [,2] [,3]
[1,]   -2  -38    2
[2,]   10  -26   14
[3,]   58   22   62

ちなみにx, yには補間した位置情報が入ってる。また、interp/interpp関数ともにlinear引数があって、これをFALSE(デフォルトはTRUE)にするとスプライン補間にかわるはずなんだけど、このデータだとなぜか全部NAになった、なぜ???関数を掘っていったら、"sdsf3p"なるFortran関数がよばれていたので、これを読まないとだめだな・・・

fieldsパッケージの場合

fieldsパッケージは空間統計のためのツールが格納されたパッケージとのことだ。このパッケージでは

という2つの関数が2次元線形補間を実行してくれる関数だ。これもakimaパッケージのinterp & interppの関係に近くて、点の補間かグリッドまるごと(直積集合)の補間かで関数が異なっているだけ。

> library(fields)
> #点(2,3)(4,7)(8,2)での線形補間
> interp.surface(ls, cbind(x, y))
[1]  -2 -26  62
> #直積集合(2,4,8)×(3,7,2)での線形補間
> interp.surface.grid(ls, list(x=x, y=y))
$x
[1] 2 4 8

$y
[1] 3 7 2

$z
     [,1] [,2] [,3]
[1,]   -2  -38    2
[2,]   10  -26   14
[3,]   58   22   62

結果はakimaパッケージのものとちゃんと一致している。どちらのパッケージを使うのかは、データ形式に応じて使い分けるか、変換の手間を惜しまず、どちらか一方だけ使っててもいい感じかな。

rvestで単元株数を取得したい

| 06:35 | rvestで単元株数を取得したいを含むブックマーク

やほぉの単元株数、こんなんで取得OKでござんした。

を見て、適当にcodeを変えれば他社でも動く(はず)。

以下のコードを実行すると"100"となるので、単元株数は100だなとわかる。

library(rvest)
library(stringr)
html <- html("http://m.finance.yahoo.co.jp/stock/additionaldata?code=4689.T")
node <- html %>% html_node("#additionaldata") %>% html_text
unit_stock <- str_match(node, "単元株数.*株\\n") %>% c
str_match(unit_stock, "\\d+") %>% as.numeric

参考

*1:均一なグリッドじゃないという意味かと

2015-01-28

日本国債(JGB, Japanese Government Bond)ちゃんが息してない

| 07:07 | 日本国債(JGB, Japanese Government Bond)ちゃんが息してないを含むブックマーク

「低金利ここに極まり」みたいな水準の国債市場、みなさんいかがお過ごしですか?

どのくらいキテるのかを見るために、直近1年程度の10年ものの日本国債(JGB, Japanese Government Bond)の推移を可視化してみました。

2015-01-12

リカレント・ニューラル・ネットワーク(Recurrent Neural Network, RNN)っていう、時系列に対するDeep Learningやりたい

| 21:20 | リカレント・ニューラル・ネットワーク(Recurrent Neural Network, RNN)っていう、時系列に対するDeep Learningやりたいを含むブックマーク

はじめに

前々からDeep Learning系、特にその時系列への応用に興味があって、それにはリカレント・ニューラル・ネットワーク(Recurrent Neural Network, 以下RNN)ってのを理解する必要があるのは知っていたが、@yamano氏から

というありがたいアドバイスを頂戴した。それに基づいてググりはじめると、Stuttgart(シュトゥットガルト)大学というドイツにある、たぶん日本でいう東工大的なポジションの大学が開発している、ニューラルネット関連の実装をたくさん詰め込んだライブラリ

が元々あって、それをラップしてRから叩けるようにしたものが、

というライブラリということがわかった。

RSNNSパッケージを使ってみる

上に書いたように、このパッケージ(ライブラリ)に入っている実装は、RNNに限らないようだが、興味あるのはRNNなのでそれをやる。Rのパッケージにくっついてくるマニュアルを読み込むとどうやら「jordan」関数ってのがあって、これを使うとよさそうだ。似たようなものとしてElman networkというのもあるらしい。詳細は原著

  • Jordan, M. I. (1986), ’Serial Order: A Parallel, Distributed Processing Approach’, Advances in Connectionist Theory Speech 121(ICS-8604), 471-495.
  • Elman, J. L. (1990), ’Finding structure in time’, Cognitive Science 14(2), 179–211.

を漁るか、Stuttgart Neural Network Simulatorのマニュアルを読まないとだめだな。どちらも「partially recurrent networks」と記述されていて、「なんでpartially?」ってのが気になってる。まぁ、とりあえず

  • 興味持つ⇒動かす⇒ロジックを調べる⇒より興味を持つ⇒より興味を持って動かす⇒・・・

ってことで、パッケージに入っている例のコードを少し改変して動かしてみる。

#パッケージ入れる
install.packages("RSNNS")
library(RSNNS)
library(dplyr)
#データ読み出す
data(snnsData)
#snnsDataデータを入力と出力で分解
inputs  <- snnsData$laser_1000.pat[,1]
outputs <- snnsData$laser_1000.pat[,2]
#更にトレーニング・
patterns <- splitForTrainingAndTest(inputs, outputs, ratio=0.15)
#size:隠れ層のユニット数
#learnFuncParams:学習関数に対するパラメーター
#linOut:活性化関数をリニア(TRUE)にするかロジスティックにするか(FALSE)
modelJordan <- patterns %>% with({
  jordan(inputsTrain, targetsTrain,
    size=c(8), learnFuncParams=c(0.1), maxit=100,
    inputsTest=inputsTest, targetsTest=targetsTest, linOut=FALSE)})

ここで、デフォルトの学習関数

によると、

  • JE_BP: Standard Backpropagation for partial recurrent networks

となっているが、これが単純なバックプロパゲーションでやってますよーという意味だと思われる。この辺を理解するにもRのパッケージじゃなくて直接Stuttgart Neural Network Simulatorのマニュアルを読まないとだめだな。で、この学習に使われるパラメーターがlearnFuncParamsで設定していると。

また、ここで使っている元々のデータは「カオス状態にあるアンモニア(NH3)レーザーの強度のデータ」だそうな。これは、公式サイトからソースを落としてきて、その中に入ってる「laser.README」というファイルに記述がある。

計算結果(modelJordan変数)を描画する用のplot関数も入っているが、中のコードを見た感じ通常のPLOT関数+αな印象。また、フィッティング結果が結構上方にバイアスかかってるように見えるんだが、こんなもんなのだろうか?

#適当に結果を描画
plotIterativeError(modelJordan) #1
plotRegressionError(patterns$targetsTrain, modelJordan$fitted.values) #2
plotRegressionError(patterns$targetsTest, modelJordan$fittedTestValues) #3
plot(inputs[1:100], type="l") #4
lines(outputs[1:100], col="red")
lines(modelJordan$fitted.values[1:100], col="green")

#1の結果(繰り返し毎に加重残差二乗和が減っていきますよー的なもの)

f:id:teramonagi:20150112210734p:image

#2の結果(トレーニングデータ(実データ) v.s. その予測値。黒線は45度線(100%HITの場合これ)、赤線は(予測~実データの単回帰直線))

f:id:teramonagi:20150112210735p:image

#3の結果(テストデータ(実データ) v.s. その予測値。黒線は45度線(100%HITの場合これ)、赤線は(予測~実データの単回帰直線))

f:id:teramonagi:20150112210736p:image

#4の結果(黒線:入力値、赤線:実データ(出力)、緑線:予測値)

f:id:teramonagi:20150112210737p:image

2015-01-11

やまかつ氏の滑り続ける<del>リプライへの統計的考察、及びその実装 -Deedleを添えて- </del>タイムライン取得

| 15:59 | やまかつ氏の滑り続ける<del>リプライへの統計的考察、及びその実装 -Deedleを添えて- </del>タイムライン取得を含むブックマーク

はじめに

こういう話があった。

これをF#を使って確かめようといろいろ調べていたら、とりあえず

を使えばいいだろうなってのを理解したんだけど、メンションを取得するAPIである

に対応するものが無かったので、頑張って開発した結果、以下のようにMerge待ち状態なんだけど

これよくよく見ると「認証ユーザーへのつぶやきを示す」メンションしかとれなくて、俺ではやまかつ氏のメンションは取れねぇよ!ってことであきらめた。

しょうがないのでタイムラインだけ取ってみた

以下のコードで、やまかつ氏の呟きが取得可能である。

への参照は張っておくこと。また、キーは

にあるから適当に使わせてもらってもいいのかもね。

//必要な名前空間のオープン
open System
open FSharp.Data
open FSharp.Data.Toolbox.Twitter
open FSharp.Data.JsonExtensions
//接続
let key = "Your Key"
let secret = "Your Secret"
let twitter = Twitter.AuthenticateAppOnly(key, secret)
//呟きの取得(直近5件)
let yamakatu = twitter.Timelines.Timeline("yamakatu", 5)
yamakatu |> Array.iter (fun y -> printfn "%s" y.Text)

結果は↓。

あ、今日はブリテロの日か
@teramonagi @aad34210 みのださんが元気かどうかは知らないけど、俺なら元気だよ―っ!
RT @dichika: 「MikonはRubyにおけるpandasのような位置づけで」 http://t.co/pTJEskBTnL #メモン
おはようございます! #パクった
RT @debugordie: じゃばにゃん というのが姉から送られてきたw http://t.co/ZowTb1goNE
val it : unit = ()

2015-01-10

クロージャー(Closure, 閉包):2.0 YOU CAN (NOT) ADVANCE

| 20:20 | クロージャー(Closure, 閉包):2.0 YOU CAN (NOT) ADVANCEを含むブックマーク

結論とまとめ

クロージャー(Closure, 閉包)について、1年経っても何も学習しておらず、また同じ罠にハマった。ググラビリティのための再掲。一年前あたりの記事はこの辺。

どうしたらいいのかって話は、以下の@tyatsuta氏の資料をちゃんと読めってことだ。

ハマり所の内容

以下のような

を作る。

make_add <- function(val)
{
  function(x){
    x + val
  }
} 

使用例はこんな感じだ。ただし、これはsapplyにかませるとうまくいかない。

> add_three <- make_add(3)
> add_three(1)
[1] 4
> add_five <- make_add(5)
> add_five(1)
[1] 6
> 
> add_list <- sapply(1:3, make_add)
> add_list[[1]](1)
[1] 4
> add_list[[2]](1)
[1] 4
> add_list[[3]](1)
[1] 4

これを直すには一度関数内で強制的に評価するようforce関数をかませると良い。あるいは、要するに強制的に評価させたいだけなんで、単にvalとだけ書いても良い。force関数もやってることは同じで、単に引数返しているだけ。

make_add <- function(val)
{
  force(val)
  function(x){
    x + val
  }
} 

これでうまくいく。

> add_three <- make_add(3)
> add_three(1)
[1] 4
> add_five <- make_add(5)
> add_five(1)
[1] 6
> 
> add_list <- sapply(1:3, make_add)
> add_list[[1]](1)
[1] 2
> add_list[[2]](1)
[1] 3
> add_list[[3]](1)
[1] 4

別な対策とメモリ使用量

「forceって書くよりも、もう別な変数に突っ込んでおいたらいんじゃね?」ってやってみた。元の関数と区別するために適当な別名つけて、元の関数関数名改める。

make_add1 <- function(val)
{
  .val <- val
  function(x){
    x + .val
  }
} 
make_add2 <- function(val)
{
  force(val)
  function(x){
    x + .val
  }
} 

結果はどちらでもOK

> sapply(1:3, make_add1)[[2]](1)
[1] 3
> sapply(1:3, make_add2)[[2]](1)
[1] 3

さらに、object.size関数でこのクロージャーのサイズを調べてみると

> object.size(make_add1)
5488 bytes
> object.size(make_add2)
5432 bytes

やはり余計な変数作ってる分、ちょっと重くなる感じか。ただ微々たるものなので、可読性を優先してやるってのもありだな。