0

I want to obtain the Hessian of the following function:

def llik_scalars(param_vector, *args):
    Fsc = param_vector[0]
    Qsc = param_vector[1]
    Rsc = param_vector[2]

    y = args[0]
    burnin = args[1]

    F = np.matrix(Fsc)
    Q = np.matrix(Qsc)
    R = np.matrix(Rsc)

    predstate, predp, _, _ = kalmanfilter(F=F, Q=Q, R=R, y=y, plot = False)
    T = len(predp)
    predstate = np.array([predstate[t].item() for t in range(len(predstate))])
    predp = np.array([predp[t].item() for t in range(len(predp))])

    Sigmat = predp + Rsc
    Mut = predstate

    LL = 0
    for t in range(burnin, T):
        exponent = -0.5 * (y[t]-Mut[t])**2 / Sigmat[t]
        cc = 1 / math.sqrt(2*math.pi*Sigmat[t])
        LL -= math.log(cc*math.exp(exponent))
    return LL

I am trying to do this with the Hessian function of the numdifftools package. In the documentation, I found the following information. If you want for instance the hessian of the rosenbrock function, which is defined as Rosen, The hessian is calculated in the following way:

> H = nd.Hessian(rosen)([1, 1])

Where the Hessian is calculated in the point [1,1]

Following the documentation, it should be possible to give arguments to the Hessian function:

class Hessian(f, step=None, method=’central’, order=2, full_output=False, **step_options)
Parameters
fun [function] function of one array fun(x, *args, **kwds)

I tried this in the following way:

hess = nd.Hessian(kf.llik_scalars(themin.x, (y,burnin)))(themin.x)

themin.x is the point where i want to evaluate the Hessian.

themin.x
Out[49]: array([0.67605231, 0.7457089 , 0.72205726])

The error I get when I run the above code:

burnin = args[1]

IndexError: tuple index out of range

I dont understand how the tuple is out of range

Rens
  • 177
  • 9

1 Answers1

0

It is a bit difficult to debug your call without having the whole functions and parameters, that's why I just made a simple example using the rosen function.

I added two faked additional parameters and named the new function rosen_ext.

import numdifftools as nd import numpy as np

def rosen(x):
    return (1-x[0])**2 + 105.*(x[1]-x[0]**2)**2

def rosen_ext(x, *args):
    y = args[0]
    burnin = args[1]

    print(y)
    print(burnin)

    return (1-x[0])**2 + 105.*(x[1]-x[0]**2)**2

I could reproduce your error when calling the Hessian for rosen_ext in your way.

The Hessian for rosen in x = [2, 2] should be

[[4202. -840.]
 [-840.  210.]]

To get this result for rosen_ext you need to call it like this:

# call of the Hessian for the modified rosen in some point X = [2, 2]
H = nd.Hessian(rosen_ext)([2, 2], 2, 3)
print(H)

It returns a bunch of 2 and 3 (helps to see if the additional parameters were put in right way) and the correct answer.

2
3
...
2
3
2
3
2
3
[[4202. -840.]
 [-840.  210.]]

Using this approach I would assume your call should be like this:

hess = nd.Hessian(kf.llik_scalars)(themin.x, y, burnin)

I cannot debug it, so could you please check if it works?

Anton
  • 4,544
  • 2
  • 25
  • 31