1

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

ashwin agrawal
  • 1,603
  • 8
  • 16
Jean_N
  • 489
  • 1
  • 4
  • 19

3 Answers3

2
library(Ryacas)

yac_str("HessianMatrix(x^2 * (y/4) - a*(3*x + 3*y/2 - 45), {x,y})")
# "{{y/2,x/2},{x/2,0}}"

yac_str("PrettyForm(HessianMatrix(x^2 * (y/4) - a*(3*x + 3*y/2 - 45), {x,y}))") %>% cat
# /              \
# | / y \ / x \  |
# | | - | | - |  |
# | \ 2 / \ 2 /  |
# |              |
# | / x \ ( 0 )  |
# | | - |        |
# | \ 2 /        |
# \              /

Ref: yacas.readthedocs

Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225
  • Thank you,Stephane, this worked for the toy example! Unfortunately when I use my real function (that's rather more complcited) I get a lot of "Deriv" expressions in in the Hessian that do not get evaluated. PrettyForm is not useful in that case as I want to process further the output and subtract from it another expression for the large sample variance of another estimator (trying to show that one estimator is more efficient than another for all values of the true (multinomial) underlying distribution). Thank you, I mark your answer as correct, it gave me some pointers on how to proceed! – Jean_N Oct 28 '19 at 09:11
  • @Jean_N: Thanks, good to know. I will also just mention the [caracas](https://cran.r-project.org/package=caracas) package for you; this is based on SymPy rather than yacas. – mikldk Jun 06 '20 at 08:53
1

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
mikldk
  • 194
  • 1
  • 4
  • I came back to this question much later as I re-started working on the problem. In the meantime I had to switch to nlsDeriv to help me out with my problem, but your solution works and is exactly what I needed, thank you! I hope to work more with Ryacas now and get to know it better. – Jean_N Jun 04 '20 at 20:33
0

The error, Error in Hessian(L) : could not find function "Hessian", is due the fact that Hessian function is not part of Ryacas package, instead it is part of numDeriv package in R.

Also below is the documentation of Hessian function in R:

https://www.rdocumentation.org/packages/numDeriv/versions/2016.8-1.1/topics/hessian

I hope this solves your problem. So please install numDeriv package first and then use the Hessian function.

I also tried to use the above example, below are the results:

library(Ryacas)
library(numDeriv)
L <- yac_symbol("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)")
Hessian(L)

Output:

{}

The output is empty set, but I guess you can refer to the documentation of Hessian function and it will definitely help.

ashwin agrawal
  • 1,603
  • 8
  • 16
  • Thank you, but I am afraid this is not what I am looking for. Here is why: if you scroll down to the bottom of the link to the Ryacas documentation, you will see the "Hessian" function. Morever, the link to the function you give shows that the function is differently named (with a small h and not big), and, lastly, the hessian function from the numDeriv package gives a numerical approximation - not what I am looking for here (as the title states, i am looking for a symbolic calculation). – Jean_N Oct 28 '19 at 06:40