Hatena::ブログ(Diary)

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

2014-11-28

二次元カーネル密度推定

| 09:46 | 二次元カーネル密度推定を含むブックマーク

Rで二次元カーネル密度推定してみる。別に2次元が好きなわけではなく、これはたまたまだ。1次元の場合は、Rにデフォルトで組み込まれているstatsパッケージのdensity関数を使えばヨロシ。

デフォルトで入ってくるパッケージの調べ方とインストール

MASSパッケージにあるkde2d関数を使いたいので、このパッケージを入れるんだが、MASSってデフォルトで入っていたような気もする。気もするが、どうだったかは覚えていない。こんな時、実際に入っているかいないかを調べてみる。RでロードされるデフォルトのパッケージはgetOption関数にdefaultPackagesを文字列として与えると取得できる。

> getOption("defaultPackages")
[1] "datasets"  "utils"     "grDevices" "graphics"  "stats"     "methods"  

あぁ、入ってない。というわけで、install.packages関数インストールする必要がある。そしてMASSパッケージもロードだ。

install.packages("MASS")
library(MASS)

Rでやる

適当に2次元混合正規分布を作ってカーネル密度推定してみる。pipeRパッケージをちゃんと使いこなせるようになるのがいい大人の条件だと強く信じて疑わないので、pipeRしてます。

library(MASS)
library(dplyr)
library(pipeR)

#データ数
N <- 10^3
#中心位置
mu <- list(mu1=c(1,1), mu2=c(5, 8))
#データ生成
x <- lapply(mu, function(mu)data.frame(x=rnorm(N, mu[1]), y=rnorm(N, mu[2]))) %>>% rbind_all
plot(x)

#二次元カーネル密度推定
x %>>% with({
  #バンド幅の計算
  h <- c(bandwidth.nrd(x),bandwidth.nrd(y))
  #密度推定
  kde2d(x,y,h,n=80)
}) %>>%
(~ res) %>>% #結果をresに格納
image

kde2d関数で密度推定した結果をres変数に突っ込みつつ、そのまま画像出力(image関数)させているわけですpipeR素晴らしい。

f:id:teramonagi:20141128094542p:image

ちなみに、kde2d関数で使えるカーネルはガウシアンカーネル限定だ。その他のパッケージについては

PDFを見て、どんなのがあるのか調べていくのがよさげ。気がついたらこの資料の著者の1人も神Hadleyでびっくりだ。そしてこのPDF一次元の方法(univariate)のみに特化してるそうで、二次元の場合はまた別なソースを探す必要がある。。。どこ。。。

ksパッケージなるものがあって、これを使えばよさげか?

このパッケージだと1〜6次元まで対応とのこと。ただ、マニュアルを読む限り、こいつもガウシアンカーネル限定だそうだ。そういうもんなのか?

バンド幅の選択など

kde2d関数引数でいうところのhで指定するのがカーネルのバンド幅だ。こいつの選択方法・その良し悪しは、あまりよくわかってないが、必要に応じて以下のサイトの日本語役と手法名を照らし合わせて、テキスト/原著論文読んで〜ってやらないとだめだろうな。

上のコードで使ったバンド幅計算関数であるbandwidth.nrd関数はHELPで確認した限り

  • Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Springer, equation (5.5) on page 130.

に詳細があるそうで、その実装も簡単でHELPに記載があり、データの25% & 75%クオンタイル点の差を1.34で除した量をむにゃむにゃルート取ったり計算しているだけだった、なんだこれ。

ちゃんと確率密度関数になってる?

ここでやってるカーネル密度推定は、データを¥left¥{ x_i; i=1,¥dots,N¥right¥}として、確率密度関数p(x)を、元のデータが属している標本空間¥Omega上に定義された適当なカーネル関数K : ¥Omega ¥times ¥Omega ¥rightarrow ¥mathbb{R}を使って

  •  p(x) = ¥frac{1}{N} ¥sum_{i=1}^{N} k(x, x_i)

として表現する手法なわけだ。従って、計算結果(res変数)がちゃんと確率密度関数になっているかどうかは、xについて¥mathbb{R}、あるいは高次元(D次元)の場合は¥mathbb{R}^Dの範囲で積分して、その答えが1になるかを確認すればよい。ここでもpipeRを使う。単に結果を台形公式で数値積分しているだけですわ。

> #確率が大体1になるか?
> res %>>% 
+   (~ dx <- . %>>% (x) %>>% diff %>>% first) %>>%
+   (~ dy <- . %>>% (y) %>>% diff %>>% first) %>>%
+   (z) %>>% 
+   sum %>>% `*`(dx*dy)
[1] 0.9928492

荒々な計算で大体1になってるので、多分大丈夫だろう。

こいつって正定値性は担保されてるんでしたっけ?

機械学習の文脈でいうカーネル法カーネル関数正定値性を持たせて、再生核ヒルベルト空間を作らせて〜のノリでいくわけですが、昔からある、ここでやってるようなカーネル法では正定値性を特に意識はしてないみたい*1

参考

↑二番目の資料の「用語に関する注意」が脳の混乱を防ぐために大事。

*1:参考LINKの2番目、福水先生の資料より、パルツェン窓関数〜って話がここでやってる話に相当

2014-11-27

IfSharp (notebook)サイコオオオオオオオオオオオオオオオオオオオオオオオオオ

| 00:11 | IfSharp (notebook)サイコオオオオオオオオオオオオオオオオオオオオオオオオオを含むブックマーク

はじめに

IPython notebook、昔々触ったときは「コンセプトは良いけどまだまだですナ。ハハハ」みたいな感じで流してたのだが、今ここでちゃんと触ってみようとおもったけれども、最近scikit-learnを必要に応じて触るくらいでpythonあまりいじってないし、俺にはF#あるし、iPython NotebookのF#版であるifSharp使った方が楽しそうだったので、未来の自分が別PCに導入することも考えてメモしておく。なんだかよくわからないけど、64bitマシンじゃないとだめだってさ。

とりあえず、Visual Studioが入っていないお友達は

からVisual Studio Community版いれとけ。

ifSharpの公式

ifSharpの公式は以下。ifSharp NotebookじゃなくてifSharpが正式名称っぽい。

Documentation

Github

インストール&簡単な例

公式(Github)のが間違ってる気がするのでメモる。上記のサイトで言う「Automatic Installation」でいこう、マニュアルインストールはめんどい。まずPython周りは大人しくAnacondaいれておく。

次いで、以下から最新のifSharpのsetup.exeを落としてYESYESしてインストール

これでインストール自体は完了で、あとはインストールしたifsharp.exeを起動する。デスクトップショートカット作ってればそれダブルクリックでOKだ。そうすると以下のようにiPythonが立ち上がる。そして、下図の赤丸で囲ってある「New NoteBook」を押下。

f:id:teramonagi:20141127234252p:image

そうするとIP[y]: NoteBookのIPの部分がIFへと変わり、ifSharpとしてNotebookを作成することができる。あとはチャッチャとF#のコードを書けばいい。

f:id:teramonagi:20141127234253p:image

出来上がったipynbファイルをGithubにアップすれば、iPython Notebook同様、nbviewerで見ることができる。

そもそも、iPython Notebookの使い方がよくわかってないので、そちらも調べる必要があるな・・・

その他特筆事項

rvestでリアルタイムな為替レートを取得したい

| 19:04 | rvestでリアルタイムな為替レートを取得したいを含むブックマーク

USD/JPYレート、こんなんで取得OKでござんした。html関数以下を再実行すれば直近のレートに置き換わる。

#パッケージ読み込み&html取得
library(rvest)
html <- html("http://info.finance.yahoo.co.jp/fx/detail/?code=USDJPY=FX")
#BID, ASK取得
bid <- html %>% html_node('#USDJPY_detail_bid') %>% html_text() %>% as.numeric()
ask <- html %>% html_node('#USDJPY_detail_ask') %>% html_text() %>% as.numeric()
#出力
print(c(BID=bid, ASK=ask))

このrvestライブラリの使い方は、もっと勉強する必要がある。

参考

RubyからRを叩きたい

| 15:00 | RubyからRを叩きたいを含むブックマーク

はじめに

F#に魂を売り渡したとしていても、そういうRuby欲、あると思います。ざっと調べたところ、RinRubyなるライブラリを使うのがよさそうだった。ただ2年前の情報なので、これが一番いいのかは自信あまりない。参考リンクにあるように、愛ちゃんパイセン(@dritoshi)が4年も前に試してるところをみると、今は他の何かがあるかも?

RinRubyのサイト

公式サイト

Github

このGithubを見たところ、別のRubyからRを使うライブラリとして、RSRubyなるものもある

こちらは日本語の記事がRinRubyよりひっかかるので、そっちでもいいかも?・・・なんて思ってたら、こっちは@berobero11氏がすでにきれいにまとめてくださっていた、ありがたい。

使い方

まず、Ruby側でのパッケージのインストール

gem install rinruby

インストール。適当に走らせたサンプルは以下。実行すると、正規乱数の結果を表示したり、ランダムフォレストの結果(元はテーブル、ただしRuby上ではRubyの行列)が出る。

# -*- coding: utf-8 -*-
require "rinruby"      
# R.xxxxでR上の変数にアクセス可能。ない場合は変数の作成
# R上で10個正規乱数を生成して、R上のxという変数に格納
R.x = R.pull "rnorm(3)"
# あるいはevalを用いて以下のようにもかける
R.eval "y <- rnorm(3)"
# 作った乱数配列の和・標準偏差を計算
x_sum = R.pull "sum(x)"
x_sd = R.pull "sd(x)"
# Ruby で出力する
puts x_sum
puts x_sd
puts (R.pull "y")
puts (R.y) #1つ上と等価な表現

# 以下のdata.frameの受け渡しはだめだった
#z = R.pull "iris"
#puts z

# ヒアドキュメントを使ってそれっぽい分析も(以下の例はランダムフォレスト)
# ライブラリの読み込みも可能
R.eval <<EOS
library(randomForest)
train <- sample(nrow(iris), 0.8*nrow(iris))
test  <- setdiff(1:nrow(iris), train)
iris.rf   <- randomForest(Species~., data=iris[train,])
iris.pred <- predict(iris.rf, iris[test,])
res = table(iris[test, 5], iris.pred)
EOS
puts R.res
# Rを終了
R.quit

data.frame渡したい

RinRubyでdata.frameをRubyとやりとりしようとすると難しいくさい。

なので、上のStackoverflowの回答では、またまた、別のライブラリとして

使うのがいんじゃね?って話が出ていた。このGithubにはここで挙げた3つのRubyからRを使うライブラリのプロコンが載っているので、それを見て、自分のやりたいことと照らし合わせるのがよさげっぽ。あるいは上記の@berobero11氏のBLOGにあるようにRSRubyかなあ。

参考

未来の俺へのメモ

インストールはRのバージョン指定したりでめんどいけど、大人しく@berobero11氏の

にならって、RSRubyが速度・使える型の多さからよさげ。

*1:Math Jaxか?

2014-11-13

円内に一様分布する乱数を生成する時は、俺、絶対忘れないよヤコビアンのこと

| 09:42 | 円内に一様分布する乱数を生成する時は、俺、絶対忘れないよヤコビアンのことを含むブックマーク

頭出し

ある半径Rの円内に一様分布する乱数を生成する時には注意しないといけないことがありますよというお話。所謂「一度はやってしまうミス」系でもある。この手の話は円に限ったわけではなく、円の高次元版である球、あるいは超球(次元>3)、あるいは任意の座標変換をかませてそこにヤコビアンが出て来るときでも同じ。

使うライブラリは以下の2つ。なければinstall.packages関数インストールしておく。

library(ggplot2)
library(dplyr)

また、以下のように定数を2つ定義しておく。意味はコメントにある通りだ。

#サンプル数
N <- 10^4
#半径のサイズ
R <- 4

本題

さて、問題のある半径Rの円内に一様分布する乱数を生成するにはどうしたらいいのかというと、非常に単純に考えた場合、以下のように(俺は)思考する。

  • X方向の成分として[-R, R]の間の一様乱数を生成する
  • Y方向の成分として[-R, R]の間の一様乱数を生成する
  • X^2 + Y^2 < R^2を満たすようなデータだけを残す

こうすることで、半径R内に一様に分布した乱数が得られそうだぞと、で、実際にやってみる。

#一様乱数から生成
df <- data.frame(x=runif(N, min=-R, max=R), y=runif(N, min=-R, max=R)) %>% filter(sqrt(x^2+y^2)<R)
ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")

f:id:teramonagi:20141113093258p:image

散布図を見る限り、確かに一様に分布していそうだぞと、これでいいぞと。

ただ、このやり方はちょっとだけ無駄があって、「X^2 + Y^2 < R^2」を満たすデータのみを取得するというフィルタリングの処理を噛ませているために、発生させた乱数の全てを使っているわけではなく、一部(円からはみ出てしまった部分)を捨ててしまっている。そこでもうちょっと効率的なやり方を考えたい。高校数学でお勉強した"極座標"を思い出すとx, y座標と動径方向r、角度Θは

  •  x = r¥cos(¥theta)
  •  y = r¥sin(¥theta)

という関係で結びつけられていた。従って、円内に一様に分布する乱数を得るために

  • 動径方向rの成分として[0, R]の間の一様乱数を生成する
  • 角度Θ方向の成分として[0, 2π]の間の一様乱数を生成する
  • XとYを上記の数式に従って変換し、生成する

というようにしてやれば半径R内に一様に分布した乱数が得られそうだぞと、で、実際にやってみる。

#動径・角度方向に分けてサンプリング
theta <- runif(N, min=0, max=2*pi)
r <- runif(N, min=0, max=R)
df <- data.frame(x=r*cos(theta), y=r*sin(theta))
ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")

f:id:teramonagi:20141113093257p:image

あ、あれ?この結果の図を見ると明らかに中心部分の方にサンプルが集中していて、外側が手薄になっている、つまり、一様に分布していないぞと、これはよくないぞということが散布図からわかる。

この原因は今度は大学の数学で習う、重積分変数変換公式における面積(要)素の変換、

  •  dxdy = r dr d¥theta

のrを完全に見ていない、つまりそのヤコビアンを無視しているからだ*1。従って、動径方向に対しては半径rに比例するように乱数を引いてやらないといけなくて、[0, R]の一様な乱数uを逆関数数を用いて、

  •  r = ¥sqrt{2u}

と変換して動径方向の乱数rを生成すればよいことがわかる。これは動径方向の累積分関数P(R)が

  •  P(R) = ¥int_0^{R} r dr = ¥frac{1}{2}R^2

とかける事から、このP(R)の逆関数をとることで証明される。

これを実際にやってみると

theta <- runif(N, min=0, max=2*pi)
r <- sqrt(2*runif(N, min=0, max=R))
df <- data.frame(x=r*cos(theta), y=r*sin(theta))
ggplot(df, aes(x=x, y=y)) + geom_point(color="blue")

f:id:teramonagi:20141113093256p:image

散布図より、一様に分布してるのがわかる。

*1:そもそも"一様"というのは考えている空間における抽象的な意味での面積素上に一様に分布するということなんで、こいつを無視してはいかん

2014-11-12

F#でプレゼンテーション/スライドを作りたい

| 00:25 | F#でプレゼンテーション/スライドを作りたいを含むブックマーク

ぼけーっとFSharpのハッシュタグを眺めていたらFsRevealというプロジェクトがあることを知った。

こいつを使うと驚くほど簡単にF#でプレゼンテーション/スライドを作ることができるそうだ。すげぇ!やってみたい!プロジェクトのサイトは以下。

使い方

まず、ここからzipファイルをダウンロード以下からして、そのzipファイルを解凍し、解凍先のフォルダのslides/input.mdをMarkDownぽく適当に編集する。

そしてその後、同フォルダにあるbuild.cmdを叩く*1

build.cmd

すると、必要なパッケージのダウンロードビルドが行われてoutput/input.htmlに結果が出力される。そしてそのhtmlファイルがプレゼンになっているわけだ!なんと簡単!更に、F#に関してはF#のコード上にマウスオーバーするとその型情報がポップアップされるのも素敵過ぎる。

そういや、F# Formattingでレポート書こうと思って諦めた時期があったので、今度はそれをやってみよう。。。

*1Windowsの場合

2014-11-11

geom_pathがまだよくわからない

| 19:27 | geom_pathがまだよくわからないを含むブックマーク

geom_pathというか、ggplot2をよく理解できていないメモ。