附录 C 符号计算

相比于数值计算,符号计算可以无限精度,包括微分、积分运法,求解线性、非线性方程(组),常微分、偏微分方程(组)等,R 自带几个函数如 deriv()D() 等可以做一些简单的微分运算,扩展包 Ryacas 提供 Yacas 核心计算引擎,symengine 引入 C++ 符号计算库SymEngine,相比于 Ryacas[116]symengine 不会和 Base R 函数冲突。Python 的符号计算模块 sympy [117] 不仅支持简单的四则运算,还支持微分、积分、解方程等,详见官方文档 https://sympy.org/

16年在统计之都灌水符号计算与R语言,相应的 Rmd 源文件放在Github上。

# 多元函数求偏导
ft <- deriv(expression(sin(x1) + sin(x2) + cos(3 * x1 * x2) + x1^2 + x2^2),
  namevec = c("x1", "x2"), function.arg = TRUE
)
# 隐函数求偏导
deriv(y ~ sin(cos(x) * y), namevec = c("x","y"), function.arg = TRUE)
## function (x, y) 
## {
##     .expr1 <- cos(x)
##     .expr2 <- .expr1 * y
##     .expr4 <- cos(.expr2)
##     .value <- sin(.expr2)
##     .grad <- array(0, c(length(.value), 2L), list(NULL, c("x", 
##         "y")))
##     .grad[, "x"] <- -(.expr4 * (sin(x) * y))
##     .grad[, "y"] <- .expr4 * .expr1
##     attr(.value, "gradient") <- .grad
##     .value
## }

下面以标准正态分布的密度函数为例,

NormDensity <- expression(1 / sqrt(2 * pi) * exp(-x^2 / 2))
# 递归的方法求高阶倒数
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)))

Deriv 可以将 R 表达式简化

library(Deriv)
Simplify(DD(NormDensity, "x", 3))
## x * (3 - x^2) * exp(-(x^2/2))/sqrt(2 * pi)

\(x (3 - x^2) \mathrm{e}^{-x^2/2}/\sqrt{2 \pi}\)eval() 将表达式转为函数,代入数值运算。

\[\tau(x) = \frac{(-1)^{j-1}}{\sqrt{j!}}\phi^{(j)}(x)\]

Tetrachoric <- function(x, j) {
  (-1)^(j - 1) / sqrt(factorial(j)) * eval(DD(NormDensity, "x", j))
}
Tetrachoric(2, 3)
## [1] -0.04408344
Tetrachoric 函数

图 C.1: Tetrachoric 函数

表达式转函数

ExpToFun<-function(x) eval(Simplify(DD(NormDensity, "x", 4)))
ExpToFun(2)
## [1] -0.2699548

函数求积分

integrate(ExpToFun, 0, pi)
## -0.06192048 with absolute error < 5.8e-12

对函数求微分

fun <- function(x) x * pi / sqrt(x)
# D(body(fun),'x')
Simplify(D(body(fun), "x"))
## 0.5 * pi/sqrt(x)
Dfun <- function(x) {
}
body(Dfun) <- Simplify(D(body(fun), "x"))
Dfun
## function (x) 
## 0.5 * pi/sqrt(x)
Dfun(4)
## [1] 0.7853982

下面简单介绍 symengine 的符号计算能力

library(symengine)
# 声明几个符号变量
use_vars(x, y, z)
# 表达式展开
expr <- (x + y + z) ^ 2L - 42L
expand(expr)
## (Add)    -42 + 2*x*y + 2*x*z + 2*y*z + x^2 + y^2 + z^2

变量替换

a <- S("a")
# z 用 a 替换
expr <- subs(expr, z, a)
# y 用 x^2 替换
expr <- subs(expr, y, x^2L)
expr
## (Add)    -42 + (a + x + x^2)^2

表达式求 2 阶偏导

d1_expr <- DD(expr, "x", 2)
expand(d1_expr)
## (Add)    2 + 4*a + 12*x + 12*x^2

求解带参数 \(a\) 的一元二次方程

solutions <- solve(d1_expr, "x")
solutions
## VecBasic of length 2
## V( -1/2 + (-1/2)*sqrt(1 + (-1/3)*(2 + 4*a)), -1/2 + (1/2)*sqrt(1 + (-1/3)*(2 + 4*a)) )

参考文献

[116]
M. M. Andersen and S. Højsgaard, Ryacas: A computer algebra system in R,” Journal of Open Source Software, vol. 4, no. 42, 2019,Available: https://doi.org/10.21105/joss.01763
[117]
A. Meurer et al., SymPy: Symbolic computing in python,” PeerJ Computer Science, vol. 3, p. e103, Jan. 2017, doi: 10.7717/peerj-cs.103.