53

As detailed in our book with George Casella, Monte Carlo statistical methods, these methods are used to produce samples from a given distribution, with density $f$ say, either to get an idea about this distribution, or to solve an integration or optimisation problem related with $f$. For instance, to find the value of $$\int_{\mathcal{X}} h(x) f(x)\text{d}x\...


27

Firstly, let me note [somewhat pedantically] that There are several different kinds of MCMC algorithms: Metropolis-Hastings, Gibbs, importance/rejection sampling (related). importance and rejection sampling methods are not MCMC algorithms because they are not based on Markov chains. Actually, importance sampling does not produce a sample from the ...


22

The bibliography states that if q is a symmetric distribution the ratio q(x|y)/q(y|x) becomes 1 and the algorithm is called Metropolis. Is that correct? Yes, this is correct. The Metropolis algorithm is a special case of the MH algorithm. What about "Random Walk" Metropolis(-Hastings)? How does it differ from the other two? In a random walk, the ...


16

Well, if you are looking "for any pointers"... The (scaled)(inverse)Wishart distribution is often used because it is conjugate to the multivariate likelihood function and thus simplifies Gibbs sampling. In Stan, which uses Hamiltonian Monte Carlo sampling, there is no restriction for multivariate priors. The recommended approach is the separation strategy ...


14

This is a most interesting question, which relates to the issue of approximating a normalising constant of a density $g$ based on an MCMC output from the same density $g$. (A side remark is that the correct assumption to make is that $g$ is integrable, going to zero at infinity is not sufficient.) In my opinion, the most relevant entry on this topic in ...


13

the main rationale behind using the Metropolis-algorithm lies in the fact that you can use it even when the resulting posterior is unknown. For Gibbs-sampling you have to know the posterior-distributions which you draw variates from.


12

1) You could think about this method as a random walk approach. When the proposal distribution $x \mid x^t \sim N( x^t, \sigma^2)$, it is commonly referred to as the Metropolis Algorithm. If $\sigma^2$ is too small, you will have a high acceptance rate and very slowly explore the target distribution. In fact, if $\sigma^2$ is too small and the ...


11

There seem to be some misconceptions about what the Metropolis-Hastings (MH) algorithm is in your description of the algorithm. First of all, one has to understand that MH is a sampling algorithm. As stated in wikipedia In statistics and in statistical physics, the Metropolis–Hastings algorithm is a Markov chain Monte Carlo (MCMC) method for obtaining a ...


11

A1: Indeed the Gaussian distribution is probably the most used proposal distribution primarily due to ease of use. However, one might want to use other proposal distributions for the following reason Heavy Tails: The Gaussian distribution has light tails. This means that $N(x_{t-1}, \sigma^2)$ will possibly only suggest values between $(x_{t-1} - 3\sigma, ...


11

The acceptance rate depends largely on the proposal distribution. If it has small variance, the ratio of the probabilities between the current point and the proposal will necessarily always be close to 1, giving a high acceptance chance. This is just because the target probability densities we typically work with are locally Lipschitz (a type of smoothness) ...


11

One of the reasons why the original construction of Hamiltonian Monte Carlo can be tricky to understand is that it is more restrictive than necessary, if only to simplify the theoretical proofs. In particular, the negation of the momenta in the deterministic update is indeed practically irrelevant because of the full momenta resampling*. If we include it ...


10

Here you go - three examples. I've made the code much less efficient than it would be in a real application in order to make the logic clearer (I hope.) # We'll assume estimation of a Poisson mean as a function of x x <- runif(100) y <- rpois(100,5*x) # beta = 5 where mean(y[i]) = beta*x[i] # Prior distribution on log(beta): t(5) with mean 2 # (...


10

We don't use MCMC to calculate the $p(\theta | y)$ for each value (or many values) of $\theta$. What MCMC (or the special case of Gibbs sampling) does is generate a (large) random sample from $p(\theta | y)$. Note that $p(\theta | y)$ is not being calculated; you have to do something with that vector (or matrix) of random numbers to estimate $p(\theta)$. ...


10

This type of fine-tuned (Gibbs) MCMC is appropriate for cases when one conditional distribution is most "sticky" than other conditional distributions in the problem. For instance, updating only one [random] part of $\beta$ may be profitable when updating the whole vector results in high rejection rates or in very small moves. An early reference on mixing ...


10

There is a lot of confusion. You want to evaluate the posterior $$ f(\theta|\mathbf{y}) =\frac{f(\mathbf{y}|\theta)f(\theta)}{f(\mathbf{y})} $$ where I use $f()$ to indicate a density, to be as general as possible. You have to decide the likelihood of your model, i.e., $f(\mathbf{y}|\theta)$, and the prior over $\theta$, i.e., $f(\theta)$. In other words, ...


10

The conditional density kernels are: $$\begin{equation} \begin{aligned} f(x|y) &\propto \exp(-|x|-a \cdot |x-y|), \\[6pt] f(y|x) &\propto \exp(-|y|-a \cdot |x-y|). \\[6pt] \end{aligned} \end{equation}$$ The difficulty here is to derive the actual densities that go with these kernels, which takes a bit of algebra. Integration over the full range of ...


9

It is indeed a very poor idea to start learning a topic just from an on-line code with no explanation. Better read a book (like our Introduction to Monte Carlo methods with R!) or an introductory paper and write your own code. As written, this code proposes a random walk on the parameter $(a,b)$ which in your case could be $$(\text{logit}(\mu),\log(\tau))=\...


9

1) the Normal and Uniform are symmetric probability density functions themselves, is this notion of "symmetry" the same as the "symmetry" above? Both distributions are symmetric around their mean. But the symmetry in Metropolis-Hastings is that $q(x|y)=q(y|x)$ which makes the ratio cancel in the Metropolis-Hastings acceptance probability. If one uses a ...


9

$\mathrm{d}q$ is uniform across the entire space and that's the problem! Unfortunately as we consider higher-dimensional spaces out intuition of uniform starts failing us and we end up in conceptual difficulties like this. Yes, the volume of the neighborhood around any given point remains the same size as we increase the dimensionality of our space. But ...


8

1. The problem is not about ergodicity No, this is not related to ergodicity. In the chain without cycling around, one can still move from any island to any other island, and (provided that there are differences in the populations!) the chain is not periodic (because one sometimes stays put), so the chain is ergodic, hence it has a unique stationary ...


8

In step 4, you don't have to reject the proposal $x,\theta$ every time its new likelihood is lower; if you do so, you are doing a sort of optimization instead of sampling from the posterior distribution. Instead, if the proposal is worse then you still accept it with an acceptance probability $a$. With pure Gibbs sampling, the general strategy to sample ...


8

Here is an excerpt from our Monte Carlo Statistical Methods book: 10.3.3. Metropolizing the Gibbs Sampler Hybrid MCMC algorithms are often useful at an elementary level of the simulation process; that is, when some components of the Gibbs sampler conditionals cannot be easily simulated. Rather than looking for a customized algorithm such as Accept--Reject ...


8

An easy example of acceptance probability equal to one is when simulating from the exact target: in that case $$\dfrac{\pi(x')q(x',x)}{\pi(x)q(x,x')}=1\qquad\forall x,x'$$ While this sounds like an unrealistic example, a genuine illustration is the Gibbs sampler, which can be interpreted as a sequence of Metropolis-Hastings steps, all with probability one. ...


8

Disclaimer: although there is nothing to complain about Ben's answer (!), except maybe that the normalising constant of the conditional is not of direct use, here is what I wrote while being off-line, so I may as well post it! The full conditional of $X$ given $Y$ has a density that is proportional to \begin{align} f(x|y) &\propto \exp\{ -|x|-a|y-x|\}\...


8

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 ...


8

Preliminaries: The book Introducing Monte Carlo methods with R (no exclamation mark in the title, even though the Springer book series is called Use R!) was co-authored by my late friend George Casella. Here is a full reproduction of the paragraph containing the quote: The study of independent Metropolis-Hastings algorithms is certainly interesting, but ...


7

The linked algorithm seems to actually be slightly erroneous (I explain the error in the end of the answer), however, the Hastings ratio you ask about is computed correctly Goal of the algorithm The algorithm described is a Metropolis-Hastings sampling on one-dimensional $x$ targeting the unit-exponential distribution. To use this to obtain the samples on ...


7

Are you sure the joint density$$f(x_1,x_2)=\left(\dfrac{x_1}{x_2}\right)\left(\dfrac{\alpha}{x_2}\right)^{x_1-1}\exp\left\{-\left(\dfrac{\alpha}{x_2}\right)^{x_1} > \right\}\mathbb{I}_{\mathbb{R}^*_+}(x_1,x_2)$$ is integrable? When I consider the conditional$$f(x_2|x_1)=\dfrac{1}{{x_2}^{x_1}}\exp\left\{-\dfrac{\gamma}{{x_2}^{x_1}} \right\}$$ it should be ...


7

The slice sampler does not "sample from the log-density". It can, however, use the log density in the calculations to obtain a dependent sequence of observations from the density. The basic idea of a slice sampler is in terms of the density itself, but for various reasons (computational accuracy, primarily) it's usually more convenient to work with the log-...


7

It's not entirely correct to say that inverse methods are impossible to compute. There are perfectly good numerical approximations to the inverse Gaussian CDF. As far as I'm aware, plenty of methods use it to generate gaussian random variables. There are of course plenty of other, possibly simpler methods of generating Gaussians. Concerning rejection ...


Only top voted, non community-wiki answers of a minimum length are eligible