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 '19 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 '19 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 '19 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."

| cite | improve this answer | |
$\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 '19 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 '19 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 '19 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.