Is empirical Bayes an iterative scheme? From the Wikipedia article on the empirical Bayes method (emphasis mine):

In, for example, a two-stage hierarchical Bayes model, observed data $y=\{y_{1},y_{2},\dots ,y_{n}\}$ are assumed to be generated from an unobserved set of parameters $\theta =\{\theta _{1},\theta _{2},\dots ,\theta _{n}\}$ according to a probability distribution $p(y \mid \theta)$. In turn, the parameters $\theta$ can be considered samples drawn from a population characterised by hyperparameters $\eta$, according to a probability distribution $p(\theta \mid \eta )$. In the hierarchical Bayes model, though not in the empirical Bayes approximation, the hyperparameters $\eta$, are considered to be drawn from an unparameterized distribution $p(\eta )$.

However, later in the article, it is written that

Alternatively, the expression can be written as
$$
p(\theta \mid y)=\int p(\theta \mid \eta ,y)p(\eta \mid y)\;d\eta =\int {\frac {p(y\mid \theta )p(\theta \mid \eta )}{p(y\mid \eta )}}p(\eta \mid y)\;d\eta
$$
and the term in the integral can in turn be expressed as
$$
p(\eta \mid y)=\int p(\eta \mid \theta )p(\theta \mid y)\;d\theta
$$
These suggest an iterative scheme, qualitatively similar in structure to a Gibbs sampler, to evolve successively improved approximations to $p(\theta \mid y)$ and $p(\eta \mid y)$. First, calculate an initial approximation to $p(\theta \mid y)$ ignoring the $\eta$ dependence completely; then calculate an approximation to $p(\eta \mid y)$ based upon the initial approximate distribution of $p(\theta \mid y)$; then use this $p(\eta \mid y)$ to update the approximation for $p(\theta \mid y)$; then update $p(\eta \mid y)$; and so on.

My questions are:

*

*If $\eta$ is not considered be drawn from an unparameterized distribution $p(\eta)$ in the empirical Bayes method, then how/why is it treated as a random variable in the expressions above for $p(\theta \mid y)$ and $p(\eta \mid y)$?

*Earlier in the article, it is written that


Using Bayes' theorem,
$$
p(\theta \mid y)={\frac {p(y\mid \theta )p(\theta )}{p(y)}}={\frac {p(y\mid \theta )}{p(y)}}\int p(\theta \mid \eta )p(\eta )\,d\eta
$$
In general, this integral will not be tractable analytically or symbolically and must be evaluated by numerical methods.

However, I do not see how the integrals in the expressions above for $p(\theta \mid y)$ and $p(\eta \mid y)$ are any more tractable.


*The article does not specify how $p(y \mid \theta),p(\theta \mid \eta),p(y \mid \eta),$ and $p(\eta \mid \theta)$ are computed. How are they computed in the empirical Bayes method?


*It is my understanding that the empirical Bayes method essentially involves estimating the (non-random) parameter $\eta$ for the prior on $\theta$ using some samples of $\theta$. Intuitively, I would assume that this would be done once using maximum likelihood estimation, as was done in the book Introduction to Empirical Bayes: Examples from Baseball Statistics by David Robinson. However, in this case, we are not given samples of $\theta$, but we are instead given samples of $y$. How is this dealt with?

Update
Here is an example of an empirical Bayes estimator that more closely fits my understanding. It is part of the Bayes estimator Wikipedia entry.

First, it is correct that, in the (parametric) empirical Bayes approach, the hyperparameter is estimated from the data $y$, by $\hat\eta(y)$, replacing the prior $p(\theta \mid \eta ,y)$ with the empirical "prior" $p(\theta \mid \hat\eta(y) ,y)$ and failing to account for the plugged-in $\hat\eta(y)$. There is therefore no prior distribution on $\eta$ involved in this setting.
The following paragraph in the Wikipedia article is thus puzzling to me as well, as it assumes the existence of an hyperprior $p(\eta)$. It would be more appropriate for the hierarchical Bayes article. In this setting, since
$$p(\theta \mid y)=\int p(\theta \mid \eta ,y)p(\eta \mid y)\;d\eta =\int {\frac {p(y\mid \theta )p(\theta \mid \eta )}{p(y\mid \eta )}}p(\eta \mid y)\;d\eta$$
one can rewrite
$$p(\theta \mid y)=\mathbb E_{\eta|y}\left[\frac {p(y\mid \theta )p(\theta \mid \eta )}{p(y\mid \eta )}\right]$$
and, similarly,
$$p(\eta \mid y)=\mathbb E_{\theta|y}\left[p(\eta|\theta)\right]$$
The numerical approximation would then start from a sample
$$\theta_1,\ldots,\theta_T\sim \pi(\theta|\hat\eta(y),y)$$
then generate
$$\eta_1,\ldots,\eta_T\sim \hat\pi(\eta|y)\propto\sum_{t=1}^T p(\eta|\theta_t)$$
then redraw a new sample
$$\theta_1,\ldots,\theta_T\sim \hat\pi(\theta|y) \propto\sum_{t=1}^T p(y\mid \theta )p(\theta \mid \eta^t )$$
until some stability is achieved in the samples.
This is however a rewording of the Gibbs sampler (which is the special case when $T=1$). The method is thus exact sampling from the hierarchical posterior and unrelated with empirical Bayes, as far as I understand.
Using the alternative representation
$$p(\theta \mid y)={\frac {p(y\mid \theta )p(\theta )}{p(y)}}={\frac {p(y\mid \theta )}{p(y)}}\int p(\theta \mid \eta )p(\eta )\,d\eta$$
makes little sense as well if using the implementation:

*

*generate $$\eta_1,\ldots,\eta_T\sim p(\eta)$$

*generate $$\theta_1,\ldots,\theta_T\sim \hat\pi(\theta|y) \propto p(y|\theta) \sum_{t=1}^T p(\theta|\eta_t)$$
which is not a particularly recommended approximation (as the $\eta_t$'s are simulated from the prior).
