In different sources there is an algorithm how to calculate the variance of MLE in R. To keep it short:
construct the negative log likelihood function.
minimize it via nlm or optim with hessian=TRUE
invert the Hessian and read out the diagonal components.
I have simulated 1000 Bernoulli random variables and did this algorithm. Then I took the closed formula for MLE and for its variance. Values were almost identical, so yes, this works. But there is a problem, when I come from theoretical point of view. According to asymptotic normality we have \begin{align*} \sqrt{n}(\hat{\theta}_n - \theta_0) \overset{d}{\rightarrow} \mathcal{N}(0,\mathcal{I}^{-1}(\theta_0)) \end{align*} and there for $n$ large enough \begin{align*} \hat{\theta}_n \sim \mathcal{N}(\theta_0,\frac{\mathcal{I}^{-1}(\theta_0)}{n}) \end{align*} I ask myself, why the Hessian matrix after it's inversion isn't divided by $n$ (the number of observations)? If I got it correct, then $\mathcal{H}^{-1}(\hat{\theta})$ is only the approximation for $\mathcal{I}^{-1}(\theta_0)$.