Symbolically calculate the Hessian in (r)yacas

Background: I am using the Ryacas package, trying to figure out a symbolic expression for a the large-sample variance of an MLE estimator.

To do that, I need the (inverse of the) Hessian matrix of the likelihood function. I don't have mathematica (and the online version seems too cumbersome to use to me) and hence I am trying with the Ryacas package which is an interface to the YACAS computer algebra system.

Question: I can't seem to figure out how to compute the Hessian, however. Using the guide here: https://cran.r-project.org/web/packages/Ryacas/vignettes/high-level.html gives me an error for this. Here is a minimal reproducible example (taken from that guide)

L <- yac_symbol("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)")
Hessian(L)
Error in Hessian(L) : could not find function "Hessian"

When I try in another way, namely using the (new?) interface

L <- yac_symbol("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)")
y_fn(L, "HessianMatrix")

I can't seem to get a usable answer either but just HessianMatrix((x^2*y)/4-a*(3*x+(3*y)/2-45))

Would anybody have an idea how to fix this? Would be greatly appreciated!

Thanks

Solution

Here is what you can do with Ryacas v1.1.1:

> library(Ryacas)
> packageVersion("Ryacas")
[1] ‘1.1.1’
> x <- ysym("x")
> y <- ysym("y")
> a <- ysym("a")
> L <- x^2 * (y/4) - a*(3*x + 3*y/2 - 45)
> H <- Hessian(L, c("x", "y", "a"))
> H
{{   y/2,    x/2,     -3},
 {   x/2,      0, (-3)/2},
 {    -3, (-3)/2,      0}} 
> as_r(H)
expression(rbind(c(y/2, x/2, -3), c(x/2, 0, -3/2), c(-3, -3/2, 
    0)))
> eval(as_r(H), list(x = 2, y = 2, a = 2))
     [,1] [,2] [,3]
[1,]    1  1.0 -3.0
[2,]    1  0.0 -1.5
[3,]   -3 -1.5  0.0