3
$\begingroup$

I'm thoroughly confused by the Metropolis Algorithm as defined in Casella and Berger's Statistical Inference. Namely, here's the definition (p.254):

Let $Y \sim f_Y(y)$ and $V \sim f_V(v)$, where $f_Y$ and $f_V$ have common support. To generate $Y \sim f_Y(y):$

  1. Generate $V \sim f_V$. Set $Z_0=V$.

Then for $i=1,2,...:$

  1. Generate $U_i \sim \text{uniform}(0,1)$, $V_i \sim f_V$, and calculate $\rho_i=\text{min}(\frac{f_Y(V_i)}{f_V(V_i)} \frac{f_V(Z_{i-1})}{f_Y(Z_{i-1})},1).$

  2. Set $$Z_i =\begin{cases} V_i & \text{if }U_i\leq \rho_i \\ Z_{i-1} & \text{if }U_i > \rho_i \end{cases}.$$

Then $Z_i \to Y$ in distribution as $i \to \infty$.

But with such a definition $\rho_1 = 1$ and $Z_1 = V$, so that all $\rho_i=1$ and all $Z_i = V$. Am missing something or is there an error in the definition?

Is it that even though $V_1$ and $Z_0$ have the same distribution, you're drawing twice and then not necessarily $v_1 = z_0$?

$\endgroup$
  • 1
    $\begingroup$ Why do you think that $p_i = 1$ for all $i$? $\endgroup$ – jbowman Mar 10 at 18:36
  • $\begingroup$ @jbowman, see my last paragraph. Is that's how it's supposed to be understood? Because then I don't think that $\rho_i = 1$ for all $i$ anymore. But if not, then I thought the two terms would just be reciprocals of one another, because $V_i = Z_{i-1} = V$. $\endgroup$ – Ryker Mar 10 at 20:34
  • 1
    $\begingroup$ Your last paragraph does not appear to say anything about why you think $p_i = 1$ for all $i$.... $\endgroup$ – jbowman Mar 11 at 2:01
8
$\begingroup$

The confusion stems from a misunderstanding of the notation $$V \sim f_V$$ which means both (a) $V$ is a random variable with density $f_V$ and (b) $V$ is created by a PRNG algorithm that reproduces a generation of a random variable with density $f_V$. Each time a generation $V_i\sim f_V$ occurs in the algorithm from Casella and Berger, a new realisation of a random variable with density $f_V$ occurs, which is independent from all previous realisations, hence different from these previous realisations. (Setting aside discussions about the pseudo-randomness of these generators and the resulting approximation in the independence statement.) Equivalently, stating that the $𝑉_𝑖$'s are all identically distributed from the same distribution $𝑓_𝑉$ does not mean that their realisations all are numerically identical.

The starting point of the Metropolis-Hastings algorithm is arbitrary, either fixed $Z_0=0$ for instance or random, for instance $Z_0\sim f_V$ [a notation meaning that $Z_0$ is distributed from $f_V$]. This starting value is always accepted. For $i=1$, one generates $V_1\sim f_V$ [meaning that $V_1$ is distributed from $f_V$, independently and thus different from $Z_0$] $$Z_1 =\begin{cases} V_1 & \text{if }U_1\leq \rho_1=\min\left(\frac{f_Y(V_1)}{f_V(V_1)} \frac{f_V(Z_{0})}{f_Y(Z_{0})},1\right) \\ Z_{0} & \text{if }U_1 > \rho_1 \end{cases}$$ and $\rho_1\ne 1$ in general. Hence sometimes $V_1$ is accepted and sometimes not. The same applies to the following steps.

To make a toy illustration on how the algorithm applies, take $f_V$ to be the density of a $\mathcal N(0,1)$ distribution and $f_Y$ to be the density of a $\mathcal N(1,1)$ distribution. A sequence of iid generations from $f_V$ is for instance [by a call to R nrorm] $$0.45735433,-0.99178415,-1.08312586,-0.85762451,0.92186197,-0.50442298,...$$ [note that they are all different] and a sequence of generations from $\mathcal U$ is for instance [by a call to R nunif] $$0.441328,0.987837,0.386258,0.316593,0.195910,0.2772669,...$$ [note that they are all different]. Applying the algorithm with starting value $Z_0=0$ means considering $$\frac{f_Y(V_1)}{f_V(V_1)} \frac{f_V(Z_{0})}{f_Y(Z_{0})}=0.9582509\big/ 0.6065307=1.579889>1$$ which implies that $Z_1=V_1=0.45735433$. Then $$\frac{f_Y(V_2)}{f_V(V_2)} \frac{f_V(Z_{1})}{f_Y(Z_{1})}= 0.2249709 \big/ 0.9582509 = 0.2347724 < U_2=0.987837$$ which implies that $Z_2=Z_1$. The algorithm can be applied step by step to the sequences provided above, which leads to \begin{align*} \frac{f_Y(V_3)}{f_V(V_3)} \frac{f_V(Z_{2})}{f_Y(Z_{2})}&=0.2053581 \big/0.9582509 = 0.2143051 < U_3 \qquad Z_3 = Z_1\\ \frac{f_Y(V_4)}{f_V(V_4)} \frac{f_V(Z_{3})}{f_Y(Z_{3})}&=0.2572712 \big/0.9582509 = 0.2684800 < U_4 \qquad Z_4 = Z_1\\ \frac{f_Y(V_5)}{f_V(V_5)} \frac{f_V(Z_{4})}{f_Y(Z_{4})}&=1.5247980 \big/0.9582509 = 1.591230 > 1 \qquad\quad Z_5 = V_5\\ &\qquad\qquad\vdots\\ \end{align*} producing a sequence as, e.g., below (notice the flat episodes in the graph, which correspond to a sequence of rejections).

$\qquad\qquad\qquad$enter image description here

As a last remark, the only potentially confusing part in the description of Casella and Berger is the very first sentence where the random variables $Y$ and $V$ are not needed. It could have been clearer to state "Let $f_Y$ and $f_V$ be two densities with common support."

$\endgroup$
  • $\begingroup$ What's unclear about it? The definition above states that $Z_0 =V \sim f_V$ and that all subsequent $V_i \ \sim f_V$ also. This is the definition given in the book and I was or am confused about how the algorithm works given that. $\endgroup$ – Ryker Mar 10 at 18:10
  • $\begingroup$ At the very end, when you say that $Z_2 = Z_1$ and so on, you mean that $Z_1 = Z_2 = Z_3 = \dots$? $\endgroup$ – Ryker Mar 11 at 12:42
  • 1
    $\begingroup$ Not at all: I mean that you can take over and repeatedly apply the same steps for the subsequent values of $V_3$, $V_4$, ..., to see when the next value of $V_i$ is accepted and when it isn't. $\endgroup$ – Xi'an Mar 11 at 12:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.