47

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


22

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


20

I'm not an expert in any of these, but I thought I'd put them out there anyway to see what the community thought. Corrections are welcome. One increasingly popular method, which is not terribly straightforward to implement, is called Hamiltonian Monte Carlo (or sometimes Hybrid Monte Carlo). It uses a physical model with potential and kinetic energy to ...


19

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


13

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


12

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.


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

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


10

In order to get this, and to simplify the matters, I always think first in just one parameter with uniform (long-range) a-priori distribution, so that in this case, the MAP estimate of the parameter is the same as the MLE. However, assume that your likelihood function is complicated enough to have several local maxima. What MCMC does in this example in 1-D ...


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

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


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


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

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


9

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


8

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


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

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

$\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

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


7

There are two basic assumptions that lead to this relationship: The stationary distribution $\pi(\cdot)$ doesn't change too quickly (i.e. it has a bounded first derivative). Most of the probability mass of $\pi(\cdot)$ is concentrated in a relatively small subset of the domain (the distribution is "peaky"). Let's consider the "small $\sigma^2$" case first. ...


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

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


7

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


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