symbolic-computation.Rmd

# 符号计算 {#chap-symbolic-computation}

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

16年在统计之都灌水[符号计算与R语言](https://cosx.org/2016/07/r-symbol-calculate)，相应的 Rmd 源文件放在[Github](https://github.com/XiangyunHuang/Symbolic-Computing)上。

```{r}
# 多元函数求偏导
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)
```

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

```{r}
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)
```
Deriv 可以将 R 表达式简化

```{r}
library(Deriv)
Simplify(DD(NormDensity, "x", 3))
```
即 $x (3 - x^2) \mathrm{e}^{-x^2/2}/\sqrt{2 \pi}$，`eval()` 将表达式转为函数，代入数值运算。

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


```{r}
Tetrachoric <- function(x, j) {
  (-1)^(j - 1) / sqrt(factorial(j)) * eval(DD(NormDensity, "x", j))
}
Tetrachoric(2, 3)
```

```{r Tetrachoric, fig.cap="Tetrachoric 函数", echo=FALSE, fig.width=5.5, fig.height=5}
t <- seq(-4, 4, length.out = 80)
plot(c(-4, 4), c(-0.3, 0.3),
  xlab = expression(x), ylab = "",
  main = expression(Tetrachoric(x)), type = "n"
)
abline(v = 0, h = 0, lty = 2)
n <- 8
for (i in 1:n) {
  lines(t, Tetrachoric(t, i),
    lty = 1, lwd = 1.5, type = "l",
    col = terrain.colors(9)[i]
  )
}
legend(-4, .3,
  legend = paste("n=", seq(n)),
  col = terrain.colors(9), lwd = 1
)
```

表达式转函数

```{r}
ExpToFun<-function(x) eval(Simplify(DD(NormDensity, "x", 4)))
ExpToFun(2)
```

函数求积分

```{r}
integrate(ExpToFun, 0, pi)
```
对函数求微分

```{r}
fun <- function(x) x * pi / sqrt(x)
# D(body(fun),'x')
Simplify(D(body(fun), "x"))
```


```{r}
Dfun <- function(x) {
}
body(Dfun) <- Simplify(D(body(fun), "x"))
Dfun
```


```{r}
Dfun(4)
```


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

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

变量替换

```{r}
a <- S("a")
# z 用 a 替换
expr <- subs(expr, z, a)
# y 用 x^2 替换
expr <- subs(expr, y, x^2L)
expr
```

表达式求 2 阶偏导

```{r}
d1_expr <- DD(expr, "x", 2)
expand(d1_expr)
```

求解带参数 $a$ 的一元二次方程

```{r}
solutions <- solve(d1_expr, "x")
solutions
```