I had a look at it. It's not very hard to compute the Hessian, but since newpem uses a non-minimal parameterization the Hessian will be indefinite and interpreting it as an inverse covariance matrix is not possible. One could possible transform the model to a minimal canonical form and compute the covariance of the parameters in that formulation, but it seems a bit messy. One could also add regularization to the parameter vector, but doing this in a nice way is non-trivial. structured_pem would not suffer from this problem as long as the user-chosen parameterization is minimal