caracas: Computer Algebra in R via SymPy

Quick start

x <- symbol('x')
eq <- 2*x^2 - x
eq
## [caracas]:    2    
##            2⋅x  - x
as.character(eq)
## [1] "2*x^2 - x"
as_r(eq)
## expression(2 * x^2 - x)
tex(eq) # $$`r tex(eq)`$$
## [1] "2 x^{2} - x"

\[2 x^{2} – x\]

solve_sys(eq, x)
## Solution 1:
##   x =  0 
## Solution 2:
##   x =  1/2
der(eq, x)
## [caracas]: 4⋅x - 1
subs(eq, x, "y")
## [caracas]:    2    
##            2⋅y  - y

Linear algebra

A <- matrix(c("x", 2, 0, "2*x"), 2, 2)
B <- as_symbol(A)
B
## [caracas]: ⎡x   0 ⎤
##            ⎢      ⎥
##            ⎣2  2⋅x⎦
determinant(B)
## [caracas]:    2
##            2⋅x
Binv <- inv(B) # or solve_lin(B)
Binv
## [caracas]: ⎡ 1      ⎤
##            ⎢ ─    0 ⎥
##            ⎢ x      ⎥
##            ⎢        ⎥
##            ⎢-1    1 ⎥
##            ⎢───  ───⎥
##            ⎢  2  2⋅x⎥
##            ⎣ x      ⎦
tex(Binv)
## [1] "\\left[\\begin{matrix}\\frac{1}{x} & 0\\\\- \\frac{1}{x^{2}} & \\frac{1}{2 x}\\end{matrix}\\right]"

\[\left[\begin{matrix}\frac{1}{x} & 0\\- \frac{1}{x^{2}} & \frac{1}{2 x}\end{matrix}\right]\]

eigen_val(Binv)
## [[1]]
## [[1]]$eigval
## [caracas]: 1
##            ─
##            x
## 
## [[1]]$eigmult
## [1] 1
## 
## 
## [[2]]
## [[2]]$eigval
## [caracas]:  1 
##            ───
##            2⋅x
## 
## [[2]]$eigmult
## [1] 1

Matrix-vector multiplication with %*%; subsetting with [, diag() etc. also works.

Maximising the multinomial likelihood

p <- as_symbol(paste0("p", 1:3))
y <- as_symbol(paste0("y", 1:3))
a <- as_symbol("a")
l <- sum(y*log(p))
l
## [caracas]: y₁⋅log(p₁) + y₂⋅log(p₂) + y₃⋅log(p₃)
L <- -l + a*(sum(p) - 1)
L
## [caracas]: a⋅(p₁ + p₂ + p₃ - 1) - y₁⋅log(p₁) - y₂⋅log(p₂) - y₃⋅log(p₃)
g <- der(L, c(p, a))
g
## [caracas]: ⎡    y₁      y₂      y₃                  ⎤
##            ⎢a - ──  a - ──  a - ──  p₁ + p₂ + p₃ - 1⎥
##            ⎣    p₁      p₂      p₃                  ⎦
sol <- solve_sys(g, c(p, a))
sol
## Solution 1:
##   p1 =       y₁     
##        ────────────
##        y₁ + y₂ + y₃ 
##   p2 =       y₂     
##        ────────────
##        y₁ + y₂ + y₃ 
##   p3 =       y₃     
##        ────────────
##        y₁ + y₂ + y₃ 
##   a  =  y₁ + y₂ + y₃
sol[[1L]]$p1
## [caracas]:      y₁     
##            ────────────
##            y₁ + y₂ + y₃
tex(sol[[1L]]$p1)
## [1] "\\frac{y_{1}}{y_{1} + y_{2} + y_{3}}"

caracas objects can be turned into R objects using as_r():

l
## [caracas]: y₁⋅log(p₁) + y₂⋅log(p₂) + y₃⋅log(p₃)
l_e <- as_r(l)
l_e
## expression(y1 * log(p1) + y2 * log(p2) + y3 * log(p3))
eval(l_e, list(p1 = 0.2, p2 = 0.3, p3 = 0.5, 
               y1 = 18, y2 = 31, y3 = 51))
## [1] -101.6435

More information and examples

Please refer to one or more of the following:

• https://cran.r-project.org/package=caracas
• https://github.com/r-cas/caracas/
• https://r-cas.github.io/caracas/