今まで R パッケージユーザでしたが, 2018年は何かしら R パッケージを開発できたらと考えています。開発を行う上でRにおけるメタプログラミングの知識は必要と思い調べてみました。
- 表現式
- 表現式の捕捉
- 関数呼び出しの構築
- 表現式の評価
- stats::lm() の仕組み
- モデルの更新
表現式
Rプログラミング本格入門 から R におけるメタプログラミングについて言及している箇所を引用する。
メタプログラミングによって, 言語そのものに手を加えて, 特定のシナリオで使われる特定の構文を簡単に作ることができる.
メタプログラミングによって柔軟性の高いプログラミングを可能にする。R でメタプログラミングが活用されるのは主に非標準評価 （Non-Standard Evaluation; NSE） である。
まず, R でメタプログラミングを行う上で基本となる 表現式 (expression) について調べてみた。
表現式は R で実行できる処理を表すオブジェクトで, そのデータ構造は抽象構文木 (abstract syntax tree; AST) である。表現式を構成する要素は, 定数・名前・呼び出し・ペアリストの4つである。表現式は評価が成功すると値となる。
表現式が AST であることを pryr パッケージを使い確認してみる。 pryr パッケージは R の言語としての側面を深く理解するために便利なツール。
> pryr::ast(sqrt(1 + x^2))
\- ()
\- `sqrt
\- ()
\- `+
\- 1
\- ()
\- `^
\- `x
\- 2
表現式の捕捉
表現式が評価される前に, 表現式自体を変数に格納することを表現式の捕捉 (capture) という。表現式を捕捉する関数に quote() や substitute() がある。
R ではあらゆる操作が関数呼び出しとなるので, quote() や substitute() で全てを捕捉できることができる。
base::quote
quote() は関数呼び出しが捕捉された場合は呼び出しオブジェクトを返す。つまり評価される前の処理内容を捕捉できる。
> call1 <- quote(sqrt(1 + x^2))
call1 のオブジェクトを調べると, オブジェクトの型が language でクラスが call と出力されることから, 捕捉された表現式自体を格納した変数は言語オブジェクトで呼び出しオブジェクトであることがわかる。
> typeof(call1)
[1] "language"
> class(call1)
[1] "call"
> is.call(call1)
[1] TRUE
> pryr::call_tree(call1)
\- ()
\- `sqrt
\- ()
\- `+
\- 1
\- ()
\- `^
\- `x
\- 2
呼び出しオブジェクトの内部構造は list に変換することで確認できる。
> as.list(call1)
[[1]]
sqrt
[[2]]
1 + x^2
言語オブジェクトには呼び出しオブジェクト以外にも, 表現式オブジェクト, シンボルオブジェクトがある。
base::substitute
substitute() は入力された表現式を捕捉し表現式中のシンボルを環境中の値と置き換えた表現式を返す。
> func1 <- function(x) substitute(x+1)
> func1(sqrt(2))
sqrt(2) + 1
> eval(func1(sqrt(2)))
[1] 2.414214
pryr::modify_call でも呼び出しを変更できる。
> x <- c(1, 2 ,3 ,NA ,NA)
> call2 <- quote(mean(x, na.rm = TRUE))
> call2
mean(x, na.rm = TRUE)
> eval(call2)
[1] 2
> pryr::modify_call(call2, list(na.rm = FALSE))
mean(x = x, na.rm = FALSE)
> eval(pryr::modify_call(call2, list(na.rm = FALSE)))
[1] NA
substitute() は単一の呼び出しを変更できるが, より複雑なタスクに対しては AST を扱える codetools パッケージが向いている。
関数呼び出しの構築
base::call
call() は関数名を表す文字列 (シンボル) とその引数を指定することで呼び出しオブジェクトを生成する。
> call3 <- call("sqrt", 2)
> call3
sqrt(2)
> typeof(call3)
[1] "language"
> class(call3)
[1] "call"
call() はそのまま位置で引数を組み立てるため, 名前付き引数を指定できる match.call() の方が安全。pryr::standardise_call では内部で match.call() を使い引数をすべて名前付き引数として並べ直している。
> call("read.csv", header = FALSE, "", "sample.csv")
read.csv(header = FALSE, "", "sample.csv")
> match.call(read.csv, call("read.csv", header = FALSE, "", file = "sample.csv"))
read.csv(file = "sample.csv", header = FALSE, sep = "")
> pryr::standardise_call(call("read.csv", header = FALSE, "", file = "sample.csv"))
read.csv(file = "sample.csv", header = FALSE, sep = "")
表現式の評価
base::eval
eval() は表現式を評価する。eval() が quote() の皮を一枚剥ぐことを確認してみる。
> 1:10
[1] 1 2 3 4 5 6 7 8 9 10
> quote(1:10)
1:10
> quote(quote(1:10))
quote(1:10)
> eval(quote(quote(1:10)))
1:10
> eval(eval(quote(quote(1:10))))
[1] 1 2 3 4 5 6 7 8 9 10
stats::lm() の仕組み
stats::lm のように formula を指定してモデリングを行えるのが R のひとつの特徴である。
model <- lm(mpg ~ wt, data = mtcars)
lm.R の lm() の流れを少し追ってみる。まず, 以下で model の初期化を行う。
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action", "offset"),
names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
## need stats:: for non-standard evaluation
mf[[1L]] <- quote(stats::model.frame)
mf <- eval(mf, parent.frame())
続いて method のチェックなどを行った後, lm.fit() または lm.fwit() で最小二乗法によりパラメータを推定する。
x <- model.matrix(mt, mf, contrasts)
z <- if(is.null(w)) lm.fit(x, y, offset = offset,
singular.ok=singular.ok, ...)
else lm.wfit(x, y, w, offset = offset, singular.ok=singular.ok, ...)
lm.fit() は z <- .Call(C_Cdqrls, x * wts, y * wts, tol, FALSE) でC言語の関数を呼ぶ。呼ばれる関数は SEXP Cdqrls(SEXP x, SEXP y, SEXP tol, SEXP chk) でQR分解で最小二乗解を計算する Fortran の dqrls() の Wrapper 関数である。(lm.c, dqrls.f)
SEXP はC言語で R オブジェクト (S表現式) を表すデータ型で, C言語に渡す引数や返り値はこの SEXP 型である必要がある。
F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
REAL(coefficients), REAL(residuals), REAL(effects),
&rank, INTEGER(pivot), REAL(qraux), work);
さらに dqrls() の中では LINPACK の機能を使っている。
最終的に lm() は得られた z オブジェクトに属性を追加し返す。
class(z) <- c(if(is.matrix(y)) "mlm", "lm")
z$na.action <- attr(mf, "na.action")
z$offset <- offset
z$contrasts <- attr(x, "contrasts")
z$xlevels <- .getXlevels(mt, mf)
z$call <- cl
z$terms <- mt
if (model)
z$model <- mf
if (ret.x)
z$x <- x
if (ret.y)
z$y <- y
if (!qr) z$qr <- NULL
z
モデルの更新
lm() で得られたモデルは stats::update で更新することができる。update.R の update() から更新の仕組みを調べてみる。
getCall <- function(x, ...) UseMethod("getCall")
getCall.default <- function(x, ...) getElement(x, "call")
## Using getCall() instead of x$call renders update.default() more
## generally applicable.
update.default <-
function (object, formula., ..., evaluate = TRUE)
{
if (is.null(call <- getCall(object)))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$...
if (!missing(formula.))
call$formula <- update.formula(formula(object), formula.)
if(length(extras)) {
existing <- !is.na(match(names(extras), names(call)))
## do these individually to allow NULL to remove entries.
for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
if(any(!existing)) {
call <- c(as.list(call), extras[!existing])
call <- as.call(call)
}
}
if(evaluate) eval(call, parent.frame())
else call
}
update.formula <- function (old, new, ...)
{
tmp <- .Call(C_updateform, as.formula(old), as.formula(new))
out <- formula(terms.formula(tmp, simplify = TRUE))
return(out)
}
上記は, まず getCall() で object の呼び出しを抽出する。 getCall(x) は x$call と同様。また match.call() でモデルを生成するのに使われた … 部分を捕捉する。
次に, 新たな formula が空でないことを確認した上で update.formula() で formula を更新する。最後に更新された呼び出しを eval() により現在の呼び出し環境で評価し結果を返す。(evaluate = FALSE の場合は更新された表現式を返す)
update() を使ってみる。以下の例では getCall() で lm(formula = mpg ~ wt, data = mtcars) が抽出され, この呼び出しが lm(formula = mpg ~ wt + cyl, data = mtcars) に更新・評価されていることがわかる。
> model <- lm(mpg ~ wt, data = mtcars)
>
> model
Call:
lm(formula = mpg ~ wt, data = mtcars)
Coefficients:
(Intercept) wt
37.285 -5.344
>
> model2 <- update(model, formula = . ~ . + cyl)
> model2
Call:
lm(formula = mpg ~ wt + cyl, data = mtcars)
Coefficients:
(Intercept) wt cyl
39.686 -3.191 -1.508
R言語徹底解説 でも紹介されているが, pryr を用いて update() に似た関数 update_model() を作ってみる。
> update_call <- function(object, formula., ...) {
+ call <- object$call
+
+ if (!missing(formula.)) {
+ call$formula <- update.formula(formula(object), formula.)
+ }
+
+ pryr::modify_call(call, pryr::dots(...))
+ }
>
> update_model <- function(object, formula., ...) {
+ call <- update_call(object, formula., ...)
+ eval(call, parent.frame())
+ }
>
> model3 <- update_model(model, formula = . ~ . + cyl)
> model3
Call:
lm(formula = mpg ~ wt + cyl, data = mtcars)
Coefficients:
(Intercept) wt cyl
39.686 -3.191 -1.508
identical() で update() と update_model() の結果を比較すると等価であることが確認できる。
> identical(model2, model3)
[1] TRUE
環境を捕捉すると環境内のオブジェクトも捕捉するため確保されるメモリに注意する必要がある。
おわりに
参考文献は『R言語徹底解説』と『Rプログラミング本格入門』の2冊です。
[1] こわくないRパッケージ開発！2016
[2] Raw quotation of an expression
[3] Rのメタプログラミングとモナド