How to get asymptotic covariance matrix when observed information matrix is singular I'm fitting different models by Maximum Likelihood. To do this I'm using a stochastic version of Newton-Raphson algorithm, where both the gradient and the Hessian of the likelihood are estimated at each step. 
In most cases I reach convergence, but then I often encounter the following problem: the estimate Hessian $\hat H$ at convergence is negative definite. This is a problem because I can't invert it to get an estimate of the observed Fisher information.
This happens because one or more parameters are weakly identified. On the other hand other parameters seem to be well identified and their corresponding entries in the estimated Hessian seem to make sense. Identifiability is difficult to assess beforehand for the dynamical models I'm trying to fit.
What I would like is an approach that (starting from the $\hat H$) points out what set of parameters are not identifiable and that gives variances for the remaining parameters.
I tried to do this by tilting the smallest eigenvalue of $\hat H$ in order to get a better conditioning number. The results seem arbitrary: depending on the conditioning number I want to achieve I get different variances for the parameters.
EDIT: 
This is a typical example of an Hessian I can end up with:
H <- matrix(c( 67.23586, 10.477815, 138.696877,
               10.47782, -3.238982,   2.592774,
              138.69688,  2.592774, 473.161347 ), 3, 3, byrow = TRUE)

You see that I have a negative entry in the main diagonal: the second parameter is weakly identifiable, while the other well identified. The other situation I often encounter is that of sub-group of parameters that are redundant - highly correlated.
 A: It is less general than your problem, but I think this article may help: Ridge Estimators in Logistic Regression. Basically, instead of maximizing the (log-)likelihood $L$, you maximize $$ L_\lambda(\beta) = L(\beta) - \lambda \|\beta\|_2^2 .$$ 
This regularization increases all the eigenvalues of the Hessian matrix, so it is not so far away from what you are doing right now.
A: I post here my solution to this problem even though it more of a quick fix than real solution. It seems to work quite well when there are one or more uncorrelated parameters that are not identifiable, not sure how well it does when there is a set of highly correlated parameters.
Start with an initial Hessian $H_0$ and index $k = 0$ and than follow this procedure:


*

*At iteration $k$ verify whether $H_k$ is positive definite "enough" by checking whether the ratio between smallest and biggest eigenvalues is above a certain tolerance. If it is then stop otherwise go to step to 2.

*Increase the small eigenvalues so that the matrix is PD according to the above criterion.

*Invert the modified Hessian in order to get the corresponding covariance matrix $\Sigma_k$.

*Find the parameter with the highest variance (the largest diagonal entry of $\Sigma_k$) and eliminate the corresponding row and column in $H_k$ in order to obtain $H_{k+1}$. Go back to step 1.
In step 4 the idea is that the parameter with the highest variance is the least identifiable. Given that the eigenvectors corresponding to smallest eigenvalues lay mostly in its direction, if we get rid of it the conditioning number of the Hessian should improve.  
In the simple case of the above $3 \times 3$ matrix we get:
CH <- cleanHessian(H)

CH
$hessian
      [,1]     [,2]
[1,]  67.23586 138.6969
[2,] 138.69688 473.1613

$badParam
[1] 2

The function output are a Hessian where the entries corresponding to the unidentifiable parameters have been removed, and a vector containing the indexes of these parameters. In this case the second parameter was not identified. 
I can now invert the reduced Hessian to get the asymptotic covariance matrix: 
solve(CH$hessian)
            [,1]         [,2]
[1,]  0.03762240 -0.011028182
[2,] -0.01102818  0.005346114   

According to my experience with the model I'm fitting here these covariances
are credible.
Here is the R code for the procedure:
 cleanHessian <- function(hessian, tol = 1e-8)
{
  badParam <- numeric()

  repeat{

eDec <- eigen(hessian)

# Find small eigenvalues
lowEigen <- which(eDec$values < eDec$values[1] * tol)

# If there are none, I consider the hessian to be PD and I exit
if( length(lowEigen) == 0 ) break

# Increase the small eigenvalues to a threshold
eDec$values[lowEigen] <- eDec$values[1] * tol

# Invert the modified hessian to get the covariance
COV <- eDec$vectors%*%(t(eDec$vectors) / eDec$value)

# I identify the parameter with the highest variance, I remove the corresponding
# elements from the Hessian and I store its index
bad <- which.max( diag(COV) )
hessian <- hessian[-bad, -bad]
offset <- sum( badParam <= bad )

badParam <- c(badParam, bad + offset)

  }

  return( list("hessian" = hessian, "badParam" = badParam) )
} 

