Hot answers tagged

15

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


11

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


10

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$. The most relevant entry on this topic in regard to your suggestion is a paper by Gelfand and Dey (1994, JRSS B), where the authors develop a very similar approach to find$$\int_\mathcal{X} ...


9

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


9

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

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.


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

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


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

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


5

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


5

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


5

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


4

First of all, for your 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 ...


4

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


4

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


4

Given that the Metropolis-Hastings algorithm is based on the ergodic theorem, i.e., on forgetting the initial condition, the way one picks the initial value is of minor importance. In particular, if some information is available about regions of high probability, the starting point may be chosen in one of those regions. According to Metropolis-Hasting ...


4

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


4

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


4

This is a very good question! Although a much less common situation than likelihood underflow. If the likelihood at a given proposed value $\theta'$ is exactly $+\infty$, then the chain should move there and not move except to other values of $\theta$ with infinite likelihood. I do not know of such cases. If the likelihood at a given proposed value ...


4

$P$ has a density against a (reference) measure made of the Lebesgue measure plus the counting measure on $\{0,1\}$. The later measure gives weights of $1$ to the atoms $0$ and $1$. This means that the density at atoms like 0 and 1 is equal to the weight against the counting measure and only the counting measure!, hence it is 1/6 and 2/6 for 0 and 1, ...


4

The MLE estimate $\hat{\alpha}$ could be used as location parameter of a Normal distribution with scale parameter $\sigma$ used as prior distribution of the scaling parameter. Then, such a prior distribution can be updated into a posterior via Metropolis-Hastings algorithm, i.e. a Markov Chain Monte Carlo method used to obtain a sequence of random samples ...


4

No. There are two sorts of asymptotics here: $N$ items of data, and $T$ samples from the posterior. For most priors, as $N$ gets large the posterior distributions should become close*. For large $T$ you simply get a more accurate picture of the posterior under a particular prior for a fixed amount of data. *Asymptotic consistency of Bayesian models has ...


4

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


4

The probability density function (PDF) for the posterior distribution is proportional to the product of the likelihood and the prior distribution: $$ p(\theta|y) \propto p(y|\theta)\,p(\theta), $$ where $y$ denotes the observed data and $\theta$ denotes the unobserved parameter. If the prior is flat, $$ p(\theta) \propto 1, $$ and the likelihood is ...


3

When the (conditional) proposal distributions are continuous, you compute the ratio of the densities, \begin{equation} \frac{q(y\mid x)}{q(x \mid y)} = \frac{g(y)}{g(x)}. \end{equation} Note also that the relevant asymmetry for needing to compute the proposal ratio (instead of plain Metropolis where the proposal ratio is $1$) is whether $q(y\mid x)=q(x\mid ...


3

There are several coding issues: 1.The line of code ββ[i+1,] =rmvt(n=1,delta=ββ[i+1,],ΣΣ=sig,df=3) seems to show that the proposed value ββ[i+1,] is centred at itself: delta=ββ[i+1,] instead of being centred at the previous value delta=ββ[i,] Similar corrections apply to the lines log_prop_curr= log(dmvt(ββ[i,],delta=ββ[i,],ΣΣ=sig,df=3)) # log ...



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