遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

Rcpp sugarで統計分布の関数を使う

R Rcpp sugar
  • Rcppのsugarというのは、Rcppパッケージで提供される、cppコードとして使う関数なんだけれども、その名称がうまいことR本体の関数と同じように作ってあって、Rの関数を思い出しながらcppを書けるようにしてくれている一群の「Rユーザを甘やかす関数」のこと
  • その一覧は以下のような感じ
+ *, -, /, pow, <, <=, >, >=, ==, !=, !.
Math functions: abs(), acos(), asin(), atan(), beta(), ceil(), ceiling(), choose(), cos(), cosh(), digamma(), exp(), expm1(), factorial(), floor(), gamma(), lbeta(), lchoose(), lfactorial(), lgamma(), log(), log10(), log1p(), pentagamma(), psigamma(), round(), signif(), sin(), sinh(), sqrt(), tan(), tanh(), tetragamma(), trigamma(), trunc().
Scalar summaries: mean(), min(), max(), sum(), sd(), and (for vectors) var().
Vector summaries: cumsum(), diff(), pmin(), and pmax().
Finding values: match(), self_match(), which_max(), which_min().
Dealing with duplicates: duplicated(), unique().
d/q/p/r for all standard distributions.
  • このリストの末尾にあるように、確率分布関数のd/p/q/rが一通りsugar化提供されている
  • でも、ちょっと使うにはコツがいるようなので、メモ
  • まず、Rで言うところのpbeta()が使いたいとする
  • 調べものはこちらのRcpp本家（らしい）サイトにて、"pbeta"で検索
  • そうするとこちらに飛ぶ。さらに、"Definition at line 61 of file Rmath.h."との記載があり、Rmath.hの61行目を見ればよいことがわかる。Rmath.hにリンクがあるので飛んでいけば、確かにコードが見られて、その61行目が確認できる
  • beta分布関係の行は次の通り
     /* Beta Distribution */
   60     inline double dbeta(double x, double a, double b, int lg)         { return ::Rf_dbeta(x, a, b, lg); }
   61     inline double pbeta(double x, double p, double q, int lt, int lg) { return ::Rf_pbeta(x, p, q, lt, lg); }
   62     inline double qbeta(double a, double p, double q, int lt, int lg) { return ::Rf_qbeta(a, p, q, lt, lg); }
   63     inline double rbeta(double a, double b)                           { return ::Rf_rbeta(a, b); }
   64
  • pbetaに関して言えば、double x, double p, double q, int lt, int lgの５引数であることがわかる
  • また、cppでのpbeta関数というのは、さらに別途定義された::Rf_pbetaのことであることがわかる。このRf_は"Rの関数で言うところの"という書き方ルール（Rcppのパッケージを作った人のルール）で、"::"もcppとRとをつなぐための記法。cxxfunction()を使うとき（inline)では、その記法やらルールやらをいちいち書かなくてもよいように"sugar"化しているからcppでpbetaと打てばよくなっている
  • ただし、引数の具合は違う
  • Rのpbeta()は、そのヘルプに以下のように書かれているように、必須引数が３つ、オプション引数が３つであって、「５個」ではない
pbeta(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
  • cppの４つの関数d,p,q,rと、その引数の組み合わせを見てみる
    • d : x,a,b,lg
    • p : x,p,q,lt,lg
    • q : a,p,q,lt,lg
    • r : a,b
  • rには２引数しかない。Rのrbetaには、発生させる乱数の個数の引数と、ベータ分布を決める２引数を与えるが、cppのrには２つしかない。２つはベータ分布を決めるのに必須だから、cppのrbeta関数は、乱数を１個返す関数で、a,bはベータ分布の形を決める引数であることがわかる
    • 確認してみる。cppのrbeta関数は発生乱数の個数もとるんだな、これが…
// testbeta.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;  // use the Armadillo library for matrix computations
using namespace Rcpp;


// [[Rcpp::export]]
double testRbeta(double a, double b){
  double ret = Rf_rbeta(a,b);
  return ret;
}

// [[Rcpp::export]]
vec testRbeta2(int n, double a, double b){
//これってどうして動くんだろう…
  vec ret = rbeta(n,a,b);
  return ret;
}

// [[Rcpp::export]]
double testRbeta3(double a, double b){
  double ret = ::Rf_rbeta(a,b);
  return ret;
}

Rcpp::sourceCpp("testbeta.cpp")
testRbeta(3,5)
testRbeta2(100,3,5)
testRbeta3(3,5)
  • pbetaでやってみよう
// testbeta.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;  // use the Armadillo library for matrix computations
using namespace Rcpp;

// [[Rcpp::export]]
double testPbeta(double x,double p, double q,int lt=0,int lg=0){
  double ret = Rf_pbeta(x,p,q,lt,lg);
  return ret;
}
  • 以下の様子からltはlower.tailのT/Fに、lgはlog.pのT/Fに対応していることがわかる
> testPbeta(0.6,3,5,1)
[1] 0.903744
> testPbeta(0.6,3,5)
[1] 0.096256
> testPbeta(0.6,3,5,0)
[1] 0.096256
> pbeta(0.6,3,5,lower.tail=FALSE)
[1] 0.096256
> testPbeta(0.6,3,5,1)
[1] 0.903744
> pbeta(0.6,3,5,lower.tail=TRUE)
[1] 0.903744
> pbeta(0.6,3,5)
[1] 0.903744
> pbeta(0.6,3,5,log.p=TRUE)
[1] -0.1012091
> testPbeta(0.6,3,5,1,1)
[1] -0.1012091
  • qchisq()でやってみると
// testbeta.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;  // use the Armadillo library for matrix computations
using namespace Rcpp;

// [[Rcpp::export]]
double testQchisq(double p, double df, int lt, int lg){
  double ret = Rf_qchisq(p,df,lt,lg);
  return ret;
}

> testQchisq(0.99,1,0,0)
[1] 0.0001570879
> testQchisq(0.01,1,0,0)
[1] 6.634897
> qchisq(0.01,1,lower.tail=FALSE)
[1] 6.634897
> qchisq(0.01,1,lower.tail=TRUE)
[1] 0.0001570879
> testQchisq(0.01,1,1,0)
[1] 0.0001570879
    • cppのRf_qcihsq()の最後のint lgは0以外の値はとれないようだ…