转载

R语言做符号计算

本文作者：黄湘云，2011-2015年在中国矿业大学（北京）的数学与应用数学专业获得学士学位，并从2015年至今在中国矿业大学（北京）统计学专业硕士在读，主要研究方向为复杂数据分析。

引言

谈起符号计算，大家首先想到的可能就是大名鼎鼎的Maple，其次是Mathematica，但是他们都是商业软件，除了其自身昂贵的价格外，对于想知道底层，并做一些修改的极客而言，这些操作也很不可能实现。自从遇到R以后，还是果断脱离商业软件的苦海，R做符号计算固然比不上Maple，但是你真的需要Maple这样的软件去做符号计算吗？我们看看R语言的符号计算能做到什么程度。

符号计算

1.符号微分

在R中能够直接用来符号计算的是表达式，下面以Tetrachoric函数为例，
$$/tau(x)=/frac{(-1)^{j-1}}{/sqrt{j !}}/phi^{(j)}(x)$$
其中
$$/phi(x)=/frac{1}{/sqrt{2/pi}}e^{-/frac{x^2}{2}}$$
在R里，声明表达式对象使用expression函数

NormDensity <- expression(1/sqrt(2 * pi) * exp(-x^2/2)) class(NormDensity) ## [1] "expression"

计算一阶导数

D(NormDensity, "x") ## -(1/sqrt(2 * pi) * (exp(-x^2/2) * (2 * x/2))) deriv(NormDensity, "x") ## expression({ ##      .expr3 <- 1/sqrt(2 * pi) ##      .expr7 <- exp(-x^2/2) ##      .value <- .expr3 * .expr7 ##      .grad <- array(0, c(length(.value), 1L), list(NULL,c("x"))) ##      .grad[, "x"] <- -(.expr3 * (.expr7 * (2 * x/2))) ##      attr(.value, "gradient") <- .grad ##     .value ## })  deriv3(NormDensity, "x")  ## expression({ ##     .expr3 <- 1/sqrt(2 * pi) ##     .expr7 <- exp(-x^2/2) ##     .expr10 <- 2 * x/2 ##     .expr11 <- .expr7 * .expr10 ##     .value <- .expr3 * .expr7 ##     .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) ##     .hessian <- array(0, c(length(.value), 1L, 1L), list(NULL,  ##     c("x"), c("x"))) ##     .grad[, "x"] <- -(.expr3 * .expr11) ##     .hessian[, "x", "x"] <- -(.expr3 * (.expr7 * (2/2) - .expr11   ##     ＊.expr10)) ##     attr(.value, "gradient") <- .grad ##     attr(.value, "hessian") <- .hessian ##     .value ## })

计算 n 阶导数

DD <- function(expr, name, order = 1) {     if (order < 1)         stop("'order' must be >= 1")     if (order == 1)         D(expr, name) else DD(D(expr, name), name, order - 1) DD(NormDensity, "x", 3) ## 1/sqrt(2 * pi) * (exp(-x^2/2) * (2 * x/2) * (2/2) + ((exp(-x^2/2) *  ## (2/2)-exp(-x^2/2)*(2*x/2)*(2*x/2))*(2*x/2)+ ## exp(-x^2/2) * (2 * x/2) * (2/2)))