# Tag Info

13

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

8

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

7

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

7

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

6

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

6

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

6

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

6

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

5

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

5

I think that this paper from Heikki Haario et al. will give you the answer you need. The markovianity of the chain is affected by the adaptation of the proposal density, because then a new proposed value depends not only of the previous one but on the whole chain. But it seems that the sequence has still the good properties if great care is taken.

5

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

4

If your Gibbs sampler produces a non-irreducible Markov chain, it means there are parts of the space with positive mass under the target distribution that are never visited by the chain (for some starting values). In order to produce an irreducible chain, you need adding moves to the original Gibbs sampler, for instance by using re-parametrisations that link ...

4

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

3

I think the order is correct, but the labels assigned to p(x) and p(y|x) were wrong. The original problem states p(y|x) is log-normal and p(x) is Singh-Maddala. So, it's Generate an X from a Singh-Maddala, and generate a Y from a log-normal having a mean which is a fraction of the generated X.

3

The scheme you propose is strictly ascending the probability gradient. Think what would happen: you would end up at $x$ such that $P(x)$ is maximal and you would stay there for ever. In the end, the sample you generate would consist mostly of $x$, whatever the size of the sample, so the distribution of that sample would not be close to $P$. Also bear in ...

3

You are rather looking for a simulated annealing, which is easier to understand when formulated in the original, physics way: Having $x$ is a state of the system $f(x)$ is an energy of the system; energy is defined up to addition of a constant, so there is no problem with it being negative or positive -- the only constraint is the-lower-the-better $T$ is a ...

3

I don't have references on hand, unfortunately, but generally the more dependence there is between $X$ and $Y$, the more efficient it will be to update them jointly. In your example $X$ and $Y$ are independent; if, for example, $X$ and $Y$ were still jointly Gaussian but had a correlation of .99, updating them separately would be much worse. Updating them ...

3

Here is an R code in the univariate case for the above Metropolis-within-Gibbs approach drafted by @alberto. No indication of the chain getting stuck: the acceptance rate for the $x$ component is close to 50%. First, I picked some pseudo-values to run the algorithm: #observation from N(x,1) y=3.081927 #latent x from t(nu,theta,1) nu=3 Second, I simulated ...

3

I don't have a great example off the top of my head, but MH is easy compared to direct sampling whenever the parameter's prior is not conjugate with that parameter's likelihood. In fact this is the only reason I have ever seen MH preferred. A toy example is that $p \sim \text{Beta}(\alpha, \beta)$, and you wanted to have (independent) priors $\alpha, \beta ... 3 First of all, for you model to be hierarchical, you need hyperpriors for$\alpha$and$\beta\$ (as already explained by Procrastinator). For the sake of simplicity, lets assume uniform priors on the positive part of the real axis. So that have a hierarchical model as follows: $$y_{i}| \lambda_{i}\sim Poisson(\lambda_{i})$$ $$\lambda_{i}|\alpha, \beta \sim ... 3 You do not need the \alpha since it is a parameter. The change of variables formula applies to the variable with respect to which you are "integrating". It is x in your case. So MH is right to demand that you remove the excess factor. So what you really have is:$$ p(X|\alpha) = \frac{\exp(-\exp(\alpha))\exp(\alpha x)}{x!} $$had you applied some ... 3 When doing M-H it's a good idea to propose from a distribution that is "similar" to the target distribution. Hence, you may try to use the generated samples from a previous simulation as proposals for the simulation of a new posterior that you consider to be "close / similar" to the previous one. You may sample with replacement the previous parameters ... 2 It will give valid samples, provided you include the correct Hastings factor in the acceptance ratio. But using the inverse Hessian directly as the variance can be dangerous, and some additional corrections are recommended. For more details, see the paper Hessian-based Markov Chain Monte-Carlo Algorithms. The algorithm presented there is a bit different ... 2 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). ... 2 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. best regards 2 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 ... 2 The multivariate normal is symmetric insofar as it applies to being a proposal distribution. One can just examine the ratio of PDFs:$$ \begin{eqnarray} \frac{\mathcal{N}(X \mid Y, \Sigma)}{\mathcal{N}(Y \mid X, \Sigma)} &=& \frac{(2\pi)^{-p/2} |\Sigma|^{-1/2} \exp\left( -\frac{1}{2} (X - Y)'\Sigma^{-1}(X-Y) \right)}{(2\pi)^{-p/2} |\Sigma|^{-1/2} ...

2

If you want to extend the example you link to to a multivariate regression, take the code as it is and: Add one more predictor in the code chunk generating the data Add one more parameter in likelihood, as in pred = a1*x1 + a2*x2 + b Add the additional parameter in the prior specification Adjust the MCMC and plots to deal with 4 instead of 3 parameters

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