Tag Info

Hot answers tagged

8

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


6

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


4

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.


4

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


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

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


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

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

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


2

If you want to find the global minimum of a function, simulated annealing would be the algorithm to look at, in which case there is no need to view the function as an unnormalised density of any kind and no need to transform the function.


2

MH sampling is used when it's difficult to sample from the target distribution (e.g., when the prior isn't conjugate to the likelihood). So you use a proposal distribution to generate samples and accept/reject them based on the acceptance probability. The Gibbs sampling algorithm is a particular instance of MH where the proposals are always accepted. Gibbs ...


2

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


2

You always need some form of a likelihood. Usually knowledge of the science leads to the likelihood, things like: will the data be discrete or continuous? what is the valid range of the data? what is the likely range of the data? What are reasonable shapes for the data? will contribute to choosing a likelihood function. The likelihood function only need ...


2

You could push the MCMC machinery one step further by simulating the $\mathbf{c}_i$'s conditional on $x_i$ and $w$. Indeed, your equation $$f(x_i|w)=\int_{\Delta(xi)}f(\mathbf{c}_i|w)d\mathbf{c}_i$$ states that $x_i$ is obtained by marginalisation over the $\mathbf{c}_i$'s. In statistical terms, the $\mathbf{c}_i$'s are latent variables: they are not ...


2

You can improve the acceptance rate using delayed rejection as described in Tierney, Mira (1999). It is based on a second proposal function and a second acceptance probability, which guarantees the Markov chain is still reversible with the same invariant distribution: you have to be cautious since "it is easy to construct adaptive methods that might seem to ...


2

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


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

MCMC is a strategy for generating samples $x(i)$ while exploring the state space $X $using a Markov chain mechanism. These are irreducible and aperiodic Markov chains that have $P_{target}(\theta)$ as the invariant distribution. This mechanism is constructed so that the chain spends more time in the most important regions. In particular, it is constructed ...


1

This is an artifact of the stochasticity of MCMC sampling. The best you can do is set a specific seed before each run, using set.seed(): set.seed(123) fls1 = bms(datafls, burn = 50000, iter=100000, g = "BRIC", mprior = "uniform", nmodel = 2000, mcmc="bd", user.int=F) set.seed(123) fls2 = bms(datafls, burn = 50000, iter=100000, g = "BRIC", mprior = ...


1

You should calculate a posterior Credible Interval about your estimate. These are usually calculated from quantiles or using the "Highest Posterior Density (HPD)" region. This is slightly different from the "precision" of the estimate, but I believe you are thinking of this from a frequentist perspective and are sort of missing the forest for the trees. The ...


1

Have you looked at the standard tests for convergence? The Geweke's convergence diagnostic works with a single chain and provides a z-score. There is also the Heidelberger and Welch's convergence diagnostic, which also works with a single chain and gives a p-value of the null hypothesis that the sampled values come from a stationary distribution. If you run ...


1

The physical origin of this algorithm is some way to get intuition -- in that case, we have energies of states instead of probabilities and a system which is in a state $x_t$ and can go or not in a state $x'$. Now, isolated system will always go into a lower energy and never into higher and thus the overall distribution would be trivial -- system is stuck ...


1

As far as for the code, I think there is only one error which is the generation of random numbers according to Laplace distribution. Replace y = laplace(randn(1), b, x(i-1)); by (See the Laplace distribution wiki to learn how to generate random Laplace distributed numbers before seeing my answer.) u = rand()-0.5; uu = b / sqrt(2); y = mu - uu * ...


1

From your question: It is possible that the transition function can suggest a move to an invalid part of the state space, and the jump should be obviously rejected. The jump should not be rejected This is a great article on the problem with immediately rejecting an errant jump. As a consequence, your problem vanishes.


1

An acceptance rate of $1\%$ is extremely low. You have to aim for an acceptance rate of around $15\%-30\%$. In the normal case the optimal rate is the famous $0.234$. However, this is a tricky diagnostic tool since this "optimal rate" can also be achieved despite a terrible sampling (for example when you sample well from an entry and terribly from another ...


1

Whether the proposal should depend on the temperature is problem-specific. If you mean can the proposal depend on the temperature, then the answer is yes. The only important things is that within each temperature, the chain is a convergent Markov chain. Note that it's not only within-chain proposals that can be temperature-dependent. The proposals to move ...



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