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.
 A: 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).
