Investigate hessian instability

[guilhermebodin]

No description provided.

[guilhermebodin]

Currently, we calculate the standard error by the formula

numerical_hessian = Optim.hessian!(func, opt_hyperparameters)
std_err = diag(inv(numerical_hessian))

This is wrong, the correct should be

numerical_hessian = Optim.hessian!(func, opt_hyperparameters)
std_err = sqrt.(diag(inv(numerical_hessian)))

but sometimes the values are negative which leads to a domain error when computing the sqrt.

We should calculate using the outer product of gradients, this is not trivial.
https://discourse.julialang.org/t/inconsistency-between-standard-errors-from-hessian-and-gradient-in-forwarddiff-optim/16769