Hatena::ブログ(Diary)

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

2015-02-28

RmdファイルからRコードを抜くのはpurl関数

| 10:18 | RmdファイルからRコードを抜くのはpurl関数を含むブックマーク

例えば、以下のような内容のR Markdownファイル(hoge.Rmd)を用意していたとする。

げっへっへ、本文だぜー!!
```{r, cache=TRUE}
#Rのコード
x <- 1:10
```
また、本文だぜー!

ここから、Rのコードだけを抽出するためにはknitrパッケージにあるpurl関数を用いると良い。documentation引数に応じて出力が変わる(チャンク・オプション、本文のコメント化など)ので、以下にその結果をメモる。

documentation=0

コマンド

purl("hoge.Rmd", documentation=0)

結果


#Rのコード
x <- 1:10

documentation=1(デフォルト)

コマンド

purl("hoge.Rmd", documentation=1)

結果


## ----, cache=TRUE--------------------------------------------------------
#Rのコード
x <- 1:10

documentation=2

コマンド

purl("hoge.Rmd", documentation=2)

結果

#' げっへっへ、本文だぜー!!
## ----, cache=TRUE--------------------------------------------------------
#Rのコード
x <- 1:10

#' また、本文だぜー!

2015-02-13

無敵(不死)タイムのある指数分布っぽい確率分布からのサンプリング

| 19:18 | 無敵(不死)タイムのある指数分布っぽい確率分布からのサンプリングを含むブックマーク

こいつ

の応用で、

  • 時点1〜5までは無敵(不死)

になるような確率分布からのサンプリングをやってみた。

具体的には

lambda <- function(t){ifelse(t > 1 & t < 5, 0, 0.75)}

というように、デフォルト強度(ポアソン分布の平均値に相当する量) ¥lambdaを一時的に0と設定する。こうすると単位時間当たりの死亡率が0になるので、この間は無敵となるのだ。んで、この形からサンプリングすると、以下のように絶対死なない”無敵タイム付”の指数分布からのサンプリングになる。ちゃんと1〜5の間は確率密度が0になるので、この間は死なないよと。

f:id:teramonagi:20150213173220p:image

コードは以下、コピペすれば同じ結果を再現できるぞい。

#int_0^t \lambda(u)du - \xi
#な関数。無理やり根を探す形で書くのでこうした
objective <- function(t, lambda, xi){
  integrate(lambda, 0, t)$value - xi
}
#適当な強度の関数形
lambda <- function(t){ifelse(t > 1 & t < 5, 0, 0.75)}
#描画範囲
range_plot <- c(0, 15)
#サンプリング
set.seed(100)
x <- sapply(
  1:10^3, 
  function(i){
    uniroot(objective, range_plot, tol=10^(-3), lambda=lambda, xi=rexp(1))$root
  }
)

library(ggplot2)
qplot(x, geom="blank") + 
  geom_histogram(aes(y=..density..), fill="lightpink2") + 
  scale_x_continuous(limits=range_plot)

2015-02-11

非同次/確率的な強度λを持つ指数分布に従う確率変数τ(死亡/消滅時刻)の生成

| 19:34 | 非同次/確率的な強度λを持つ指数分布に従う確率変数τ(死亡/消滅時刻)の生成を含むブックマーク

非同次な指数分布

非同次なポアソン分布ってのがある。日本語で定義がちゃんと載ってるページがなかったので、英語版wikipediaを参照すると

が該当する。要するにこれは、今まで、一定だと思っていたポアソン分布の強度(普通のポアソン分布だと平均だと思っていたパラメータ) ¥lambdaを時間や空間に依存させてみたりするってもんだ。あと、この拡張として、この ¥lambda自身を確率過程だと思っても以下の話は大体同じになる。正確な測度論の記述じゃないけど。

そして、

の資料に説明があるように、その非同次なポアソン分布の裏側には、非同次な指数分布が生まれてるんだろうなぁと考えることが出来るわけで、ここではそれを考えたい。具体的には、非同次な指数分布に従う乱数を生成するにはどうしたらよいのかという話だ。日本語の教科書にこの辺の記載があるのかないのかもよくわからないが、金融だと、

の「3.5 Process with Jumps」や「9.4 Credit Risk」を見ておくと良い。日本語だとたぶん信頼性工学系の書籍を漁ることになるんだとおもう、多分。

非同次な指数分布の数式展開

さて、(非)同次な指数分布に従う確率変数 ¥tauと書くことにすると、まず、同次な指数分布の分布関数(時刻 Tまでに死亡/消滅していない確率)は

 P(¥tau < T) = 1-¥exp¥left¥{-¥lambda T ¥right¥}

と書くことができるが、今、非同次なケースを考えようとすると、強度 ¥lambda変数(時間・空間)依存性を持つので、以下のように積分形式で書くことができる。

 P(¥tau < T) = 1-¥exp¥left¥{-¥int_0^T ¥lambda(t)dt ¥right¥}

これは ¥lambdaが一定値だと思うと、ちょうど同次のケースに戻るので、これで定義としては良さそうだぞと、逆に、これを持って非同次な指数分布の分布関数とするぞとそういうことです。

ここで、ここから先、逆変換法を使っていくが、逆変換法については

あたりを見てください*1

さて、非同次な指数分布に対して逆変換法を適用するために区間[0, 1]の一様乱数 uを導入し、

 u = 1-¥exp¥left¥{-¥int_0^T ¥lambda(t)dt ¥right¥}

と書く。あとは、これをTについて解いてやればよいので、ちょっと変形して

 ¥log(1-u) = ¥int_0^T ¥lambda(t)dt

とする。ここから先は解析的に解くのがキツいので、こいつはこのままで放置するとして、この数式を満たすような最小の Tが非同次な指数分布に従うというわけだ。もう少しちゃんと書くと

 ¥tau = ¥inf¥left¥{ T>0; ¥int_0^T ¥lambda(t)dt = ¥log(1-u) ¥right¥}

という形になる。更に、ここで、

 ¥log(1-u)

は逆変換法の視点から解釈するに、 ¥lambda=1とした指数分布になっているので、それを ¥xiと書くことにすると、

 ¥tau = ¥inf¥left¥{ T>0; ¥int_0^T ¥lambda(t)dt = ¥xi ¥right¥}

となる。なので、シミュレーションの手順としては

1:  ¥lambda=1(平均1)とした指数分布から一個乱数 ¥xiを取得する

2:  ¥int_0^T ¥lambda(t)dtの計算を数値積分なり解析解で頑張って計算する

3:  ¥int_0^T ¥lambda(t)dt = ¥xiを満たすような Tを求めて、こいつを非同次な指数分布からの乱数だと思う

という手順になる。

Rでやってみる

というわけで、上で書いたアルゴリズムを実際に書いてみた。ここでは、

 ¥lambda(t) = 15t^2 + 1

と設定している。深い意味はなく、解析的に計算したかったのでこうしている。この関数積分は解析的に計算できるので、全体として確率密度関数 f(T)

 f(T)dT = P(T < ¥tau < T+dt) = P(¥tau < T+dt) - P(¥tau < T) ¥sim ¥lambda(T) ¥exp¥left¥{ -¥int_0^T ¥lambda(t)dt ¥right¥} = (15T^2 + 1)¥exp¥left¥{ -(5T^3+T) ¥right¥} dT

となるんで、これとサンプリング結果を比較する。結果をPLOTすると、以下のようにほぼ一致しているので、これで良さそうだ。

f:id:teramonagi:20150211203836p:image

ソースは以下。コメントにあるように平均1の指数分布と強度 ¥lambda積分が一致する点はめんどくさいのでソルバー(uniroot)関数に投げちゃってる。あと、乱数生成は汎用的にしたかったので数値積分で処理しているので遅いぞと。

#int_0^t \lambda(u)du - \xi
#な関数。無理やり根を探す形で書くのでこうした
objective <- function(t, lambda, xi){
  integrate(lambda, 0, t)$value - xi
}
#適当な強度の関数形
lambda <- function(t){5*3*t^2+1}
#描画範囲
range_plot <- c(0, 2)
#サンプリング
set.seed(100)
x <- sapply(
  1:10^4, 
  function(i){
    uniroot(objective, range_plot, tol=10^(-3), lambda=lambda, xi=rexp(1))$root
  }
)
#解析解との比較
#上のラムダを積分(手でがんばる)
lambda_integrated <- function(t){5*t^3+t}
#解析解(確率密度関数)
analytical_solution <- function(t){
  lambda(t)*exp(-lambda_integrated(t))
}
#比較PLOT
library(ggplot2)
qplot(x, geom="blank") + 
  geom_histogram(aes(y=..density..), fill="aquamarine",color="black") + 
  stat_function(fun=analytical_solution, color="red", size=1) + 
  scale_x_continuous(limits=range_plot)

## 追記

「要するにCox回帰モデルで作る指数分布ってこと?」という突っ込みをいただきましたが、その解釈でもよいと思います。なので、Cox回帰モデルに従うサンプルデータを逆に生成できるぞと、データの水増しに使うことができるぞとそういうことです。

*1:他にいい資料がググったら出てくるかもしれないが、知らん

2015-02-08

shinyappsのデプロォーイ(deploy)周りのメモ

| 10:16 | shinyappsのデプロォーイ(deploy)周りのメモを含むブックマーク

はじめに

今、日本の国債金利を自動でデータ更新して可視化するようなもんを作ろうとしてるんですが

こいつを作るときにあわあわした点をメモっておく。

認証名のトークンパスワードは.Rprofileに書いておくのがよさげ

以下を参考に、.Rprofileファイル内に、optionなんかを使って設定しておくのがよいのかなと。これがベストプラクティスなのかは知らない。

optionを自分で追加できるってのがミソ。

> options(teramonagix=1:10)
> options()$teramonagix
 [1]  1  2  3  4  5  6  7  8  9 10

今回の場合だと

shinyapps::setAccountInfo(
  name='teramonagi', 
  token=options()$shinyapps.token, 
  secret=options()$shinyapps.secret
)

というように書いている。

ディレクトリ名を表す際の末尾のパス・セパレータはいれちゃいかん!!!on Windows

デプロイをバッチ叩くだけでいいようにしたくて、

deployApp(paste0(getwd(), "/../"))

ってバッチの中で書いてたんだけど、これはだめだった on Windows。正しくは

deployApp(paste0(getwd(), "/.."))

これは多分

に書いてある

Trailing path separators are invalid for Windows file paths apart from ‘/’ and ‘d:/’ (although some functions/utilities do accept them), so as from R 3.1.0 a trailing / or \ is removed.

ってことっぽくて、例えば以下のコードで似たような事例を確認することができる。

> path <- "C:/temp/moge"
> file.exists(path)
[1] TRUE
> file.exists(paste0(path, "/../"))
[1] FALSE
> file.exists(paste0(path, "/.."))
[1] TRUE

今回のケースでは、上記のdeployApp関数の中で呼ばれているfile.exists関数がFALSEを返してきやがるので、処理が途中で終了してバッチが上手く走らなかったというお話。

httpsアクセス内でのread.csv

引数のファイル(URL)は"http://~~~"じゃないとだめだった。はじめhttpsにしていたのだが、アプリが起動しなかった。この辺のエラーログはshinyappsの管理画面のLOGから取れるので、ありがたい&ググりやすい。

Rの環境設定周りのお話

R_HOMEだの.RProfileだのはめんどいので使っていなかったのだが、これがあると便利なので、今回は設定した。以下を見ておけば大体OK。この資料内で言及されているPDFも良い。

ID名など

output系のIDに、ハイフン"-"を使ってはいけない(アンダースコア"_"はOK)。これはRの変数名のルールとしてハイフンが使えないからで、具体的にはこのIDに対して、処理を割り振る時に、

output$hoge-hoge

と各必要があり、ここで処理が落ちるから。

hoxo_mhoxo_m 2015/02/08 10:38 RStudio の上の方(Debug, Tools の一段下)に Publish ボタンがありませんか?
それを押すと幸せになれるかもしれません。

hoxo_mhoxo_m 2015/02/08 11:09 ご参考まで。
http://qiita.com/hoxo_m/items/141282a5d1ca0fc4f827

teramonagiteramonagi 2015/02/08 17:27 ありがとうございます!灯台下暗しってやつですな。。。

2015-02-07

F#の配列(Array)を辞書(Dictionary)のキーにしちゃいかん!!代わりにタプル(tuple)を使おう

| 19:27 | F#の配列(Array)を辞書(Dictionary)のキーにしちゃいかん!!代わりにタプル(tuple)を使おうを含むブックマーク

すごいハマった。

let dic = new System.Collections.Generic.Dictionary<int option [], string>()
dic.[[|Some 1; Some 1|]] <- "A"
dic.[[|Some 1; Some 2|]] <- "B"
dic.[[|Some 1; Some 3|]] <- "C"

というコードをを実行しても

System.Collections.Generic.KeyNotFoundException: 指定されたキーはディレクトリ内に存在しませんでした。
   場所 System.Collections.Generic.Dictionary`2.get_Item(TKey key)
   場所 <StartupCode$FSI_0141>.$FSI_0141.main@()
エラーのため停止しました

となって死ぬ。多分、キーとして値丸コピーしてるんじゃなくて、アドレスだけメモってるんだろうな。代わりにtuple使っとけってことっぽい。こっちはちゃんと取れる。

let dic = new System.Collections.Generic.Dictionary<(int option*int option), string>()
dic.[(Some 1, Some 1)] <- "A"
dic.[(Some 1, Some 2)] <- "B"
dic.[(Some 1, Some 3)] <- "C"
dic.[(Some 1, Some 3)]

Enumの全要素をゲッツしたい

| 13:05 | Enumの全要素をゲッツしたいを含むブックマーク

こういうEnumがあったとする。

type Direction = 
    | Up = 0
    | Right = 1 
    | Down = 2
    | Left = 3

その要素を全ゲッツしたければ

System.Enum.GetValues(typeof<Direction>)

とすればよい。System.Enumを叩けば、Enum周りのものはだいたい出てくると覚えよう。これも.NETの機能か。

System.Arrayからのキャスト

上の結果はSystem.Arrayという型になっているが、これだと使いにくいので、以下のようにcastするのがよさげ。

System.Enum.GetValues(typeof<Direction>) |> Seq.cast |> Seq.map unbox<Direction>;;

参考