# Understanding the MCMC learning process in a complex statistical model

I am currently reading Rasmussen & Ghahramani's paper on Infinite Mixtures of Gaussian Process Experts. (you can find it at this link)

Unfortunately, I lack experience with understanding how complex statistical models are learnt, especially when multiple techniques are combined. As a result, there are a couple of points which I do not understand about the algorithm specified in the paper.

I've listed out the algorithm's steps at the bottom for convenience. My questions are:

• Step 3 is to "Do Hybrid Monte Carlo". Given that Monte Carlo is sampling, how is this different from doing "Sample", as specified in steps 5 and 6?
• More on steps 5 and 6, do the steps imply that each parameter is sampled once per algorithm iteration (one cycle from steps 2 to 6), or does this imply sampling multiple times and then taking the sample mean to compute the expectation of the parameter? (if yes, why are we taking the expectation in the first place?)
• Step 4 asks to "Optimize" the hyper-hypers - does this imply something like maximum likelihood estimation, and in the first place what exactly is the quantity that we are aiming to optimize here?
• Step 6 has a footnote specifying that the proposal is asymmetric, even though it is a Gaussian proposal (aren't Gaussian proposals symmetric?). The kernel widths are positive, so I suppose the proposal is asymmetric to prevent sampling from outside of the positive reals - but how is this done? Is the "Gaussian proposal" actually supposed to be log normal instead?
• Step 6's footnote also mentions that the Gaussian fit uses the derivative and Hessian of the log posterior - what does this mean?

Your answer is of course very much appreciated even if you are addressing not all of the points above - any clarification at all would already be very helpful!

As mentioned above, the algorithm for learning the data, as outlined by the paper, is:

1. Initialize indicator variables $c_i$ to a single value.
2. Do a Gibbs sampling sweep over all indicators.
3. Do Hybrid Monte Carlo for hyperparameters of the GP covariance function, $v_0; v_1; w_d$, for each expert in turn.
4. Optimize the hyper-hypers, $a$ & $b$, for each of the variance parameters.
5. Sample the Dirichlet process concentration parameter, $\alpha$ using ARS.
6. Sample the gating kernel widths $\phi$, using the Metropolis method to sample from the pseudo-posterior with a Gaussian proposal fit at the current $\phi$.
7. Repeat from 2 until the Markov chain has adequately sampled the posterior.

I've managed to come up with answers of my own to most of the issues. As a result, I've unstruck the issues previously solved, and this answer will contain my findings.

Step 3 is to "Do Hybrid Monte Carlo". Given that Monte Carlo is sampling, how is this different from doing "Sample", as specified in steps 5 and 6?

Hybrid Monte Carlo (aka Hamiltonian Monte Carlo) is, in fact, also sampling. While the wording of "do" vs "sample" initially confused me (leading to the next question), the paper does not write anything to suggest that this procedure is different from the typical Hamiltonian Monte Carlo.

More on steps 5 and 6, do the steps imply that each parameter is sampled once per algorithm iteration (one cycle from steps 2 to 6), or does this imply sampling multiple times and then taking the sample mean to compute the expectation of the parameter? (if yes, why are we taking the expectation in the first place?)

Nowhere in the paper is sampling multiple times per iteration mentioned. As far as I know, the typical approach to MCMC for a hierarchical model is to sample each parameter once, to construct the new state (which comprises of parameters at all levels) of the Markov chain.

Step 6 has a footnote specifying that the proposal is asymmetric, even though it is a Gaussian proposal (aren't Gaussian proposals symmetric?). The kernel widths are positive, so I suppose the proposal is asymmetric to prevent sampling from outside of the positive reals - but how is this done? Is the "Gaussian proposal" actually supposed to be log normal instead?

Upon a closer read of the paper, the gating kernel function is of the form $$K_\phi(x_i,x_i')=\exp\Big(-\frac{1}{2}\sum_{d}(x_{id}-x_{i'd})^2/\phi_d^2\Big)$$ where $x_{id}$ is the $d$-th coordinate of $x_i$, and $\phi_d$ is the d-th coordinate of $\phi$. In step 6, we aim to sample a new value of $\phi$.

Since $\phi_d$ is squared in the kernel, even negative reals are acceptable values (though exploration of the negative reals is essentially meaningless).

The paper does limit exploration of the negative reals by using a log normal; however the log-normal distribution is used as the prior, instead of the proposal. The proposal itself is still Gaussian - the next question answers why it is considered asymmetric.

Step 6's footnote also mentions that the Gaussian fit uses the derivative and Hessian of the log posterior - what does this mean?

This is Hessian-based Metropolis-Hastings, or Stochastic Newton sampling, where the proposal density is: $$q(.|x) = \mathcal{N}(x − H^{-1}(x)g(x), −H^{−1}(x))$$ where $H(x)$ and $g(x)$ are the Hessian and gradient at $x$ respectively. In our case, we are sampling for $\phi$, so we replace $x$ by $\phi$ in the above formula.

The Gaussian proposal for $\phi$ being asymmetric is, thus, because the parameters for the proposal are dependent on the Hessian & gradient of the current state of $\phi$. Hence, while step 6 is written to use the Metropolis method, to be exact, we are using Metropolis-Hastings.