# Monte Carlo approximation to find expected value of gradient square

I need to to calculate this term:

$$\mathbb{E}\left[S(Y, L,\theta)S(Y,L,\theta)^\prime\right]$$

Where

$$S(Y,L,\theta) =\frac{\partial}{\partial\theta} l(Y,L,\theta)$$

With $$\theta$$ = maximum likelihood estimates

And:

$$\mathcal{l}(Y, L,\theta) = \sum_{i = 1}^n-\frac{ni}{2}\log(2\pi) -\frac{1}{2}\log(det(\Sigma)) - \frac{1}{2}(L_i - X_i\beta)\Sigma^{-1}(L_i - X_i\beta)^t$$

Each observation $$n$$ is a family composed of $$ni$$ individuals. Y is the observed data, L is the unobserved liability variable distributed as $$N(X\beta, \Sigma)$$, and $$\theta$$ a vector of two parameters ($$\beta, h^2$$), where $$h^2$$ is variance component of $$\Sigma=h^2K + (1-h^2)I$$. Observed data Y are used to define liability variable $$Li$$ based on the prevalence of the disease.
Thus, I need to calculate the expected value of the square gradient of the conditional expectation of complete log-likelihood, i.e. my E-Step in EM algorithm.

I need this therm as it is present in Louis Methid to calculate observed Fisher information.

I've concluded that the only way to estimate this term is to use a stochastic approximation with a Monte Carlo simulation.

So once obtained the MLE estimates, I can generate $$N$$ Monte Carlo samples. But here's where my doubts start: after generating the sample should I calculate the score $$S(Y,L,\theta)$$ using always the original MLE estimate or should I calculate for each sample firstly the MLE and then the score function evaluated at the new MLE? I think the former but I want to be sure.

Thanks.

• There is no visible $h$ or $Y_i$ in your likelihood function, which also seems to be missing a $1/2$ in its last term. The first time involves a $ni/2$ factor that does not seem right either. Commented Sep 5, 2021 at 5:08
• Thanks, I've edited the question to give some clarifications. Commented Sep 5, 2021 at 12:21

## 1 Answer

To be clear, I will outline steps. First, examine a problem for which you know precisely all the parameter values.

You generate random estimates from an expression with the specified parameters and an error term. This can be accomplished by inverting (that is solving the mathematical function) for the error term. If you assume this is from a uniform distribution, then no more work. If exponentially distributed, for example, apply an appropriate transform (same log) to make it uniformly distributed.

Then, generate random deviates by inserting all values into your derived formula that specifies your function of interest.

Use all the data to test your method to uncover (reverse estimate) the true parameter values. As you know the right answer, construct an error estimate to assist in the best estimate path.

More advanced, use select data, here is some advice courtesy of Wikipedia, to quote:

While the naive Monte Carlo works for simple examples, an improvement over deterministic algorithms only can be accomplished with algorithms that use problem-specific sampling distributions. With an appropriate sample distribution it is possible to exploit the fact that almost all higher-dimensional integrands are very localized and only small subspace notably contributes to the integral.[6] A large part of the Monte Carlo literature is dedicated in developing strategies to improve the error estimates. In particular, stratified sampling—dividing the region in sub-domains—and importance sampling—sampling from non-uniform distributions—are two examples of such techniques.

Note: The usual naive Monte Carlo approach is to sample points uniformly over sample space given N uniform random starting deviates. This works, in general, due to the law of large numbers. However, it does not imply necessarily good efficiency especially for non-linear functions of interest.

So, to answer your question use "original MLE estimate or should I calculate for each sample firstly the MLE", this is answered since you know the right answer and can determine which works best over repeated simulations!

• Thanks for your detailed question, now I've read about Monte Carlo integration and the concept is very clear to me, while I have still some doubts about random estimates generation because my function is not trivial. I cannot generate L from a uniform distribution because these are distribute as $N(X\beta, \Sigma)$, thus $L = X\beta + U$ where $U$ is the error term distributed as $N(0, \Sigma)$. So, can I simply generate L from a random normale multivariate distribution and evaluate my function of interest with many MC samples? Commented Sep 5, 2021 at 13:48
• Spreadsheets (like Google Sheets) have an inverse Normal function. Also, FYI an approximation is the sum of 12 Uniforms. Commented Sep 5, 2021 at 21:56
• I've succeded in producing what I needed, thanks a lot for your help. Commented Sep 6, 2021 at 10:39