Episode #125 of the Stack Overflow podcast is here. We talk Tilde Club and mechanical keyboards. Listen now

Hot answers tagged

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


21

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


18

Problem Suppose $Y \sim \text{N}(\text{mean} = \mu, \text{Var} = \frac{1}{\tau})$. Based on a sample, obtain the posterior distributions of $\mu$ and $\tau$ using the Gibbs sampler. Notation $ \mu$ = population mean $ \tau$ = population precision (1 / variance) $n$ = sample size $\bar{y}$ = sample mean $s^2$ = sample variance Gibbs sampler [...


15

You tried to show detailed balance for the Markov chain that is obtained by considering one transition of the Markov chain to be the 'Gibbs sweep' where you sample each component in turn from its conditional distribution. For this chain, detailed balance is not satisfied. The point is rather that each sampling of a particular component from its conditional ...


15

The algorithm that is now called Gibbs sampling forms a Markov-chain and uses Monte-Carlo simulation for its inputs, so it does indeed fall within the proper scope of MCMC (Markov-Chain Monte-Carlo) methods. Historically, the method can be traced back at least to the mid-twentieth century, but it was not well-known and was only later popularised by the ...


13

For those who find this question down the road: there's now also Stan. Stan may one day replace OpenBUGS and JAGS, but it does not yet support all the analyses that these other packages do.


12

I find this document GIBBS SAMPLING FOR THE UNINITIATED by Resnik & Hardisty very useful for non-statistics background folks. It explains why & how to use Gibbs sampling, and has examples demonstrating the algo. Seems I cannot comment yet. Gibbs sampling is not a self-contained concept. It requires some prerequisite knowledge. Below is the ...


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

Look at this case first. Dropping terms that do not depend on $x_1$, we have. $$ \pi(x_1\mid x_2,\dots,x_d) = \frac{\pi(x_1,x_2,\dots,x_d)}{\pi(x_2,\dots,x_d)} \propto e^{x_1 x_2} $$ $$ P(X_1=-1\mid X_2 = x_2, \dots, X_n=x_n) = \frac{e^{-x_2}}{C} $$ $$ P(X_1=1\mid X_2 = x_2, \dots, X_n=x_n) = \frac{e^{x_2}}{C} $$ $$ \frac{e^{-x_2}}{C} + \frac{e^{x_2}}...


11

This characterisation of the joint by the conditionals is the Hammersley-Clifford theorem. Unpublished by the authors but later established under the positivity condition by Julian Besag in 1974. The positivity condition is essentially the constraint that the support of the joint is equal to the product of the supports of the conditionals. In that case, for ...


10

Ok, what you need to do is compute the joint posterior up to a constant, i.e. $f(y_1,...,y_n|\mu,\tau)p(\mu,\tau)$. Then to compute the conditional posterior $\pi(\mu|\mathbf{y},\tau)$ you just treat the $\tau$ terms as fixed and known, so that some of them can be cancelled out. Then you do the same thing with $\mu$ to get $\pi(\tau|\mathbf{y},\mu)$. So ...


10

This discrepancy arises because there are two different parameterizations of the Gamma distribution and each relate differently to the Inverse Gamma distribution. On Wikipedia, the two parameterizations for the Gamma distribution are differentiated by using $(k,\theta)$ and $(\alpha, \beta)$. $$\text{If } X \sim \text{Gamma}(k, \theta) , \,\,\,\, f(x) = \...


9

Assuming I take the mean of the posterior distribution rather than a random sample from it, is this what is commonly referred to as Rao-Blackwellization? I am not very familiar with stochastic volatility models, but I do know that in most settings, the reason we choose Gibbs or M-H algorithms to draw from the posterior, is because we don't know the ...


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

Computing a joint distribution from conditional distributions in general is very difficult. If the conditional distributions are chosen arbitrarily, a common joint distribution might not even exist. In this case, even showing that the conditional distributions are consistent is generally difficult. One result that might be used for deriving a joint ...


8

Those distributions you call "marginal" are not marginal. They are conditional distributions because you wrote $x \mid y$. The marginal distribution of $X$, for example, is necessarily independent of the value of $Y$. To see how the conditional distribution is gamma, all you have to do is write $$f_{X \mid Y}(x) = \frac{f_{X,Y}(x,y)}{f_Y(y)} \propto f_{X,...


8

As you show in the reproduction what is written in this book, the solution is incorrect for the simple reason that the quantity $(X^{T}X)^{-1}$ is a $p\times p$ matrix, not a scalar. Hence you cannot divide by $(X^{T}X)^{-1}$. (This is a terrible way of explaining this standard derivation!) What you can write instead is $$ \beta^{T}X^{T}X\beta - 2\beta^{T}X^...


8

What you are actually doing with the two-step process you've outlined is sampling from the joint distribution $p(x_{new}, \mu \thinspace | \thinspace x_1, \dots, x_n)$, then ignoring the sampled values of $\mu$. It's not altogether intuitive, but, by ignoring the sampled values of $\mu$, you are integrating over it. A simple example may make this clear. ...


8

The Gibbs sampler can then be used to improve efficiency of (say) samples from a marginal posterior, call it $\pi_2(\theta_2|y)$. Note \begin{eqnarray*} \pi_2(\theta_2|y)&=&\int \pi(\theta_1,\theta_2|y)d\theta_1\\ &=&\int \pi_{2|1}(\theta_2|\theta_1,y)\pi_1(\theta_1|y)d\theta_1\\ &=&E(\pi_{2|1}(\theta_2|\theta_1,y)) \end{eqnarray*} ...


8

The likelihood for $n$ iid observations looks like: $ f(x_1,...x_n|\lambda,\mu) \propto \frac{1}{\lambda^n} exp(-\frac{1}{\lambda}\sum_{i=1}^n|x_i-\mu|)$ Hence a conjugate prior for $\lambda$ with $\mu, x$ known must (thinking only about the algebra) look like: $ f(\lambda) \propto \frac{1}{\lambda^a} exp(-\frac{b}{\lambda})$ As suggested by marmle, this ...


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


7

There you go- Gibbs Sampler: The burning period is to reach some stationarity in the samples burning_period=5000 iterations=1000 y=matrix(nrow=(burning_period+iterations),ncol=3) a=matrix(nrow=iterations,ncol=3) y[1,1]=0.5 #Initial Sample y[1,2]=0.6 y[1,3]=0.2 for(i in 2:(burning_period+iterations)){ #I have put 3,3,4 as my theta's. You ...


7

Much faster than what? Univariate Gibbs sampling? The two stage Gibbs sampling is certainly the most studied type of Gibbs sampling starting with Tanner and Wong (1987, JASA). There is in particular a very achieved paper by Liu, Wong and Kong (1994, Biometrika), which shows that the correlation between the $X_t$'s (and the $Y_t$'s) is (a) positive and (b) ...


7

I would use a latent variable approach, since that's what x an y are. However, its not clear that all four parameters would be identifiable in this case. It would be helpful if you had some prior information for one or two of them. Here's an example: import pymc as pm # Priors mu_x = pm.Normal('mu_x', 0, 0.001, value=0) sigma_x = pm.Uniform('sigma_x', 0, ...


7

Since we are calculating the joint distribution, we'll assume that our initial sample is $ x = P(D=0,I=0,G=0,L=0,S=0) $ . To calculate the next sample, we'll need to sample each variable from the conditional distribution. $ P(D\mid G,I,S,L) $,from the conditional independencies in the Bayes net, simplifies to just sampling $ P(D)$. We sample and get the ...


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


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