I am trying do inference on a problem with a posterior distribution with two parameters (not in closed form).
I have built a grid approximation and have been able to create a contour plot of this posterior distribution around reasonable values for the 2 parameters (computed with linear regression). I am also able to generate samples out of this distribution.
I would like to approximate this contour with a normal (from BDA):
$$ p(\theta|y) \approx N(\hat{\theta}, [I(\hat{\theta})]^{-1}) $$
where
$$ I(\theta) = -\frac{d^2}{d\theta^2}\log p(\theta|y) $$
Now, I don't really know how to compute $I(\hat{\theta})$ from my $\theta_{1..n}$ samples. I guess I can use the mean of my samples as the mean for the normal approximation, but I don't know how to compute the variance.
Any ideas how to do this?