このブログの更新は Twitterアカウント @m_hiyama で通知されます。
Follow @m_hiyama

メールでのご連絡は hiyama{at}chimaira{dot}org まで。

はじめてのメールはスパムと判定されることがあります。最初は、信頼されているドメインから差し障りのない文面を送っていただけると、スパムと判定されにくいと思います。

参照用 記事

R言語メタプログラミングのサンプル

昨日の記事「R言語メタプログラミングの基礎」では、タイトルどおりに基礎的事項を説明しましたが、実例がありませんでした。メタプログラミング機能を使ったサンプルを紹介します。

内容:

  1. 2変数関数の3Dグラフを描きたい
  2. surf()関数の使用例
  3. surf()関数の実装

2変数関数の3Dグラフを描きたい

Rの組み込み関数にcurve()があります。curve(x^2 + x, xlim=c(-2, 2)) のように呼び出すだけで、1変数関数 f(x) = x^2 + x のグラフを描いてくれます。

2変数関数のグラフを描くには、persp()関数が使えますが、curve()のようにお手軽ではありません。事前に関数の値を計算したマトリックスを作っておく必要があります。例えば、「ベルトラン/ボレルのパラドックスから見える確率の本音と建前」でpersp()を使ってグラフを描いていますが、次のような準備をしています。

d <- function(theta1, theta2) {
    abs(theta2 - theta1)
}
N <- 100
x <- seq(0, 2*pi, by=(2*pi / N))
y <- x
z <- outer(x, y, d)

persp()の呼び出しもけっこう複雑です。

persp(x, y, z,
      xlab = "θ1", ylab = "θ2", zlab = "|θ2 - θ1|",
      xlim = c(0, 2*pi), ylim = c(0, 2*pi), zlim = c(0, 2*pi),
      theta = 10, phi = 20, expand = 0.7, col = "skyblue")

描画を制御するパラメータの指定は致し方ないとしても、事前のマトリックス作成は省きたいところです。curve()と同じように、第1引数に関数の式を直接指定するだけでOKなら便利です。curve()の2変数関数版に相当する関数をsurf()(surfaceから)として、surf(x^2 + y^2) のように書くだけで2変数関数の3Dグラフを描けるようにしましょう。

surf()関数の使用例

先に、surf()を使った描画の例をいくつか示しておきます。surf(x^2 + y^2) は次のような図を描きます。

2変数関数の独立変数の名前はデフォルトでx, yですが、変更することができます。例えば、x^2 + y^2 の代わりに、s^2 + t^2 としても同じグラフを描きますが、軸のラベルは自動的に変更されます。surf(s^2 + t^2, xname=s, yname=t) とすると次の図になります、軸ラベルに注目。

2次元の正規分布(密度関数)のグラフは、例えば surf(exp(-4*(x^2 + y^2)) で描けます。

見栄えを良くするためにパラメータを調整した例は次です。


> surf(exp(-4*(x^2 + y^2)),
+ xlim=c(-2, 2), ylim=c(-2, 2), xdiv=100, ydiv=100,
+ zlab="Normal Distribution", theta=20, phi=20, col="pink")

ベルトラン/ボレルのパラドックスから見える確率の本音と建前」と同様な図は、次のようにして一発で描けます。事前の準備は不要です。


> surf(abs(theta2 - theta1), xname=theta1, yname=theta2,
+ xlim = c(0, 2*pi), ylim = c(0, 2*pi), zlim = c(0, 2*pi),
+ xdiv = 100, ydiv = 100, use.outer = TRUE,
+ theta = 10, phi = 20, expand = 0.7)

surf()関数の実装

R言語メタプログラミングの基礎」で定義した2つのユーティリティ関数を使います。念の為に再掲すると:

# デフォルト値なしを表す空な名前を返す
noDefault <- function() {
    f <- function(x){}
    formals(f)[[1]]
}

# 指定された仮引数リスト、本体、環境から関数を作成する
makeFunction <- function(fpars, fbody, fenvir = .GlobalEnv) {
    f <- function(){}
    formals(f) <- fpars
    body(f) <- fbody
    environment(f) <- fenvir
    f
}

メイン関数surf()と、関数の値を格納したマトリックスを計算する補助関数makeZMatrix()は次のようです。要点の説明はソースコードの後にあります。

[追記]ソースコードはgistにも貼り付けてあります。

以下のコードは変更する気はありませんが、gistのほうは変更するかも知れません。[/追記]

# 2変数関数の値(zの値)を成分とする行列を作成する
makeZMatrix <-
  function (
    # 2変数関数の式、コールオブジェクトを渡す
    expr,
    # 2変数関数の第1変数、第2変数の名前、文字列で指定する
    xname = "x", yname = "y",
    # 2変数関数の第1変数、第2変数の描画域
    xlim = c(-1, 1), ylim = c(-1, 1),
    # 描画域を幾つの少区間に分割するか(分割数)
    xdiv = 20, ydiv = 20,
    # 関数値マトリックス作成時にouterを用いるかどうか
    use.outer = FALSE)
{
  x0 <- xlim[1]
  x1 <- xlim[2]
  y0 <- ylim[1]
  y1 <- ylim[2]
  xStep <- (x1 - x0)/xdiv
  yStep <- (y1 - y0)/ydiv

  xBreaks <- seq(from = x0, to = x1, by = xStep)
  yBreaks <- seq(from = y0, to = y1, by = yStep)

  # 関数の仮引数リストの作成
  fplist <- list(noDefault(), noDefault()) # デフォルト値なし2変数
  names(fplist) <- c(xname, yname)  # 引数名を指定
  # 関数の作成
  fun <- makeFunction(fplist, expr) 

  if (use.outer) {
    zMatrix <- outer(xBreaks, yBreaks, fun)
  } else {
    zMatrix <- matrix(NA, xdiv + 1, ydiv + 1)
    for (i in 1:(xdiv + 1)) {
      for (j in 1:(ydiv + 1)) {
        zMatrix[i, j] <- fun(xBreaks[i], yBreaks[j])
      }
    }
           
  }
  list(
    x = xBreaks,
    y = yBreaks,
    z = zMatrix
    )
}

# 2変数関数の曲面(3Dグラフ)を描く
surf <-
  function (
    # 2変数関数の式、引数として式を書く
    z,
    # 2変数関数の第1変数、第2変数の名前、文字列でもシンボルでもどっちでもよい
    xname = "x", yname = "y",
    # 2変数関数の第1変数、第2変数の描画域
    xlim = c(-1, 1), ylim = c(-1, 1),
    # 描画域を幾つの少区間に分割するか(分割数)
    xdiv = 20, ydiv = 20,
    # 関数値マトリックス作成時にouterを用いるかどうか
    use.outer = FALSE,
    # ラベル文字列
    xlab = NULL, ylab = NULL, zlab = NULL,
    # 関数値の描画域
    zlim = NULL,
    # perspのオプション
    theta = 30, phi = 30, expand = 0.5, col = "skyblue",
    # perspへの引数リストを戻り値として出力するか
    output.args = FALSE
            )
{
  z <- substitute(z) # 式を評価せずに取り出す
  xname <- substitute(xname) # 文字列、または評価しない名前を取り出す
  yname <- substitute(yname) # 文字列、または評価しない名前を取り出す
  # 名前を文字列にする
  stopifnot(is.name(xname) || is.character(xname))
  stopifnot(is.name(yname) || is.character(yname))
  xname <- as.character(xname)
  yname <- as.character(yname)

  # 関数値のマトリックスを作成する
  xyz  <- makeZMatrix(z,
                      xname = xname, yname = yname,
                      xlim = xlim, ylim = ylim,
                      xdiv = xdiv, ydiv = ydiv,
                      use.outer = use.outer)

  # perspの引数リストを作成する
  args <- list(
    x = xyz$x,
    y = xyz$y,
    z = xyz$z,
    xlab = if (is.null(xlab)) xname else xlab,
    ylab = if (is.null(ylab)) yname else ylab,
    zlab = if (is.null(zlab)) deparse(z) else zlab,
    xlim = xlim,
    ylim = ylim,
    theta = theta,
    phi = phi,
    expand = expand,
    col = col
    )
  if (!is.null(zlim)) {
    args[["zlim"]] <- zlim
  }

  # 描画関数の呼び出し
  do.call("persp", args)
  # output.argsがTRUEなら引数リストを返す
  if (output.args) args
}

メタプログラミング機能は次のように使っています。

  1. surf()の第1引数zに指定された式をそのまま取り出すためにsubstitute()を使っています。
  2. surf()のxname, yname引数に文字列だけでなく名前(シンボル)も許すために、同じくsubstitute()を使っています。
  3. surf()からpersp()を呼び出すためにdo.call()を使っています。
  4. makeZMatrix()内で、関数の式(コールオブジェクト)と引数名から関数を作るためにユーティリティ関数makeFunction()を使っています。