Hatena::ブログ(Diary)

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

2015-05-19

引数があったりなかったりする処理

| 19:02 | 引数があったりなかったりする処理を含むブックマーク

今まで使ったことなかった。引数が存在しているか否かはmissing関数で判定するということだ。関数デフォルト引数との使い分けをどうしたらいいのかはわからない。

f <- function(x, y)
{
  if(missing(y)){return(x)}
  if(missing(x)){return(2*y)}
  x + y
}

これで、二変数関数でも一変数関数のように処理することができる。

> f(21)
[1] 21
> f(y=21)
[1] 42
> f(21, 32)
[1] 53

・・・っていうのを

みて気が付いた。

2015-05-18

C言語の関数ポインタを取るインターフェイスにRcpp::Functionを渡したい

| 17:34 | C言語の関数ポインタを取るインターフェイスにRcpp::Functionを渡したいを含むブックマーク

最近だとtemplateの魔力を使ったC++を使えばこんなことする必要はないですが、古いC言語で書かれたライブラリなんかに、関数ポインタ引数として渡す必要がある場合のお話。

以下のコードでは、関数ポインタ(typedef double(*Simulate)(double x))を引数にとるold_c_function関数にRcpp::Functionを投げ込んでみた例。old_c_function関数自体は適当な処理(x0 + f(x0))を実行する関数としています。

library(Rcpp)
sourceCpp(code='
  #include <Rcpp.h>
  typedef double(*Simulate)(double x);
  //Old-C function
  double old_c_function(double x0, Simulate f)
  {
    return(x0 + f(x0));
  }
  // [[Rcpp::export]]
  double new_interface(double x, Rcpp::Function f)
  {
    return old_c_function(x, &f);
  }
')

これをコンパイルしようとすると、これは無常にもというか、当たり前ではありますが「型変換できねぇよ!」と怒られる。

g++ -m64 -I"C:/R/R-32~1.0/include" -DNDEBUG     -I"C:/R/R-3.2.0/library/Rcpp/include" -I"C:/Users/teramonagi/AppData/Local/Temp/RtmpMR3pkQ"  -I"d:/RCompile/r-compiling/local/local320/include"     -O2 -Wall  -mtune=core2 -c file16514334dad9.cpp -o file16514334dad9.o
file16514334dad9.cpp: In function 'double new_interface(double, Rcpp::Function)':
file16514334dad9.cpp:12:32: error: cannot convert 'Rcpp::Function* {aka Rcpp::Function_Impl<Rcpp::PreserveStorage>*}' to 'Simulate {aka double (*)(double)}' for argument '2' to 'double old_c_function(double, Simulate)'
file16514334dad9.cpp:13:3: warning: control reaches end of non-void function [-Wreturn-type]

これをなんとかしたい、速度とか気にしてRcpp使ってるんだけど、遅くなってもいいから無理やりRcpp::Function使いたいって時は、1枚コンバータ(アダプタ?)を噛ませばいい。staticとしてactionメソッドを取っているのがミソで、普通のメンバー関数だとポインタが変換できなくてダメ。

sourceCpp(code='
  #include <Rcpp.h>
  typedef double(*Simulate)(double x);
  //Old-C function
  double old_c_function(double x0, Simulate f)
  {
    return(x0 + f(x0));
  }
  //Converter for Rcpp::Function <--> Simulate
  struct Converter
  {
    static double action(double x){return Rcpp::as<double>((*f_)(x));}
    static Rcpp::Function* f_;
  };
  Rcpp::Function* Converter::f_ = NULL;
  // [[Rcpp::export]]
  double new_interface(double x, Rcpp::Function f)
  {
    Converter::f_ = &f;
    return old_c_function(x, Converter::action);
  }
')

これはコンパイルできて、ちゃんと動く。以下は「3 + 3^2 + 2 = 14」という計算に対応し、答えもあってる。

> new_interface(3, function(x){x^2+2})
[1] 14

2015-05-14

線形代数ライブラリEigenの資料まとめ

| 04:58 | 線形代数ライブラリEigenの資料まとめを含むブックマーク

C++のテンプレートメタプログラミングを活用した線形代数ライブラリであるEigenに関連した資料へのLINKまとめ。

公式ドキュメント

とりあえずここを探す。

特にクイックリファレンスガイド良く見る。

講演資料

特に

Eigen a c++ linear algebra library

が、「Expression Template使ってるぜー」とかの説明もあって良かった。

eigen/Cookbook

短いけどCookbook。うっかりshared_ptrに突っ込んでたので、しょっぱなの「Structures containing Eigen types as members」から勉強になったわ…

Eigen ー C++線形代数を!

flyio氏のまとめ。とりあえず読んでおけ。

行列ライブラリEigenのメモ

暗黒通信団?のメモ(PDF)

Eigen - C++で使える線形代数ライブラリ

C++ Advent Calendar 2014 - Qiitaの記事

調べてた時のツイートまとめ

はじめはArmadillo使おうかなと思っていたが、

ということでEigenにした。

また、より高速なBlazeなんてのもあるらしい。こちらはBoost + Blas必須。

↑上記のLINKのコメ欄でEigen勢から結構突っ込まれてたから、Armadilloの方がいいとは言えないってのにここで気がついた。

2015-05-12

カイ二乗値を計算する

| 22:03 | カイ二乗値を計算するを含むブックマーク

今まであまり考えないで検定用の関数に突っ込んでたんで実際に計算してみんとす。

例が悪くてあれなんだが、以下のような取引戦略A・B、およびその勝敗のデータがあったとする。

> x <- matrix(c(110, 119, 200, 207), 2, 2)
> rownames(x) <- paste0("Strategy", LETTERS[1:2])
> colnames(x) <- c("WIN", "LOSE")
> x
          WIN LOSE
StrategyA 110   200
StrategyB 119   207

これは

  • 取引戦略A:試した310日中、110日勝った
  • 取引戦略B:試した326日中、119日勝った

というデータだと読めばよい。

(※コメント欄の指摘が正しいので、修正したらどっちもクソ戦略になっちゃった・・・)

で、「取引戦略Bの方が勝率が高いが、それって統計的に有意と言えるんでしたっけ?」という場合にはカイ二乗検定をして見ればよい。

> chisq.test(x, correct=FALSE)

	Pearson's Chi-squared test

data:  x
X-squared = 0.0716, df = 1, p-value = 0.789

結果、p値が0.789と大きく、「勝率が等しい」とする帰無仮説棄却することができず有意じゃないですねと。

で、ここで出てくる「X-squared = 0.0716」って量は、元のデータから

> expected <- rowSums(x)%o%colSums(x)/sum(x)
> sum((x-expected)^2/expected)
[1] 0.07163452

として計算できますね・・・というのを確認したかった。

Wikipediaに載ってる数式、雑すぎるんだよね。

akagi_paonakagi_paon 2015/05/12 23:38 WIN TRIAL ではなく WIN LOSE では?

teramonagiteramonagi 2015/05/14 08:22 うおお、本当だ。直しておきますorz ありがとうございました!

2015-05-10

可視化で理解する中心極限定理(ベルヌーイ分布編)

| 21:53 | 可視化で理解する中心極限定理(ベルヌーイ分布編)を含むブックマーク

こういう話がある。

非常に素晴らしいので、指数分布じゃなくて、ベルヌーイ分布版をアニメーションにしてみた。

下図は

となっている。

データが溜まってくると平均値の推定値の精度があがる(赤線の間隔 or 正規分布標準偏差が縮まる)ことが可視化されてわかりやすい。

f:id:teramonagi:20150510214829g:image

コードは以下で、ほぼ@hoxo_m氏の上記の記事のパクリ。ggplot2の勉強になりました。

library(animation)
library(ggplot2)
saveGIF( {
  prob <- 0.2
  mean_p <- prob
  sd_p <- sqrt(prob*(1-prob))
  x <- c()
  for(i in 1:50) {
    x <- c(x, sample(c(1,0), 10, prob=c(prob, 1-prob), replace=TRUE))
    mean_x <- mean(x)
    sd_x <- sd(x)/sqrt(length(x))
    p <- ggplot(data.frame(x=x), aes(x=x)) + 
      geom_histogram(aes(y=..density..), binwidth=0.1, alpha=0.5) + 
      stat_function(fun=dnorm, geom="density", color="blue", fill="blue", alpha=0.5, arg=list(mean=mean_p, sd=sd_p/sqrt(length(x)))) + 
      geom_vline(xintercept=mean_x,      color="red", size=2) +
      geom_vline(xintercept=mean_x+sd_x, color="red", size=1) +
      geom_vline(xintercept=mean_x-sd_x, color="red", size=1) +
      scale_y_continuous(labels = function(x) sprintf("%0.2f", x)) + 
      xlim(0, 1.1) + xlab("ベルヌーイ分布の平均値の分布")      
    print(p)
  }
}, cmd.fun = system, interval = 0.4, ani.width=600, ani.height=400)

推定した平均値のバラつき(正規分布の1σ or データから計算した1σ信頼区間)も、理論・データからの推計値がほぼ一致。

> sd_x
[1] 0.01823362
> sd_p/sqrt(length(x))
[1] 0.01788854