R による記号微分入門

calculus パッケージ入門

R
YUIMA
Author

司馬博文

Published

6/18/2024

概要

calculus は c++ を通じて数値微分・数値積分を高速に実行するパッケージである．同時に，ほとんどの演算を，純粋に記号操作により実行する機能も持つ．一般の多変数関数を，記号のまま微分，Taylor 展開することができる．

calculus リンク集

概要

R パッケージcalculusは，記号微分と数値微分・積分をシームレスに繋いだ機能を提供するパッケージである．

数値微積分はRcppパッケージとcubatureパッケージを通じてC++バックエンドを用いて提供する．記号微分は外部の記号計算ソフトウェアに依らない実装を提供している．

install.packages("calculus")
  • Einstein の縮約記法
  • Levi-Civita 記号の計算
  • 一般化 Kronecker デルタ
  • Taylor 展開
  • 高階導関数
  • 多変数 Hermite 多項式
  • 常微分方程式
  • 微分作用素
  • numDerivは数値微分による Jacobian, Hessian 計算を提供するが，それ以上の高階の微分ができず，テンソル値関数の微分もできない．
  • tensorAは Einstein の縮約記法を提供するが，２階のテンソルまでに限る．
  • mpolyは１変数の Hermite 多項式を提供するが，多変数には対応していない．
  • pracmaは１変数の Taylor 展開を提供するが，多変数には対応していない．
  • cubatureは多変数関数の数値積分に対応しているが，その他の（直交）座標には対応していない．

特に R は記号計算が苦手である．外部の記号計算ソフトウェアに繋ぐことは試みがあるかもしれないが，R 独自のパッケージ内で扱うシステムは従来なかった．

  • RyacasはYacasという外部の記号計算ソフトウェアへのインターフェースを提供する．
  • caracasは R-Python インターフェースreticulateを通じてSymPyという Python の記号計算ライブラリへのインターフェースを提供する．

これにより，R の，新たな統計推測手法を実装し論文として公開するために用いる言語という性質を活かすことが，パッケージcalculusの目的にある (Guidotti, 2022, p. 4)．

確率過程に対する漸近展開公式では，1000 を超える連立常微分方程式系を解き，その解を通じて多変数 Hermite 多項式の和を計算する必要があり，その際に YUIMA パッケージにおいて記号微分を実行する際に用いられている：

Article Image

YUIMA 入門

確率微分方程式のシミュレーションと推測のためのパッケージ`yuima`の構造と使い方をまとめます．

1 基本

1.1 概観

ベクトル，行列，テンソルはいずれも配列 array として実装されている．特に，ベクトルと行列はテンソルの特殊な場合と理解できるように，実装上も，統一的に Einstein の縮約記法が適用可能である．

すべての関数（＝数学的演算）は，数値バージョンと記号演算バージョンのディスパッチとして実装される．

1.2 演算

numeric, complex, character, expression のいずれかの型を持つ，同次元の配列に対して定義されている．

例えば１次元のcharacterに対しては：

library(calculus)
("a + b" %prod% 1i) %sum% (0 %prod% "c") %diff% (expression(d + e) %div% 3)
[1] "((a + b) * (0+1i)) - ((d + e) / 3)"

ここで，1i, 0, expression(d+e),3はいずれもcharacterに変換されてから実行されている：

typeof(.Last.value)
[1] "list"

この段階では+0などを省くのみで，評価や簡約化はされない．

\[ (a+b)(c+d) \] を \[ a+bc+d \] と解釈してしまうことなどを防ぐため，すべての変数は()で囲まれる．

この挙動はoptions(calculus.auto.wrap=FALSE)で無効にできる．

1.3 評価

関数evaluateが提供されている．

x <- array(letters[1:6], dim=c(2,3))
evaluate(x, var=c(a=1, b=2, c=3, d=4, e=5, f=6))
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

関数evaluateはベクトル化されている：

var <- data.frame(a=c(1,3), b=2:3, c=3:4, d=4:5, e=5:6, f=6:7)
evaluate(x, var=var)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    3    3    4    5    6    7

eval関数では，１つの変数の代入しか行えない．

複数与えられた場合でも，最後の１つのみ評価した結果が返される：

var_list <- list(a=1, b=2, c=3, d=4, e=5, f=6)
eval(parse(text=x), envir=var_list)
[1] 6

evaluate関数と同様の結果を得るには，成分ごとに繰り返し適用する必要がある．例えば次のように：

apply(x, c(1, 2), function(expr) eval(parse(text=expr), envir=var_list))
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

2 微積分

2.1 記号微分

関数derivativeは関数を表す文字列fと，微分する変数名を表すvarの２つの引数を取る．

derivative(f="sin(x)", var="x")
[1] "cos(x)"

多変数も同様．階数は引数orderで指定する： \[ \frac{\partial }{\partial x}\frac{\partial ^2}{\partial y^2}y^2\sin(x)=2\cos(x) \]

derivative(f = "y^2 * sin(x)", var = c("x", "y"), order = c(1, 2))
[1] "2 * cos(x)"

2.2 比較

2.2.1 基本的な構文の違い

\[ \frac{\partial }{\partial x}\sin(x)\bigg|_{x_0} \]

は次のように計算できる：

sym <- derivative(f="sin(x)", var=c(x=0))
num <- derivative(f=function(x) sin(x), var=c(x=0))
出力 
result <- data.frame(
  Method = c("Symbolic", "Numeric"),
  Value = c(sym, num)
)
print(result)
    Method Value
1 Symbolic     1
2  Numeric     1

2.2.2 ４階微分での比較

\[ \frac{\partial ^4}{\partial x^4}\sin(x)\bigg|_{x=0} \] は次のように計算できる：

sym <- derivative(f="sin(x)", var=c(x=0), order=4)
num <- derivative(f=function(x) sin(x), var=c(x=0), order=4)
出力 
result <- data.frame(
  Method = c("Symbolic", "Numeric"),
  Value = c(sym, num)
)
print(result)
    Method         Value
1 Symbolic  0.000000e+00
2  Numeric -9.767766e-12

2.2.3 多変数での比較

\[ \frac{\partial }{\partial x}\frac{\partial ^2}{\partial y^2}y^2\sin(x)\bigg|_{(x,y)=(0,0)} \]

f <- function(x, y) y^2 * sin(x)
derivative(f, var=c(x=0, y=0), order=c(1, 2), accuracy=6)
[1] 2

References

Guidotti, E. (2022). Calculus: High-dimensional numerical and symbolic calculus in r. Journal of Statistical Software, 104(5), 1–37.