0

I am trying to create a function for theoretical hessian matrix that I can then evaluate at different locations. First I tried setting expressions as values in a matrix or array, but although I could initially set an expression into a matrix I couldn't replace with the value calculated.

hessian_matrix <- function(gx, respect_to){
out_mat <- matrix(0, nrow=length(respect_to), ncol=length(respect_to))
  for(i in 1:length(respect_to)){
    for(j in 1:length(respect_to)){
      dthetad2x <- deriv(D(gx, respect_to[i]), respect_to[j], function.arg=TRUE)
      # also tried
      # dthetad2x <- as.expression(D(D(gx, respect_to[i])))
      out_mat[i,j] <- dthetad2x
    }
return(out_mat)
}

Because that didn't work, I decided to create an environment to house the indeces of the hessian matrix as object.

hessian_matrix <- function(gx, respect_to){
  out_env <- new.env()
  for(i in 1:length(respect_to)){
    for(j in 1:length(respect_to)){
      dthetad2x <- as.call(D(D(gx, respect_to[i]), respect_to[j]))
      assign(paste0(i,j), dthetad2x, out_env)
    }
  }
  return(out_env)
}

g <- expression(x^3-2*x*y-y^6)
h_g <- hessian_matrix(g, respect_to = c('x', 'y'))

This worked, and when I pass this in as a parameter to evaluate I can see the expression, but I can't evaluate it. I tried with call(), eval(), do.call(), get(), etc. and it didn't work. I also assigning the answer within the environment passed, making a new environment to return, or simply using variables.

fisher_observed <- function(h, at_val, params, sum=TRUE){
  out_env <- new.env()
  # add params to passed environment
  for(i in 1:length(at_val)){
    h[[names(at_val)[i]]] <- unname(at_val[i])
  }
  for(i in ls(h)){
    value <- do.call(i, envir=h, at_val)
    assign(i, value, out_env)
  }
  return(h)
}

fisher_observed(h_g, at_val=list(x=1,y=2))

According the code for do.call() this is how it should be used, but it isn't working when passed as a parameter in this way.

3 Answers3

1

I think this (almost) does what you're looking for:

fisher_observed <- function(h, at_val) {
  values <- numeric(length = length(names(h)))
  for (i in seq_len(length(names(h)))) {
    values[i] = purrr::pmap(.l = at_val, function(x, y) eval(h[[names(h)[i]]]))
  }
  names(values) = names(h)
  return(values)
}

This currently returns a named list of evaluated points:

$`21`
[1] -2

$`22`
[1] -480

$`11`
[1] 6

$`12`
[1] -2

you'd still need to re-arrange this into a matrix (should be fairly easy given the column names are preserved. I think the key thing is that the names must be characters when looking up values in h_g.

nrennie
  • 1,877
  • 1
  • 4
  • 14
1

You cannot have a matrix of "calls" but you can have a character matrix then evaluate it:

hessian_matrix <- function(gx, respect_to){
  out_mat <- matrix("", nrow=length(respect_to), ncol=length(respect_to))
  for(i in 1:length(respect_to)){
    for(j in 1:length(respect_to)){
      dthetad2x <- D(D(gx, respect_to[i]), respect_to[j])
      out_mat[i,j] <- deparse(dthetad2x)
    }
  }
  return(out_mat)
}

g <- expression(x^3-2*x*y-y^6)
h_g <- hessian_matrix(g, respect_to = c('x', 'y'))
h_g
#>      [,1]          [,2]              
#> [1,] "3 * (2 * x)" "-2"              
#> [2,] "-2"          "-(6 * (5 * y^4))"

apply(h_g,  1:2, \(x) eval(str2lang(x), list(x=1, y=2)))
#>      [,1] [,2]
#> [1,]    6   -2
#> [2,]   -2 -480

Ric
  • 5,362
  • 1
  • 10
  • 23
  • Thank you Ric - this is wonderful as it most closely does what I was asking. I struggled whether to check your solution or the one above as they're both excellent and I will likely merge them together. – Casey Jayne Feb 07 '23 at 01:06
1

R already has the hessian matrix function. You do not have to write one. You could use deriv or deriv3 as shown below:

g <- expression(x^3 - 2 * x * y - y^6)
eval(deriv3(g, c('x','y')),list(x=1,y=2))

[1] -67
attr(,"gradient")
      x    y
[1,] -1 -194
attr(,"hessian")
, , x

     x  y
[1,] 6 -2

, , y

      x    y
[1,] -2 -480

If you want to use a function, you could do:

hessian <- function(expr,values){
  nms <- names(values)
  f <- eval(deriv3(g, nms),as.list(values))
  matrix(attr(f, 'hessian'), length(values), dimnames = list(nms,nms))
}

hessian(g, c(x=1,y=2))
   x    y
x  6   -2
y -2 -480

Although the function is not necessary as you would do double computation in case you wanted the gradient and hessian

Onyambu
  • 67,392
  • 3
  • 24
  • 53