Highest voted questions tagged gibbs - Cross Validated most recent 30 from stats.stackexchange.com 2019-07-17T11:12:55Z https://stats.stackexchange.com/feeds/tag?tagnames=gibbs&sort=votes http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://stats.stackexchange.com/q/10213 90 Can someone explain Gibbs sampling in very simple words? [duplicate] Thea https://stats.stackexchange.com/users/4429 2011-05-01T19:37:56Z 2019-06-20T01:59:21Z <div class="question-status question-originals-of-duplicate"> <p>This question already has an answer here:</p> <ul> <li> <a href="/questions/325696/explanation-regarding-gibbs-sampling" dir="ltr">Explanation regarding Gibbs Sampling</a> <span class="question-originals-answer-count"> 1 answer </span> </li> </ul> </div> <p>I'm doing some reading on topic modeling (with Latent Dirichlet Allocation) which makes use of Gibbs sampling. As a newbie in statistics―well, I know things like binomials, multinomials, priors, etc.―,I find it difficult to grasp how Gibbs sampling works. Can someone please explain it in simple English and/or using simple examples? (If you are not familiar with topic modeling, any examples will do.) </p> https://stats.stackexchange.com/q/9202 41 OpenBugs vs. JAGS DanB https://stats.stackexchange.com/users/3700 2011-04-05T15:42:20Z 2017-11-15T11:32:58Z <p>I am about to try out a BUGS style environment for estimating Bayesian models. Are there any important advantages to consider in choosing between OpenBugs or JAGS? Is one likely to replace the other in the foreseeable future? </p> <p>I will be using the chosen Gibbs Sampler with R. I don't have a specific application yet, but rather I am deciding which to intall and learn.</p> https://stats.stackexchange.com/q/185631 34 What is the difference between Metropolis Hastings, Gibbs, Importance, and Rejection sampling? user1398057 https://stats.stackexchange.com/users/53410 2015-12-08T06:54:50Z 2018-02-01T16:32:30Z <p>I have been trying to learn MCMC methods and have come across Metropolis Hastings, Gibbs, Importance, and Rejection sampling. While some of these differences are obvious, i.e., how Gibbs is a special case of Metropolis Hastings when we have the full conditionals, the others are less obvious, like when we want to use MH within a Gibbs sampler, etc. Does anyone have a simple way to see the bulk of the differences between each of these? Thanks!</p> https://stats.stackexchange.com/q/8485 29 A good Gibbs sampling tutorials and references fabrizioM https://stats.stackexchange.com/users/3788 2011-03-19T09:07:52Z 2012-03-19T00:50:49Z <p>I want to learn how Gibbs Sampling works and I am looking for a good basic to intermediate paper. I have a computer science background and basic statistic knowledge. </p> <p>Anyone has read good material around? where did you learn it?</p> <p>Thanks</p> https://stats.stackexchange.com/q/46561 21 What are some well known improvements over textbook MCMC algorithms that people use for bayesian inference? Rafael S. Calsaverini https://stats.stackexchange.com/users/972 2012-12-26T17:55:02Z 2014-09-30T12:52:49Z <p>When I'm coding a Monte Carlo simulation for some problem, and the model is simple enough, I use a very basic textbook Gibbs sampling. When it's not possible to use Gibbs sampling, I code the textbook Metropolis-Hastings I've learned years ago. The only thought I give to it is choosing the jumping distribution or its parameters.</p> <p>I know there are hundreds and hundreds of specialized methods that improve over those textbook options, but I usually never think about using/learning them. It usually feels like it's too much effort to improve a little bit what is already working very well. </p> <p>But recently I've been thinking if maybe there aren't new general methods that can improve over what I've been doing. It's been many decades since those methods were discovered. Maybe I'm <em>really</em> outdated!</p> <p><strong>Are there any well known alternatives to Metropolis-Hastings that are:</strong></p> <ul> <li><strong>reasonably easy to implement,</strong> </li> <li><strong>as universally appliable as MH,</strong> </li> <li><strong>and always improves over MH's results in some sense (computational performance, accuracy, etc...)?</strong></li> </ul> <p>I know about some very specialized improvements for very specialized models, but are there some general stuff everybody uses that I don't know? </p> https://stats.stackexchange.com/q/104815 20 Gibbs sampling versus general MH-MCMC Luca https://stats.stackexchange.com/users/36540 2014-06-26T06:08:18Z 2019-03-13T22:17:20Z <p>I have just been doing some reading on Gibbs sampling and Metropolis Hastings algorithm and have a couple of questions. </p> <p>As I understand it, in the case of Gibbs sampling, if we have a large multivariate problem, we sample from the conditional distribution i.e. sample one variable while keeping all others fixed whereas in MH, we sample from the full joint distribution.</p> <p>One thing the document said was that the proposed sample is <em>always</em> accepted in Gibbs Sampling i.e. the proposal acceptance rate is always 1. To me this seems like a big advantage as for large multivariate problems it seems that the rejection rate for MH algorithm becomes quite large. If that is indeed the case, what is the reason behind not using Gibbs Sampler all the time for generating the posterior distribution?</p> https://stats.stackexchange.com/q/244573 16 When would one use Gibbs sampling instead of Metropolis-Hastings? ShanZhengYang https://stats.stackexchange.com/users/80118 2016-11-06T22:32:03Z 2018-09-13T09:00:20Z <p>There are different kinds of MCMC algorithms: </p> <ul> <li>Metropolis-Hastings</li> <li>Gibbs</li> <li>Importance/rejection sampling (related). </li> </ul> <p>Why would one use Gibbs sampling instead of Metropolis-Hastings? I suspect there are cases when inference is more tractable with Gibbs sampling than with Metropolis-Hastings, but I am not clear on the specifics. </p> https://stats.stackexchange.com/q/118442 16 Does the Gibbs Sampling algorithm guarantee detailed balance? Ian https://stats.stackexchange.com/users/43122 2014-10-09T14:31:36Z 2016-01-11T20:56:15Z <p>I have it on supreme authority<sup><a href="http://en.wikipedia.org/wiki/Gibbs_sampling#Implementation" rel="noreferrer">1</a></sup> that Gibbs Sampling is a special case of the Metropolis-Hastings algorithm for Markov Chain Monte Carlo sampling. The MH algorithm always gives a transition probability with the detailed balance property; I expect Gibbs should too. So where in the following simple case have I gone wrong?</p> <p>For target distribution $\pi(x, y)$ on two discrete (for simplicity) variables, the full conditional distributions are: \begin{align} q_1 (x;y) &amp; =\frac{\pi (x,y)}{\sum_z \pi (z,y)} \\ q_2 (y;x) &amp; =\frac{\pi (x,y)}{\sum_z \pi (x,z)} \end{align}</p> <p>As I understand Gibbs Sampling, the transition probability can be written: $$Prob\{(y_1, y_2) \to (x_1, x_2)\} = q_1(x_1; y_2) q_2(x_2; x_1)$$</p> <p>The question is, does $$\pi(y_1,y_2) Prob\{(y_1, y_2) \to (x_1, x_2)\} \overset{?}{=} \pi(x_1,x_2) Prob\{(x_1, x_2) \to (y_1, y_2)\},$$ but the closest I can get is \begin{align} \pi(y_1,y_2) Prob\{(y_1, y_2) &amp; \to (x_1, x_2)\} \\ &amp; = \pi(y_1, y_2) q_2(x_2; x_1) q_1(x_1; y_2) \\ &amp; = \frac{\pi(x_1, x_2)}{\sum_z \pi(x_1,z)}\frac{\pi(x_1, y_2)}{\sum_z \pi(z, y_2)}\pi (y_1, y_2) \\ &amp; = \pi(x_1, x_2) q_2(y_2; x_1) q_1(y_1; y_2) \end{align} That's subtly different, and does not imply detailed balance. Thanks for any thoughts!</p> https://stats.stackexchange.com/q/4191 15 Where do the full conditionals come from in Gibbs sampling? cespinoza https://stats.stackexchange.com/users/1795 2010-11-04T04:35:38Z 2013-11-24T16:47:46Z <p>MCMC algorithms like Metropolis-Hastings and Gibbs sampling are ways of sampling from the joint posterior distributions. </p> <p>I think I understand and can implement metropolis-hasting pretty easily--you simply choose starting points somehow, and 'walk the parameter space' randomly, guided by the posterior density and proposal density. Gibbs sampling seems very similar but more efficient since it updates only one parameter at a time, while holding the others constant, effectively walking the space in an orthogonal fashion. </p> <p>In order to do this, you need the full conditional of each parameter in analytical from*. But where do these full conditionals come from? $$P(x_1 | x_2,\ \ldots,\ x_n) = \frac{P(x_1,\ \ldots,\ x_n)}{P(x_2,\ \ldots,\ x_n)}$$ To get the denominator you need to marginalize the joint over $x_1$. That seems like a whole lot of work to do analytically if there are many parameters, and might not be tractable if the joint distribution isn't very 'nice'. I realize that if you use conjugacy throughout the model, the full conditionals may be easy, but there's got to be a better way for more general situations.</p> <p>All the examples of Gibbs sampling I've seen online use toy examples (like sampling from a multivariate normal, where the conditionals are just normals themselves), and seem to dodge this issue. </p> <p>* Or do you need the full conditionals in analytical form at all? How do programs like winBUGS do it?</p> https://stats.stackexchange.com/q/348984 15 Stan $\hat{R}$ versus Gelman-Rubin $\hat{R}$ definition Greenparker https://stats.stackexchange.com/users/31978 2018-05-30T13:11:06Z 2018-06-11T13:24:17Z <p>I was going through the Stan documentation which can be downloaded from <a href="http://mc-stan.org/users/documentation/" rel="noreferrer">here</a>. I was particularly interested in their implementation of the Gelman-Rubin diagnostic. The original paper <a href="http://www.jstor.org/stable/2246093?seq=1#page_scan_tab_contents" rel="noreferrer">Gelman &amp; Rubin (1992)</a> define the the potential scale reduction factor (PSRF) as follows:</p> <p>Let $X_{i,1}, \dots , X_{i,N}$ be the $i$th Markov chain sampled, and let there be overall $M$ independent chains sampled. Let $\bar{X}_{i\cdot}$ be the mean from the $i$th chain, and $\bar{X}_{\cdot \cdot}$ be the overall mean. Define, $$W = \dfrac{1}{M} \sum_{m=1}^{M} {s^2_m},$$ where $$s^2_m = \dfrac{1}{N-1} \sum_{t=1}^{N} (\bar{X}_{m t} - \bar{X}_{m \cdot})^2\,.$$ And define $B$ $$B = \dfrac{N}{M-1} \sum_{m=1}^{M} (\bar{X}_{m \cdot} - \bar{X}_{\cdot \cdot})^2 \,.$$</p> <p>Define $$\hat{V} = \left(\dfrac{N-1}{N} \right)W + \left( \dfrac{M+1}{MN} \right)B\,.$$ The PSRF is estimate with $\sqrt{\hat{R}}$ where $$\hat{R} = \dfrac{\hat{V}}{W} \cdot \dfrac{df+3}{df+1}\,,$$ where $df = 2\hat{V}/Var(\hat{V})$.</p> <p>The Stan documentation on page 349 ignores the term with $df$ and also removes the $(M+1)/M$ multiplicative term. This is their formula,</p> <blockquote> <p>The variance estimator is $$\widehat{\text{var}}^{+}(\theta \, | \, y) = \frac{N-1}{N} W + \frac{1}{N} B\,.$$ Finally, the potential scale reduction statistic is defined by $$\hat{R} = \sqrt{\frac{\widehat{\text{var}}^{+}(\theta \, | \, y) }{W}}\,.$$</p> </blockquote> <p>From what I could see, they don't provide a reference for this change of formula, and neither do they discuss it. Usually $M$ is not too large, and can often be as low so as $2$, so $(M+1)/M$ should not be ignored, even if the $df$ term can be approximated with 1. </p> <p>So where does this formula come from?</p> <hr> <p><strong>EDIT:</strong> I have found a partial answer to the question "<em>where does this formula come from?</em>", in that the <a href="http://www.stat.columbia.edu/~gelman/book/" rel="noreferrer">Bayesian Data Analysis book by Gelman, Carlin, Stern, and Rubin</a> (Second edition) has exactly the same formula. However, the book does not explain how/why it is justifiable to ignore those terms?</p> https://stats.stackexchange.com/q/129109 12 Marginal Likelihood from the Gibbs Output Zen https://stats.stackexchange.com/users/9394 2014-12-15T01:18:41Z 2015-01-05T14:18:34Z <p>I'm reproducing from scratch the results in Section 4.2.1 of</p> <p><strong><a href="http://apps.olin.wustl.edu/faculty/chib/papers/chib95.pdf" rel="noreferrer">Marginal Likelihood from the Gibbs Output</a></strong></p> <p><em>Siddhartha Chib</em></p> <p>Journal of the American Statistical Association, Vol. 90, No. 432. (Dec., 1995), pp. 1313-1321.</p> <p>It's a mixture of normals model with known number $k\geq 1$ of components. $$f(x\mid w,\mu,\sigma^2) =\prod_{i=1}^n\sum_{j=1}^k \mathrm{N}(x_i\mid\mu_j,\sigma_j^2) \, . \qquad (*)$$</p> <p>The Gibbs sampler for this model is implemented using the data augmentation technique of Tanner and Wong. A set of allocation variables $z=(z_1,\dots,z_n)$ assuming the values $1,\dots,k$ is introduced, and we specify that $\Pr(z_i=j\mid w)=w_j$ and $f(x_i\mid z,\mu,\sigma^2)=\mathrm{N}(x_i\mid\mu_{z_i},\sigma^2_{z_i})$. It follows that integration over the $z_i$'s gives the original likelihood $(*)$.</p> <p>The dataset is formed by velocities of $82$ galaxies from the Corona Borealis constellation.</p> <pre><code>set.seed(1701) x &lt;- c( 9.172, 9.350, 9.483, 9.558, 9.775, 10.227, 10.406, 16.084, 16.170, 18.419, 18.552, 18.600, 18.927, 19.052, 19.070, 19.330, 19.343, 19.349, 19.440, 19.473, 19.529, 19.541, 19.547, 19.663, 19.846, 19.856, 19.863, 19.914, 19.918, 19.973, 19.989, 20.166, 20.175, 20.179, 20.196, 20.215, 20.221, 20.415, 20.629, 20.795, 20.821, 20.846, 20.875, 20.986, 21.137, 21.492, 21.701, 21.814, 21.921, 21.960, 22.185, 22.209, 22.242, 22.249, 22.314, 22.374, 22.495, 22.746, 22.747, 22.888, 22.914, 23.206, 23.241, 23.263, 23.484, 23.538, 23.542, 23.666, 23.706, 23.711, 24.129, 24.285, 24.289, 24.366, 24.717, 24.990, 25.633, 26.960, 26.995, 32.065, 32.789, 34.279 ) nn &lt;- length(x) </code></pre> <p>We assume that $w$, the $\mu_j$'s, and the $\sigma^2_j$'s are independent <em>a priori</em> with $$(w_1,\dots,w_k) \sim \mathrm{Dir}(a_1,\dots,a_k) \, , \quad \mu_j \sim \mathrm{N}(\mu_0,\sigma_0^2) \, , \quad \sigma^2_j\sim\mathrm{IG}\!\left(\frac{\nu_0}{2},\frac{\delta_0}{2}\right) \, .$$</p> <pre><code>k &lt;- 3 mu0 &lt;- 20 va0 &lt;- 100 nu0 &lt;- 6 de0 &lt;- 40 a &lt;- rep(1, k) </code></pre> <p>Using Bayes' Theorem, the full conditionals are \begin{align*} w \mid \mu,\sigma^2,z,x &amp;\sim \mathrm{Dir}(a_1+n_1,\dots,a_k+n_k) \\ \mu_j \mid w, \sigma^2,z,x &amp;\sim \mathrm{N}\!\left( \frac{n_j m_j\sigma_0^2+\mu_0\sigma_j^2}{n_j\sigma^2_0+\sigma^2_j}, \frac{\sigma^2_0\sigma^2_j}{n_j\sigma^2_0+\sigma^2_j}\right) \\ \sigma_j^2 \mid w,\mu,z,x &amp;\sim \mathrm{IG}\!\left( \frac{\nu_0+n_j}{2},\frac{\delta_0+\delta_j}{2}\right) \\ \Pr(z_i=j\mid w,\mu,\sigma^2,x) &amp;\propto w_j \times \frac{1}{\sigma_j}e^{-(x_i-\mu_j)^2/2\sigma_j^2} \end{align*} in which $$n_j = |L_j| \, , \qquad m_j = \begin{cases}\frac{1}{n_j}\sum_{i\in L_j} x_i &amp;\;\mathrm{if}\; n_j&gt;0 \\ 0 &amp;\;\mathrm{otherwise.} \end{cases}\, , \qquad \delta_j = \sum_{i\in L_j} (x_i-\mu_j)^2 \, ,$$ with $L_j=\{i\in\{1,\dots,n\}:z_i=j\}$.</p> <p>The goal is to compute an estimate for the marginal likelihood of the model. Chib's method begins with a first run of the Gibbs sampler using the full conditionals.</p> <pre><code>burn_in &lt;- 1000 run &lt;- 15000 cat("First Gibbs run (full):\n") N &lt;- burn_in + run w &lt;- matrix(1, nrow = N, ncol = k) mu &lt;- matrix(0, nrow = N, ncol = k) va &lt;- matrix(1, nrow = N, ncol = k) z &lt;- matrix(1, nrow = N, ncol = nn) n &lt;- integer(k) m &lt;- numeric(k) de &lt;- numeric(k) rdirichlet &lt;- function(a) { y &lt;- rgamma(length(a), a, 1); y / sum(y) } pb &lt;- txtProgressBar(min = 2, max = N, style = 3) z[1,] &lt;- sample.int(k, size = nn, replace = TRUE) for (t in 2:N) { n &lt;- tabulate(z[t-1,], nbins = k) w[t,] &lt;- rdirichlet(a + n) m &lt;- sapply(1:k, function(j) sum(x[z[t-1,]==j])) m[n &gt; 0] &lt;- m[n &gt; 0] / n[n &gt; 0] mu[t,] &lt;- rnorm(k, mean = (n*m*va0+mu0*va[t-1,])/(n*va0+va[t-1,]), sd = sqrt(va0*va[t-1,]/(n*va0+va[t-1,]))) de &lt;- sapply(1:k, function(j) sum((x[z[t-1,]==j] - mu[t,j])^2)) va[t,] &lt;- 1 / rgamma(k, shape = (nu0+n)/2, rate = (de0+de)/2) z[t,] &lt;- sapply(1:nn, function(i) sample.int(k, size = 1, prob = exp(log(w[t,]) + dnorm(x[i], mean = mu[t,], sd = sqrt(va[t,]), log = TRUE)))) setTxtProgressBar(pb, t) } close(pb) </code></pre> <p>From this first run we get an approximate point $(w^*,\mu^*,\sigma^{2*})$ of maximum likelihood. Since the likelihood is actually unbounded, what this procedure probably gives is an approximate local MAP.</p> <pre><code>w &lt;- w[(burn_in+1):N,] mu &lt;- mu[(burn_in+1):N,] va &lt;- va[(burn_in+1):N,] z &lt;- z[(burn_in+1):N,] N &lt;- N - burn_in log_L &lt;- function(x, w, mu, va) sum(log(sapply(1:nn, function(i) sum(exp(log(w) + dnorm(x[i], mean = mu, sd = sqrt(va), log = TRUE)))))) ts &lt;- which.max(sapply(1:N, function(t) log_L(x, w[t,], mu[t,], va[t,]))) ws &lt;- w[ts,] mus &lt;- mu[ts,] vas &lt;- va[ts,] </code></pre> <p>Chib's log-estimate of the marginal likelihood is \begin{align} \log \widehat{f(x)} &amp;= \log L_x(w^*,\mu^*,\sigma^{2*}) + \log \pi(w^*,\mu^*,\sigma^{2*}) \\ &amp;- \log \pi(\mu^*\mid x) - \log \pi(\sigma^{2*}\mid \mu^*,x) - \log \pi(w^*\mid \mu^*,\sigma^{2*},x) \, . \end{align}</p> <p>We already have the first two terms.</p> <pre><code>log_prior &lt;- function(w, mu, va) { lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w)) + sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE)) + sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va)) } chib &lt;- log_L(x, ws, mus, vas) + log_prior(ws, mus, vas) </code></pre> <p>The Rao-Blackwellized estimate of $\pi(\mu^*\mid x)$ is $$\pi(\mu^*\mid x) = \int \prod_{j=1}^k \mathrm{N}\!\left(\mu_j^* \;\Bigg|\; \frac{n_j m_j\sigma_0^2+\mu_0\sigma_j^2}{n_j\sigma^2_0+\sigma^2_j}, \frac{\sigma^2_0\sigma^2_j}{n_j\sigma^2_0+\sigma^2_j}\right)\,p(\sigma^{2},z\mid x)\,d\sigma^2\,dz \, ,$$ and is readily obtained from the first Gibbs run.</p> <pre><code>pi.mu_va.z.x &lt;- function(mu, va, z) { n &lt;- tabulate(z, nbins = k) m &lt;- sapply(1:k, function(j) sum(x[z==j])) m[n &gt; 0] &lt;- m[n &gt; 0] / n[n &gt; 0] exp(sum(dnorm(mu, mean = (n*m*va0+mu0*va)/(n*va0+va), sd = sqrt(va0*va/(n*va0+va)), log = TRUE))) } chib &lt;- chib - log(mean(sapply(1:N, function(t) pi.mu_va.z.x(mus, va[t,], z[t,])))) </code></pre> <p>The Rao-Blackwellized estimate of $\pi(\sigma^{2*}\mid \mu^*,x)$ is $$\pi(\sigma^{2*}\mid \mu^*,x) = \int \prod_{j=1}^k \mathrm{IG}\!\left( \sigma^{2*}_j \;\Bigg|\; \frac{\nu_0+n_j}{2},\frac{\delta_0+\delta_j}{2}\right) \, p(z\mid\mu^*,x)\,dz \, ,$$ and is computed from a second reduced Gibbs run in which the $\mu_j$'s are not updated, but made equal to $\mu^*_j$ at each iteration step.</p> <pre><code>cat("Second Gibbs run (reduced):\n") N &lt;- burn_in + run w &lt;- matrix(1, nrow = N, ncol = k) va &lt;- matrix(1, nrow = N, ncol = k) z &lt;- matrix(1, nrow = N, ncol = nn) pb &lt;- txtProgressBar(min = 2, max = N, style = 3) z[1,] &lt;- sample.int(k, size = nn, replace = TRUE) for (t in 2:N) { n &lt;- tabulate(z[t-1,], nbins = k) w[t,] &lt;- rdirichlet(a + n) de &lt;- sapply(1:k, function(j) sum((x[z[t-1,]==j] - mus[j])^2)) va[t,] &lt;- 1 / rgamma(k, shape = (nu0+n)/2, rate = (de0+de)/2) z[t,] &lt;- sapply(1:nn, function(i) sample.int(k, size = 1, prob = exp(log(w[t,]) + dnorm(x[i], mean = mus, sd = sqrt(va[t,]), log = TRUE)))) setTxtProgressBar(pb, t) } close(pb) w &lt;- w[(burn_in+1):N,] va &lt;- va[(burn_in+1):N,] z &lt;- z[(burn_in+1):N,] N &lt;- N - burn_in pi.va_mu.z.x &lt;- function(va, mu, z) { n &lt;- tabulate(z, nbins = k) de &lt;- sapply(1:k, function(j) sum((x[z==j] - mu[j])^2)) exp(sum(((nu0+n)/2)*log((de0+de)/2) - lgamma((nu0+n)/2) - ((nu0+n)/2+1)*log(va) - (de0+de)/(2*va))) } chib &lt;- chib - log(mean(sapply(1:N, function(t) pi.va_mu.z.x(vas, mus, z[t,])))) </code></pre> <p>In the same way, the Rao-Blackwellized estimate of $\pi(w^*\mid \mu^*,\sigma^{2*},x)$ is $$\pi(w^*\mid \mu^*,\sigma^{2*},x) = \int \mathrm{Dir}(w^* \mid a_1+n_1,\dots,a_k+n_k) \, p(z\mid\mu^*,\sigma^{2*},x)\,dz \, ,$$ and is computed from a third reduced Gibbs run in which the $\mu_j$'s and the $\sigma^2_j$'s are not updated, but made equal to $\mu^*_j$ and $\sigma^{2*}_j$ respectively at each iteration step.</p> <pre><code>cat("Third Gibbs run (reduced):\n") N &lt;- burn_in + run w &lt;- matrix(1, nrow = N, ncol = k) z &lt;- matrix(1, nrow = N, ncol = nn) pb &lt;- txtProgressBar(min = 2, max = N, style = 3) z[1,] &lt;- sample.int(k, size = nn, replace = TRUE) for (t in 2:N) { n &lt;- tabulate(z[t-1,], nbins = k) w[t,] &lt;- rdirichlet(a + n) z[t,] &lt;- sapply(1:nn, function(i) sample.int(k, size = 1, prob = exp(log(w[t,]) + dnorm(x[i], mean = mus, sd = sqrt(vas), log = TRUE)))) setTxtProgressBar(pb, t) } close(pb) w &lt;- w[(burn_in+1):N,] z &lt;- z[(burn_in+1):N,] N &lt;- N - burn_in pi.w_z.x &lt;- function(w, z) { n &lt;- tabulate(z, nbins = k) exp(lgamma(sum(a+n)) - sum(lgamma(a+n)) + sum((a+n-1)*log(w))) } chib &lt;- chib - log(mean(sapply(1:N, function(t) pi.w_z.x(ws, z[t,])))) </code></pre> <p>After all this, we get a log-estimate $-217.9199$ which is bigger than the one reported by Chib: $-224.138$ with Monte Carlo error $.086$. </p> <p>To check if I somehow messed up the Gibbs samplers, I reimplemented the whole thing using RJAGS. The following code gives the same results.</p> <pre><code>x &lt;- c( 9.172, 9.350, 9.483, 9.558, 9.775, 10.227, 10.406, 16.084, 16.170, 18.419, 18.552, 18.600, 18.927, 19.052, 19.070, 19.330, 19.343, 19.349, 19.440, 19.473, 19.529, 19.541, 19.547, 19.663, 19.846, 19.856, 19.863, 19.914, 19.918, 19.973, 19.989, 20.166, 20.175, 20.179, 20.196, 20.215, 20.221, 20.415, 20.629, 20.795, 20.821, 20.846, 20.875, 20.986, 21.137, 21.492, 21.701, 21.814, 21.921, 21.960, 22.185, 22.209, 22.242, 22.249, 22.314, 22.374, 22.495, 22.746, 22.747, 22.888, 22.914, 23.206, 23.241, 23.263, 23.484, 23.538, 23.542, 23.666, 23.706, 23.711, 24.129, 24.285, 24.289, 24.366, 24.717, 24.990, 25.633, 26.960, 26.995, 32.065, 32.789, 34.279 ) library(rjags) nn &lt;- length(x) k &lt;- 3 mu0 &lt;- 20 va0 &lt;- 100 nu0 &lt;- 6 de0 &lt;- 40 a &lt;- rep(1, k) burn_in &lt;- 10^3 N &lt;- 10^4 full &lt;- " model { for (i in 1:n) { x[i] ~ dnorm(mu[z[i]], tau[z[i]]) z[i] ~ dcat(w[]) } for (i in 1:k) { mu[i] ~ dnorm(mu0, 1/va0) tau[i] ~ dgamma(nu0/2, de0/2) va[i] &lt;- 1/tau[i] } w ~ ddirich(a) } " data &lt;- list(x = x, n = nn, k = k, mu0 = mu0, va0 = va0, nu0 = nu0, de0 = de0, a = a) model &lt;- jags.model(textConnection(full), data = data, n.chains = 1, n.adapt = 100) update(model, n.iter = burn_in) samples &lt;- jags.samples(model, c("mu", "va", "w", "z"), n.iter = N) mu &lt;- matrix(samples$mu, nrow = N, byrow = TRUE) va &lt;- matrix(samples$va, nrow = N, byrow = TRUE) w &lt;- matrix(samples$w, nrow = N, byrow = TRUE) z &lt;- matrix(samples$z, nrow = N, byrow = TRUE) log_L &lt;- function(x, w, mu, va) sum(log(sapply(1:nn, function(i) sum(exp(log(w) + dnorm(x[i], mean = mu, sd = sqrt(va), log = TRUE)))))) ts &lt;- which.max(sapply(1:N, function(t) log_L(x, w[t,], mu[t,], va[t,]))) ws &lt;- w[ts,] mus &lt;- mu[ts,] vas &lt;- va[ts,] log_prior &lt;- function(w, mu, va) { lgamma(sum(a)) - sum(lgamma(a)) + sum((a-1)*log(w)) + sum(dnorm(mu, mean = mu0, sd = sqrt(va0), log = TRUE)) + sum((nu0/2)*log(de0/2) - lgamma(nu0/2) - (nu0/2+1)*log(va) - de0/(2*va)) } chib &lt;- log_L(x, ws, mus, vas) + log_prior(ws, mus, vas) cat("log-likelihood + log-prior =", chib, "\n") pi.mu_va.z.x &lt;- function(mu, va, z, x) { n &lt;- sapply(1:k, function(j) sum(z==j)) m &lt;- sapply(1:k, function(j) sum(x[z==j])) m[n &gt; 0] &lt;- m[n &gt; 0] / n[n &gt; 0] exp(sum(dnorm(mu, mean = (n*m*va0+mu0*va)/(n*va0+va), sd = sqrt(va0*va/(n*va0+va)), log = TRUE))) } chib &lt;- chib - log(mean(sapply(1:N, function(t) pi.mu_va.z.x(mus, va[t,], z[t,], x)))) cat("log-likelihood + log-prior - log-pi.mu_ =", chib, "\n") fixed.mu &lt;- " model { for (i in 1:n) { x[i] ~ dnorm(mus[z[i]], tau[z[i]]) z[i] ~ dcat(w[]) } for (i in 1:k) { tau[i] ~ dgamma(nu0/2, de0/2) va[i] &lt;- 1/tau[i] } w ~ ddirich(a) } " data &lt;- list(x = x, n = nn, k = k, nu0 = nu0, de0 = de0, a = a, mus = mus) model &lt;- jags.model(textConnection(fixed.mu), data = data, n.chains = 1, n.adapt = 100) update(model, n.iter = burn_in) samples &lt;- jags.samples(model, c("va", "w", "z"), n.iter = N) va &lt;- matrix(samples$va, nrow = N, byrow = TRUE) w &lt;- matrix(samples$w, nrow = N, byrow = TRUE) z &lt;- matrix(samples$z, nrow = N, byrow = TRUE) pi.va_mu.z.x &lt;- function(va, mu, z, x) { n &lt;- sapply(1:k, function(j) sum(z==j)) de &lt;- sapply(1:k, function(j) sum((x[z==j] - mu[j])^2)) exp(sum(((nu0+n)/2)*log((de0+de)/2) - lgamma((nu0+n)/2) - ((nu0+n)/2+1)*log(va) - (de0+de)/(2*va))) } chib &lt;- chib - log(mean(sapply(1:N, function(t) pi.va_mu.z.x(vas, mus, z[t,], x)))) cat("log-likelihood + log-prior - log-pi.mu_ - log-pi.va_ =", chib, "\n") fixed.mu.and.va &lt;- " model { for (i in 1:n) { x[i] ~ dnorm(mus[z[i]], 1/vas[z[i]]) z[i] ~ dcat(w[]) } w ~ ddirich(a) } " data &lt;- list(x = x, n = nn, a = a, mus = mus, vas = vas) model &lt;- jags.model(textConnection(fixed.mu.and.va), data = data, n.chains = 1, n.adapt = 100) update(model, n.iter = burn_in) samples &lt;- jags.samples(model, c("w", "z"), n.iter = N) w &lt;- matrix(samples$w, nrow = N, byrow = TRUE) z &lt;- matrix(samplesz, nrow = N, byrow = TRUE) pi.w_z.x &lt;- function(w, z, x) { n &lt;- sapply(1:k, function(j) sum(z==j)) exp(lgamma(sum(a)+nn) - sum(lgamma(a+n)) + sum((a+n-1)*log(w))) } chib &lt;- chib - log(mean(sapply(1:N, function(t) pi.w_z.x(ws, z[t,], x)))) cat("log-likelihood + log-prior - log-pi.mu_ - log-pi.va_ - log-pi.w_ =", chib, "\n") </code></pre> <p>My question is if in the above description there are any misunderstandings of Chib's method or any mistakes in its implementation.</p> https://stats.stackexchange.com/q/296415 12 Why does the redundant mean parameterization speed up Gibbs MCMC? Heisenberg https://stats.stackexchange.com/users/20148 2017-08-06T00:38:41Z 2018-02-06T15:37:33Z <p>In Gelman &amp; Hill (2007)'s book (Data Analysis Using Regression and Multilevel/Hierarchical Models), the authors claim that including redundant mean parameters can help speed up MCMC.</p> <p>The given example is a non-nested model of "flight simulator" (Eq 13.9):</p> <p>\begin{align} y_i &amp;\sim N(\mu + \gamma_{j[i]} + \delta_{k[i]}, \sigma^2_y) \\ \gamma_j &amp;\sim N(0, \sigma^2_\gamma) \\ \delta_k &amp;\sim N(0, \sigma^2_\delta) \end{align}</p> <p>They recommend a reparameterization, adding the mean parameters\mu_\gamma$and$\mu_\deltaas follows:</p> <p>\begin{align} \gamma_j \sim N(\mu_\gamma, \sigma^2_\gamma) \\ \delta_k \sim N(\mu_\delta, \sigma^2_\delta) \end{align}</p> <p>The only justification offered is that (p. 420):</p> <blockquote> <p>It is possible for the simulations to get stuck in a configuration where the entire vector\gamma$(or$\delta$) is far from zero (even though they are assigned a distribution with mean 0). Ultimately, the simulations will converge to the correct distribution, but we do not want to have to wait.</p> </blockquote> <p>How do the redundant mean parameters help with this problem? </p> <p>It seems to me that the non-nested model is slow mainly because of$\gamma$and$\delta$are negatively correlated. (Indeed, if one goes up, the other has to go down, given that their sum is "fixed" by the data). Do the redundant mean parameters help with reducing the correlation between$\gamma$and$\delta$, or something else entirely?</p> https://stats.stackexchange.com/q/379508 11 Is Gibbs sampling an MCMC method? Gabriel https://stats.stackexchange.com/users/10416 2018-11-29T21:52:40Z 2018-12-03T01:43:50Z <p>As far as I understand it, it is (at least, that is how <a href="https://en.wikipedia.org/wiki/Gibbs_sampling" rel="noreferrer">Wikipedia defines it</a>). But I've found this statement by Efron* (emphasis added):</p> <blockquote> <p>Markov chain Monte Carlo (MCMC) is the great success story of modern-day Bayesian statistics. <strong>MCMC, and its sister method “Gibbs sampling,”</strong> permit the numerical calculation of posterior distributions in situations far too complicated for analytic expression.</p> </blockquote> <p>and now I'm confused. Is this just a minor difference in terminology, or is Gibbs sampling something other than MCMC?</p> <p>[*]: <a href="https://doi.org/10.1080/10543406.2011.607736" rel="noreferrer">Efron 2011, "The Bootstrap and Markov-Chain Monte Carlo"</a></p> https://stats.stackexchange.com/q/94792 11 Gibbs sampling for Ising model Collin https://stats.stackexchange.com/users/31632 2014-04-23T01:16:38Z 2018-04-12T00:53:50Z <p>Homework question:</p> <p>Consider the 1-d Ising model. </p> <p>Let$x = (x_1,...x_d)$.$x_i$is either -1 or +1</p> <p>$\pi(x) \propto e^{\sum_{i=1}^{39}x_ix_{i+1}}$</p> <p>Design a gibbs sampling algorithm to generate samples approximately from target distribution$\pi(x)$. </p> <p><strong>My attempt:</strong></p> <p>Randomly choose values (either -1 or 1) to fill vector$x = (x_1,...x_{40})$. So maybe$x = (-1, -1, 1, 1, 1, -1, 1, 1,...,1)$. So this is$x^0$.</p> <p>So now we need to move on and do the first iteration. We have to draw the 40 different x's for$x^1$separately. So...</p> <p>Draw$x_1^1$from$\pi(x_1 | x_2^0,...,x_{40}^0)$</p> <p>Draw$x_2^1$from$\pi(x_2 | x_1^1, x_3^0,...,x_{40}^0)$</p> <p>Draw$x_3^1$from$\pi(x_3 | x_1^1, x_2^1, x_4^0,...,x_{40}^0)$</p> <p>Etc..</p> <p>So the part that's tripping me up is how do we actually draw from the conditional distribution. How does$\pi(x) \propto e^{\sum_{i=1}^{39}x_ix_{i+1}}$come into play? Maybe an example of one draw would clear things up.</p> https://stats.stackexchange.com/q/129956 11 How to derive Gibbs sampling? Anne van Rossum https://stats.stackexchange.com/users/27826 2014-12-22T00:16:12Z 2014-12-22T23:44:34Z <p>I'm actually hesitating to ask this, because I'm afraid I will be referred to other questions or Wikipedia on Gibbs sampling, but I don't have the feeling that they describe what's at hand.</p> <p>Given a conditional probability$p(x|y)$: $$\begin{array}{c|c|c} p(x|y) &amp; y = y_0 &amp; y = y_1 \\ \hline x = x_0 &amp; \tfrac{1}{4} &amp; \tfrac{2}{6} \\ \hline x = x_1 &amp; \tfrac{3}{4} &amp; \tfrac{4}{6} \\ \end{array}$$</p> <p>And a conditional probability$p(y|x)$: $$\begin{array}{c|c|c} p(y|x) &amp; y = y_0 &amp; y = y_1 \\ \hline x = x_0 &amp; \tfrac{1}{3} &amp; \tfrac{2}{3} \\ \hline x = x_1 &amp; \tfrac{3}{7} &amp; \tfrac{4}{7} \\ \end{array}$$</p> <p>We can uniquely come up with the joint probability$f_{unique}=p(x,y)$:</p> <p>$$\begin{array}{c|c|c|c} p(x,y) &amp; y = y_0 &amp; y = y_1 &amp; p(x) \\ \hline x = x_0 &amp; a_0 &amp; a_1 &amp; c_0 \\ \hline x = x_1 &amp; a_2 &amp; a_3 &amp; c_1 \\ \hline p(y) &amp; b_0 &amp; b_1 &amp; \\ \end{array}$$</p> <p>Because, although we have$8$unknowns, we have more ($4*2+3$) linear equations:</p> <p>$ a_0+a_1+a_2+a_3=1 \\ b_0+b_1 = 1 \\ c_0+c_1 = 1 $</p> <p>As well as:</p> <p>$ \tfrac{1}{4} b_0 = a_0 \\ \tfrac{3}{4} b_0 = a_2 \\ \tfrac{2}{6} (1-b_0) = a_1 \\ \tfrac{4}{6} (1-b_0) = a_3 \\ \tfrac{1}{3} c_0 = a_0 \\ \tfrac{2}{3} c_0 = a_1 \\ \tfrac{3}{7} (1-c_0) = a_2 \\ \tfrac{4}{7} (1-c_0) = a_3 $</p> <p>It's quickly solved by$c_0=\tfrac{3}{4}b_0$,$\tfrac{2}{3}c_0=a_1$. Namely by equating$\tfrac{2}{4}b_0=a_1$with$\tfrac{2}{6}(1-b_0)=a_1$. This gives$b_0=\tfrac{2}{5}$and the rest follows.</p> <p>$$\begin{array}{c|c|c|c} p(x,y) &amp; y = y_0 &amp; y = y_1 &amp; p(x) \\ \hline x = x_0 &amp; \tfrac{1}{10} &amp; \tfrac{2}{10} &amp; \tfrac{3}{10} \\ \hline x = x_1 &amp; \tfrac{3}{10} &amp; \tfrac{4}{10} &amp; \tfrac{7}{10} \\ \hline p(y) &amp; \tfrac{4}{10} &amp; \tfrac{6}{10} &amp; \\ \end{array}$$</p> <p>So, now we go to the continuous case. It is imaginable to go to intervals and keep the above structure in-tact (with more equations than unknowns). However, what happens when we go to (point) instances of random variables? How does sampling</p> <p>$$x_a \sim p(x|y=y_b) \\ y_b \sim p(y|x=x_a)$$</p> <p>iteratively, lead to$p(x,y)$? Equivalent to the constraint$a_0 + a_1 + a_2 + a_3=1$, how does it ensure$\int_X \int_Y p(x,y) dy dx = 1$for example? Likewise with$\int_Y p(y|x)dy=1$. Can we write down the constraints and derive Gibbs sampling from first principles?</p> <p>So, I'm not interested in how to perform Gibbs sampling, which is simple, but I'm interested in how to derive it, and preferably how to prove that it works (probably under certain conditions).</p> https://stats.stackexchange.com/q/68171 11 How do programs like BUGS/JAGS automatically determine conditional distributions for Gibbs sampling? user4733 https://stats.stackexchange.com/users/4733 2013-08-23T13:23:34Z 2014-09-26T02:04:47Z <p>Seems like full conditionals are often quite difficult to derive, yet programs like JAGS and BUGS derive them automatically. Can someone explain how they algorithmically generate full conditionals for any arbitrary model specification?</p> https://stats.stackexchange.com/q/189454 11 How to test if a cross-covariance matrix is non-zero? TomHall https://stats.stackexchange.com/users/12035 2016-01-06T00:20:18Z 2018-02-27T15:23:13Z <p><strong>The background of my study</strong>: </p> <p>In a Gibbs sampling where we sample$X$(the variable of interests) and$Y$from$P(X|Y)$and$P(Y|X)$respectively, where$X$and$Y$are$k$-dimensional random vectors. We know that the process is usually split into two stages:</p> <ol> <li>Burn-in Period, where we discard all the samples. Denote the samples as$X_1\sim X_t$and$Y_1\sim Y_t$.</li> <li>"After-Burn-in" Period, where we average the samples$\bar{X} = \frac{1}{k}\sum_{i=1}^k X_{t+i}$as our final desired result.</li> </ol> <p>However, the samples in the "after-burn-in" sequence$X_{t+1}\sim X_{t+k}$are not independently distributed. Therefore if I want to inspect the variance of the final result, it becomes</p> <p>$$\operatorname{Var}[\bar{X}] = \operatorname{Var}\left[\sum_{i=1}^k X_{t+i}\right] = \frac{1}{k^2}\left(\sum_{i=1}^k\operatorname{Var}[X_{t+i}] + \sum_{i=1}^{k-1} \sum_{j=i+1}^k \operatorname{Cov}[X_{t+i},X_{t+j}]\right)$$</p> <p>Here the term$\operatorname{Cov}[X_{t+i},X_{t+j}]$is a$k\times k$cross-covariance matrix applies to any$(i,j)$with$i&lt;j$. </p> <p>For example, I have </p> <p>$$X_{t+1} = (1,2,1)'\\ X_{t+2} = (1,0,2)'\\ X_{t+3} = (1,0,0)'\\ X_{t+4} = (5,0,-1)'$$</p> <p>then I could estimate the covariance matrix$\operatorname{Cov}[X_{t+i}, X_{t+i+1}]$with </p> <p>$$\frac{1}{3}\sum_{i=1}^3 (X_{t+i}-\mu_{t+i})(X_{t+i+1}-\mu_{t+i+1})'$$</p> <p>Now I am interested in if the resulting estimation is significantly non-zero so that I need to include it into my variance estimation of$\operatorname{Var}[\bar{X}]$.</p> <p><strong>So here comes my questions</strong>:</p> <ol> <li>We sample$X_{t+i}$from$P(X_{t+i}|Y_{t+i})$. Since$Y_{t+i}$is changing, I think$X_{t+i}$and$X_{t+i+1}$are not from the same distribution, so$\operatorname{Cov}[X_{t+i},X_{t+j}]$is not the same as$\operatorname{Cov}[X_{t+i},X_{t+i}]$. Is this statement correct?</li> <li>Suppose I have enough data to estimate$\operatorname{Cov}[X_{t+i},X_{t+i+1}]$(neighboring samples in the sequence), is there any way to test if the covariance matrix is significantly a non-zero matrix? Broadly speaking, I am interested in an indicator which guides me to some <em>meaningful</em> cross-covariance matrices that should be included in my final variance estimation.</li> </ol> https://stats.stackexchange.com/q/116179 10 Bayesian modeling using multivariate normal with covariate Robert Smith https://stats.stackexchange.com/users/2676 2014-09-18T21:14:29Z 2014-12-29T17:47:10Z <p>Suppose you have an explanatory variable${\bf{X}} = \left(X(s_{1}),\ldots,X(s_{n})\right)$where$s$represents a given coordinate. You also have a response variable${\bf{Y}} = \left(Y(s_{1}),\ldots,Y(s_{n})\right)$. Now, we can combine both variables as: </p> <p>$${\bf{W}}({\bf{s}}) = \left( \begin{array}{ccc}X(s) \\ Y(s) \end{array} \right) \sim N(\boldsymbol{\mu}(s), T)$$</p> <p>In this case, we simply choose$\boldsymbol{\mu}(s) = \left( \mu_{1} \; \; \mu_{2}\right)^{T}$and$T$is a covariance matrix that describes the relation between$X$and$Y$. This only describes the value of$X$and$Y$at$s$. Since we have more points from other locations for$X$and$Y$, we can describe more values of${\bf{W}}(s)$in the following way:</p> <p>$$\left( \begin{array}{ccc} {\bf{X}} \\ {\bf{Y}} \end{array}\right) = N\left(\left(\begin{array}{ccc}\mu_{1}\boldsymbol{1}\\ \mu_{2}\boldsymbol{1}\end{array}\right), T\otimes H(\phi)\right)$$</p> <p>You will notice that we rearranged the components of$\bf{X}$and$\bf{Y}$to get all$X(s_i)$in a column and after that, concatenate all$Y(s_i)$together. Each component$H(\phi)_{ij}$is a correlation function$\rho(s_i, s_j)$and$T$is as above. The reason we have the covariance$T\otimes H(\phi)$is because we assume it is possible to separate the covariance matrix as$C(s, s')=\rho(s, s') T$.</p> <p>Question 1: When I calculate the conditional${\bf{Y}}\mid{\bf{X}}$, what I'm actually doing is generating a set of values of$\bf{Y}$based on$\bf{X}$, correct? I already have$\bf{Y}$so I would be more interested in predicting a new point$y(s_{0})$. In this case, I should have a matrix$H^{*}(\phi)$defined as</p> <p>$$H^{*}(\phi) = \left(\begin{array}{ccc}H(\phi) &amp; \boldsymbol{h} \\ \boldsymbol{h}&amp; \rho(0,\phi) \end{array}\right)$$</p> <p>in which$\boldsymbol{h}(\phi)$is a vector$\rho(s_{0} - s_{j};\phi)$. Therefore, we can construct a vector (without rearrangement):</p> <p>$${\bf{W^{*}}} = \left({\bf{W}}(s_{1}), \ldots, {\bf{W}}(s_{n}), {\bf{W}}(s_{0})\right)^{T} \sim N\left(\begin{array}{ccc}\boldsymbol{1}_{n+1} \otimes \left( \begin{array}{ccc} \mu_{1} \\ \mu_{2} \end{array} \right)\end{array}, H(\phi)^{*}\otimes T\right)$$</p> <p>And now I just rearrange to get a joint distribution$\left(\begin{array}{ccc} {\bf{X}} \\ x(s_0) \\{\bf{Y}} \\ y(s_0)\end{array} \right)$and obtain the conditional$p(y(s_0)\mid x_0, {\bf{X}}, {\bf{Y}})$. </p> <p>Is this correct? </p> <p>Question 2: For predicting, the paper I'm reading indicates that I must use this conditional distribution$p(y(s_0)\mid x_0, {\bf{X}}, {\bf{Y}})$and obtain a posterior distribution$p(\mu, T, \phi\mid x(s_0), {\bf{Y}}, {\bf{X}})$, but I'm not sure how to obtain the posterior distribution for the parameters. Maybe I could use the distribution$\left(\begin{array}{ccc}{\bf{X}} \\ x(s_0)\\ {\bf{Y}}\end{array}\right)$that I think is exactly the same as$p({\bf{X}}, x(s_0), {\bf{Y}}\mid\mu, T, \phi)$and then simply use Bayes' theorem to obtain$p(\mu, T, \phi\mid {\bf{X}}, x(s_0), {\bf{Y}}) \propto p({\bf{X}}, x(s_0), {\bf{Y}}\mid\mu, T, \phi)p(\mu, T, \phi)$</p> <p>Question 3: At the end of the subchapter, the author says this:</p> <blockquote> <p>For prediction, we do not have${\bf{X}}(s_0)$. This does not create any new problems as it may be treated as a latent variable and incorporated into$\bf{x}'$This only results in an additional draw within each Gibbs iteration and is a trivial addition to the computational task.</p> </blockquote> <p>What does that paragraph mean? </p> <p>By the way, this procedure can be found in <a href="http://sankhya.isical.ac.in/search/64a2/64a2020.pdf?origin=publication_detail" rel="nofollow">this paper</a> (page 8), but as you can see, I need a bit more of detail.</p> <p>Thanks!</p> https://stats.stackexchange.com/q/65693 9 Confusion related to Gibbs sampling user34790 https://stats.stackexchange.com/users/12329 2013-07-26T20:58:43Z 2014-11-16T10:24:23Z <p>I came across <a href="http://web.mit.edu/~wingated/www/introductions/mcmc-gibbs-intro.pdf">this article</a> where it says that in Gibbs sampling every sample is accepted. I am a bit confused. How come if every sample it accepted it converges to a stationary distribution.</p> <p>In general Metropolis Algorithm we accept as min(1, p(x*)/p(x)) where x* is the sample point. I assume that x* points us to a position where the density is high so we are moving to the target distribution. Hence I suppose that it moves to the target distribution after a burn in period.</p> <p>However, in Gibbs sampling we accept everything so even though it may take us to a different place, how can we say that it converges to the stationary/target distribution</p> <p>Suppose we have a distribution$p(\theta) = c(\theta)/Z$. We cannot calculate Z. In metropolis algorithm we use the term$c(\theta^{new})/c(\theta^{old})$to incorporate the distribution$c(\theta)$plus the normalizing constant Z cancels out. So it's fine</p> <p>But in Gibbs sampling where are we using the distribution$c(\theta)$</p> <p>For eg in the paper <a href="http://books.nips.cc/papers/files/nips25/NIPS2012_0921.pdf">http://books.nips.cc/papers/files/nips25/NIPS2012_0921.pdf</a> its given</p> <p>so we don't have the exact conditional distribution to sample from, we just have something that is directly proportional to the conditional distribution</p> <p><img src="https://i.stack.imgur.com/8yH9o.png" alt="enter image description here"></p> https://stats.stackexchange.com/q/223952 8 Rao-Blackwellization of Gibbs Sampler mscnvrsy https://stats.stackexchange.com/users/123353 2016-07-15T14:59:21Z 2016-07-19T07:24:29Z <p>I am currently estimating a stochastic volatility model with Markov Chain Monte Carlo methods. Thereby, I am implementing Gibbs and Metropolis sampling methods.<br /><br />Assuming I take the mean of the posterior distribution rather than a random sample from it, is this what is commonly referred to as <i>Rao-Blackwellization</i>?<br /><br /> Overall, this would result in taking the mean over the means of the posterior distributions as parameter estimate.<br /><br /></p> https://stats.stackexchange.com/q/136579 8 Gibbs Sampler transition kernel user2016445 https://stats.stackexchange.com/users/68266 2015-02-06T13:41:30Z 2015-02-06T17:31:34Z <p>Let$\pi$be the target distribution on$(\mathbb{R}^d,\mathcal{B}(\mathbb{R^d}))$which is absolutely continuously wrt to the$d$-dimensional Lebesgue measure, i.e :</p> <p>$\pi$admits a density$\pi(x_1,...,x_d)$wrt to$\lambda^d$with $$\lambda^d(dx_1,...,dx_d) = \lambda(dx_1) \cdot \cdot \cdot \lambda (dx_d)$$</p> <p>Let us assume that the full conditionals$\pi_i(x_i|x_{-i})$from$\pi$are known. So the transition kernel of the Gibbs-Sampler is clearly the product of the full conditionals from$\pi$. </p> <p><strong><em>Is the transition kernel absolutely continuously wrt to the$d$-dimensional Lebesgue measure too ?</em></strong> </p> https://stats.stackexchange.com/q/392607 8 Would an "importance Gibbs" sampling method work? user118967 https://stats.stackexchange.com/users/190575 2019-02-15T05:45:33Z 2019-03-18T08:04:59Z <p>I suspect this is a fairly unusual and exploratory question, so please bear with me.</p> <p>I am wondering if one could apply the idea of importance sampling to Gibbs sampling. Here's what I mean: in Gibbs sampling, we change the value of one variable (or block of variables) at a time, sampling from the conditional probability given the remaining variables.</p> <p>However, it may not be possible or easy to sample from the exact conditional probability. So instead we sample from a proposal distribution <span class="math-container">$q$</span> and use, for example, Metropolis-Hastings (MH).</p> <p>So far, so good. But here is a divergent path: what happens if we, instead of using MH, use the same idea used in importance sampling, namely we sample from <span class="math-container">$q$</span> and keep an importance weight <span class="math-container">$p/q$</span> of the current sample?</p> <p>In more detail: assume we have variables <span class="math-container">$x_1,\dots,x_n$</span> and a factored distribution <span class="math-container">$\phi_1,\dots,\phi_m$</span> so that <span class="math-container">$p \propto \prod_{i=1}^m \phi_i$</span>. We keep the proposal probability <span class="math-container">$q_i$</span> used to sample the current value of each variable <span class="math-container">$x_i$</span>. At each step we change a subset of the variables and update <span class="math-container">$p(x)/q(x)$</span> (only the factors of <span class="math-container">$p$</span> and <span class="math-container">$q$</span> that are affected). We take the samples and their importance weight to compute whatever statistic we are interested in.</p> <p>Would this algorithm be correct? If not, any clear reasons why not? Intuitively it makes sense to me as it seems to be doing the same thing importance sampling does, but with dependent samples instead.</p> <p>I did implement this for a Gaussian random walk model and observed that the weights become smaller and smaller (but not monotonically), so the initial samples end up having too much importance and dominate the statistic. I'm pretty certain the implementation is not buggy, because at each step I compare the updated weight to an explicit brute-force computation of it. Note that the weights do not go down indefinitely to zero, because they are <span class="math-container">$p/q$</span> where both <span class="math-container">$p$</span> and <span class="math-container">$q$</span> are products of a finite number of densities, and each sample is obtained from a Normal distribution that only rarely will be zero.</p> <p>So I am trying to understand why the weights go down like that, and whether this is a consequence of this method being actually not correct.</p> <hr> <p>Here's a more precise definition of the algorithm, applied to a Gaussian random walk on variables <span class="math-container">$X_1,\dots,X_n$</span>. The code follows below.</p> <p>The model is simply <span class="math-container">$X_i \sim \mathcal N(X_{i-1}, \sigma^2), i = 1,\dots,n$</span>, with <span class="math-container">$X_0$</span> fixed to <span class="math-container">$0$</span>.</p> <p>The weight of the current sample is <span class="math-container">$\frac{\prod_i p(x_i)}{\prod_i q(x_i)}$</span>, where <span class="math-container">$p$</span> are the Gaussian densities and <span class="math-container">$q$</span> are the distributions from which the current values have been sampled. Initially, we simply sample the values in a forward manner, so <span class="math-container">$q = p$</span> and the initial weight is <span class="math-container">$1$</span>.</p> <p>Then at each step I choose <span class="math-container">$j \in \{1,\dots,n\}$</span> to change. I sample a new value <span class="math-container">$x'_j$</span> for <span class="math-container">$X_j$</span> from <span class="math-container">$\mathcal N(X_{j-1},\sigma^2)$</span>, so this density becomes the new used proposal distribution for <span class="math-container">$X_j$</span>.</p> <p>To update the weight, I divide it by the densities <span class="math-container">$p(x_j | x_{j-1})$</span> and <span class="math-container">$p(x_{j+1} | x_j)$</span> of old value <span class="math-container">$x_j$</span> according to <span class="math-container">$x_{j-1}$</span> and <span class="math-container">$x_{j+1}$</span>, and multiply by the densities <span class="math-container">$p(x'_j | x_{j-1})$</span> and <span class="math-container">$p(x_{j+1} | x'_j)$</span> of new value <span class="math-container">$x'_j$</span> according to <span class="math-container">$x_{j-1}$</span> and <span class="math-container">$x_{j+1}$</span>. This updates the numerator <span class="math-container">$p$</span> of the weight.</p> <p>To update the denominator <span class="math-container">$q$</span>, I multiply the weight by the old proposal <span class="math-container">$q(x_j)$</span> (thus removing it from the denominator) and divide it by <span class="math-container">$q(x'_j)$</span>.</p> <p>(Because I sample <span class="math-container">$x'_j$</span> from the normal centered on <span class="math-container">$x_{j-1}$</span>, <span class="math-container">$q(x'_j)$</span> is always equal to <span class="math-container">$p(x'_j | x_{j-1})$</span> so they cancel out and the implementation does not actually use them).</p> <p>Like I mentioned before, in the code I compare this incremental weight computation to the actual explicit computation just to be sure.</p> <hr> <p>Here's the code for reference.</p> <pre class="lang-java prettyprint-override"><code>println("Original sample: " + currentSample); int flippedVariablesIndex = 1 + getRandom().nextInt(getVariables().size() - 1); println("Flipping: " + flippedVariablesIndex); double oldValue = getValue(currentSample, flippedVariablesIndex); NormalDistribution normalFromBack = getNormalDistribution(getValue(currentSample, flippedVariablesIndex - 1)); double previousP = normalFromBack.density(oldValue); double newValue = normalFromBack.sample(); currentSample.set(getVariable(flippedVariablesIndex), newValue); double previousQ = fromVariableToQ.get(getVariable(flippedVariablesIndex)); fromVariableToQ.put(getVariable(flippedVariablesIndex), normalFromBack.density(newValue)); if (flippedVariablesIndex &lt; length - 1) { NormalDistribution normal = getNormalDistribution(getValue(currentSample, flippedVariablesIndex + 1)); double oldForwardPotential = normal.density(oldValue); double newForwardPotential = normal.density(newValue); // println("Removing old forward potential " + oldForwardPotential); currentSample.removePotential(new DoublePotential(oldForwardPotential)); // println("Multiplying new forward potential " + newForwardPotential); currentSample.updatePotential(new DoublePotential(newForwardPotential)); } // println("Removing old backward potential " + previousP); currentSample.removePotential(new DoublePotential(previousP)); // println("Multiplying (removing from divisor) old q " + previousQ); currentSample.updatePotential(new DoublePotential(previousQ)); println("Final sample: " + currentSample); println(); // check by comparison to brute force calculation of weight: double productOfPs = 1.0; for (int i = 1; i != length; i++) { productOfPs *= getNormalDistribution(getValue(currentSample, i - 1)).density(getValue(currentSample, i)); } double productOfQs = Util.fold(fromVariableToQ.values(), (p1, p2) -&gt; p1*p2, 1.0); double weight = productOfPs/productOfQs; if (Math.abs(weight - currentSample.getPotential().doubleValue()) &gt; 0.0000001) { println("Error in weight calculation"); System.exit(0); } </code></pre> https://stats.stackexchange.com/q/50170 8 Bayesian estimation of Dirichlet distribution parameters Solmaz https://stats.stackexchange.com/users/20939 2013-02-17T07:57:25Z 2013-11-30T09:41:07Z <p>I want to estimate parameters of Dirichlet mixture models using Gibbs sampling and I have some questions about that:</p> <ol> <li><p>Is a mixture of Dirichlet distributions equivalent to a Dirichlet process? What is their main differences if is not?</p></li> <li><p>Also, if I want to estimate a single Dirichlet distribution's parameters, which distribution for parameters should be selected as priors in Bayesian framework?</p></li> </ol> <p>In all of the papers I found an estimation of a multinomial distribution using Dirichlet priors. I need estimation of a Dirichlet distribution using multinomial priors, perhaps. </p> <p>Is the posterior function also in the form of DIRICHLET(α+N) similar to the case “estimation of multinomial distribution using Dirichlet priors”? as the multiplication of probability density function for iid samples are not considered in the definition of the likelihood function. I again cannot understand why.</p> <p>e.g. as expressed in this paper: <a href="http://www.stat.ufl.edu/~aa/cda/bayes.pdf" rel="nofollow">http://www.stat.ufl.edu/~aa/cda/bayes.pdf</a> or <a href="http://research.microsoft.com/en-us/um/people/minka/papers/minka-multinomial.pdf" rel="nofollow">http://research.microsoft.com/en-us/um/people/minka/papers/minka-multinomial.pdf</a></p> <hr> <p>so thanks for your attention</p> <p>my data is Hyperion (a kind of hyperspectral remote sensing imagery) and i want to perform hyperspectral unmixing using mixture of Dirichlet sources and i will apply Gibbs sampling method for parameter estimation. my data is in dimension (614*512*224) which is commonly available AVIRIS sensor data for Cuprite Nevada district and is almost 200MB. also this data is available via (<a href="http://aviris.jpl.nasa.gov/data/free_data.html" rel="nofollow">http://aviris.jpl.nasa.gov/data/free_data.html</a>). unfortunately i don't know how can i sent my data. </p> <p>i just ask you to help me in statistical modelling tasks for my PHD thesis. i will be so grateful if you help me to solve my confusions in modelling.</p> <p>all the best solmaz</p> https://stats.stackexchange.com/q/100591 8 Can I subsample a large dataset at every MCMC iteration? alberto https://stats.stackexchange.com/users/36312 2014-05-30T10:55:52Z 2017-09-04T07:55:39Z <p><strong>Problem:</strong> I want to perform a Gibbs sampling to infer some posterior over a large dataset. Unfortunatelly, my model is not very simple and thus sampling is too slow. I would consider variational or parallel approaches, but before going that far...</p> <p><strong>Question:</strong> I would like to know whether I could randomly sample (with replacement) from my dataset at every Gibbs iteration, so that I have less instances to learn from at every step.</p> <p>My intuition is that even if I change the samples, I would not be changing the probability density and therefore the Gibbs sample should not notice the trick. Am I right? Are there some references of people having done this?</p> https://stats.stackexchange.com/q/411247 8 What is the correct way to write the elastic net? gbd https://stats.stackexchange.com/users/185448 2019-06-02T20:53:28Z 2019-06-03T04:24:04Z <p>I am confused about the correct way to write the elastic net. After reading some research papers there seems to be three forms</p> <p>1) <span class="math-container">$\exp\{-\lambda_1|\beta_k|-\lambda_2\beta_k^2\}$</span></p> <p>2) <span class="math-container">$\exp\{-\frac{(\lambda_1|\beta_k|+\lambda_2\beta_k^2)}{\sqrt{\sigma^2}}\}$</span></p> <p>3) <span class="math-container">$\exp\{-\frac{(\lambda_1|\beta_k|+\lambda_2\beta_k^2)}{2\sigma^2}\}$</span></p> <p>I just don't understand the correct way to add <span class="math-container">$\sigma^2$</span>. Is any of the above expressions correct?</p> https://stats.stackexchange.com/q/45946 7 Generating samples from Gibbs sampling Manish https://stats.stackexchange.com/users/14433 2012-12-14T22:58:04Z 2016-06-22T13:27:50Z <p>I am quite new to sampling. I am doing Gibbs sampling for a Bayesian network. I am aware about the algorithm for the Gibbs sampling but there's one thing I am not able to understand. </p> <p>For example let's assume you have 3 variables$x_1, x_2, x_3$means you need to generate sample of the form$(x_1,x_2,x_3)$. And suppose that all three parameters can take 3 values 0,1 or 2. Now, let's take initial value of the sample as (2,1,1) and compute the$\mathbb P(x_1=0 | x_2=1,x_3=1)$,$\mathbb P(x_1=1 |x_2=1,x_3=1)$and$\mathbb P(x_1=2 |x_2=1,x_3=1)$. </p> <p>Now, I am selecting the largest probability among these 3 and assigning that value to$x_1$and using it for the next sample. I want to ask if this method for generating model is correct or I am doing something wrong.</p> https://stats.stackexchange.com/q/56543 7 PyMC: how can I define a function of two stochastic variables, with no closed-form distribution? roger_ https://stats.stackexchange.com/users/24323 2013-04-19T01:17:53Z 2016-04-04T19:23:15Z <p>I'm learning PyMC and basically I have a random variable$Z = X + Y$where (say)$X \sim \mathrm{Normal}(\theta_X)$and$Y \sim \mathrm{Lognormal}(\theta_Y)$and$Z$has no simple closed-form distribution. Now I have observations$z_i,\,i=1...N$of$Z$and I want to infer$\theta_X$and$\theta_Y$. What's the most straight-forward way of doing this with PyMC?</p> <p>If I had the distribution of$Z$available, then I think I could do:</p> <pre><code>Z = DistZ('Z', param_x=theta_x, param_y=theta_y, value=z, observed=True) </code></pre> <p>and then do inference, but I don't know <code>DistZ</code>. It's also easy to define the sum as:</p> <pre><code>@pymc.deterministic def z_sum(x=Y, Y=y): return x + y </code></pre> <p>but then I don't think I can define an observed deterministic function.</p> <p>I <em>think</em> I could do something like:</p> <pre><code>@pymc.stochastic(observed=True) def z_sum(value=z, x=X, y=Y): def logp(z, x): # return log-likelihood </code></pre> <p>but I'm not clear on the details. I do know the joint likelihood$\mathcal{L}(z, x)$, but I was hoping it wouldn't be needed.</p> <p>I was able to do this with a custom Gibbs sampler (using the joint likelihood), but I'm looking for a more "elegant" solution with PyMC.</p> <hr> <p><strong>EDIT</strong>: found a <a href="http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-faqs/#q14" rel="nofollow">similar question</a> in the BUGS FAQ that says functions of random variables aren't supported. Not sure if that applies to PyMC, and what the standard approach is.</p> https://stats.stackexchange.com/q/24209 7 Gibbs sampler from conditional distribution user1061210 https://stats.stackexchange.com/users/9177 2012-03-06T20:32:18Z 2013-05-17T17:38:10Z <p>I am trying to propose Gibbs sampling with the density below,</p> <p>$$p(y_1,y_2,y_3)\propto \exp [-({{y}_{1}}+{{y}_{2}}+{{y}_{3}}+{{\theta}_{12}{y_1}{y_2}}+{{\theta }_{13}{y_3}{y_1}}+{{\theta }_{23}{y_2}{y_3}})]$$</p> <p>where,$({{y}_{1}},{{y}_{2}},{{y}_{3}})\in R_{+}^{3}$and${{\theta }_{ij}}=i+j$</p> <p>How do I find the full conditional distribution?</p> <p>And then, I'll generate sample$\{(y_{1}^{i},y_{2}^{i},y_{3}^{i})\}$for$i=1,...n$.</p> <p>I understand Gibbs sampling, sample one variable while keep others fixed. </p> https://stats.stackexchange.com/q/185311 7 Why does sampling from the posterior predictive distribution$p(x_{new} \mid x_1, \ldots x_n)$work without having to average out the integral? user1398057 https://stats.stackexchange.com/users/53410 2015-12-06T15:34:37Z 2019-03-03T23:29:37Z <p>In a Bayesian model, the posterior predictive distribution is usually written as:</p> <p>$$p(x_{new} \mid x_1, \ldots x_n) = \int_{-\infty}^{\infty} p(x_{new}\mid \mu) \ p(\mu \mid x_1, \ldots x_n)d\mu$$</p> <p>for a mean parameter$\mu$. Then, inside most books, such as this link:</p> <p><a href="http://people.stat.sc.edu/Hitchcock/stat535slidesday18.pdf">Sampling MCMC</a></p> <p>It is claimed that it is often easier to sample from$p(x_{new} \mid x_1, \ldots x_n)$using Monte Carlo methods. Commonly, the algorithm is to:</p> <p>for$j=1 \ldots J$:</p> <p>1) Sample$\mu^{\ j}$from$p(\mu \mid x_1, \ldots x_n)$then</p> <p>2) Sample$x^{\ * j}$from$p(x_{new} \mid \mu^{\ j})$. </p> <p>Then,$x^{\ * 1}, \ldots, x^{\ * J}$will be an iid sample from$p(x_{new} \mid x_1, \ldots x_n)$. </p> <p>What confuses me is the validity of this technique. My understanding is that Monte Carlo approaches will approximate the integral, so in this case, why do the$x^{\ * j}$'s each constitute a sample from$p(x_{new} \mid x_1, \ldots x_n)$? </p> <p>Why isn't is the case that the average of all those samples instead will be distributed as$p(x_{new} \mid x_1, \ldots x_n)$? I am under the assumption that I am creating a finite partition to approximate the integral above. Am I missing something? Thanks!</p> https://stats.stackexchange.com/q/68424 7 Gibbs sampling to produce posterior pdf TeTs https://stats.stackexchange.com/users/29021 2013-08-27T11:55:51Z 2013-08-28T13:38:25Z <p>Suppose we have the following classical normal linear regression model:</p> <p>$$y_i = \beta_1 x_{1i} + \beta_2x_{2i} + \beta_3x_{3i} + e_i$$</p> <p>where$e_{i} \sim iid.N(0, \sigma^2)$for all$i = 1, 2, \cdots, n$and$x_{1i} = 1$for all$i = 1, 2, \cdots, n$.</p> <p>Assume that we have known data values for both$x_{2i}$and$x_{3i}$for all$i = 1, 2, \cdots, n$. Defining$\boldsymbol{\beta} = (\beta_1, \beta_2, \beta_3)'$and assuming a non-informative prior of the form$p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$, then we can show that the conditional posterior pdf for$\boldsymbol{\beta}$and$\sigma$, that is,$p(\boldsymbol{\beta}|\sigma, \mathbf{y})$and$p(\sigma|\boldsymbol{\beta},\mathbf{y})$are normal and inverted gamma, respectively.</p> <p>The question is: Use a Gibbs Sampler and estimate the posterior pdf of the parameter function:$\displaystyle \psi = \frac{\beta_2 + \beta_3}{\sigma^2}$.</p> <p>Now I have run a Gibbs sampler (in R) and after a burn in period of 100 draws, I have obtained 1000 draws of$(\boldsymbol{\beta}^{(i)}, \sigma^{(i)})$, that is, I have a sample of$(\beta_1^{(i)}, \beta_2^{(i)}, \beta_3^{(i)}, \sigma^{(i)})$for$i = 1, 2, \cdots, 1000$, how can I use these draws to produce an estimate of the posterior pdf of$\psi$? In other words, how can I estimate$p(\psi|\mathbf{y})$?</p> <hr> <p>EDIT: I'm new to MCMC logarithms. I do understand what you mean but I am still not too sure how to use it in the context of this question. From what I've learnt so far, say we have$\boldsymbol{\theta} = (\theta_1, \theta_2)'$and$p(\boldsymbol{\theta}|\mathbf{y}) = p(\theta_1, \theta_2|\boldsymbol{y})$is the joint posterior, then the marginal posterior of$\theta_1$is given by$p(\theta_1|\mathbf{y}) = \int_{\theta_1} p(\theta_1|\theta_2,\mathbf{y})p(\theta_2|\mathbf{y})d\theta_2$, now say I have a sample of$M$draws of$(\theta_1^{(i)}, \theta_2^{(i)})$from$p(\theta_1, \theta_2|\boldsymbol{y})$, then to estimate$p(\theta_1|\mathbf{y})$, we use the following sample mean:$\widehat{p(\theta_1|\mathbf{y})} = \frac{1}{M} \sum_{i=1}^M p(\theta_1|\theta_2^{(i)}, \mathbf{y}) $, that is, we estimate the marginal density by averaging over the conditional densities. How can I implement that here?</p> <hr> <p>EDIT2: After the help of @Tomas and @BabakP, I have tried to code the problem into R myself (I'm quite new to R, only learnt the language in the last couple of weeks). The following is my code: (Note: just for practice, I purposely went the other route of drawing from each conditional rather than @BabakP's more efficient decomposition)</p> <pre><code>########################################################### ######## Generation of data = y ####### ########################################################### true_beta = matrix(c(2,0.5,0.7),3,1) ## True value of beta (vector) used in the data generating process (dgp) true_sig = 0.3 ## True value of sig used in the data generating process (dgp) ## Setting the random number seed set.seed(123456, kind = NULL, normal.kind = NULL) ## Number of observations nobs=20 ## Generating the values for x1, x2, x3 and y. Using Gaussian distribution here for convenience. x1 = matrix(1,nobs,1) x2 = matrix(c(rnorm(nobs,0,1)),nobs,1) x3 = matrix(c(rnorm(nobs,0,1)),nobs,1) u = matrix(c(rnorm(nobs,0,1)),nobs,1)%*%true_sig x = cbind(x1,x2,x3) y = x%*%true_beta + u ########################################################### ### Specification of sample statistics ########################################################### invxx = solve(t(x)%*%x) betahat = invxx%*%t(x)%*%y df = nobs - 3 sighat = sqrt((t(y-x%*%betahat)%*%(y-x%*%betahat))/df) ############################################################### ### ### Generation of draws of beta/sigma/psi via Gibbs sampling ### ############################################################### B = 100+1 ## Extra +1 since the Gibbs sampler coded below starts at j=2, hence burn-in period is from j=2 to j=101 M = 1000 ## Number of draws after burn in period. repl = B+M ## Total number of Gibbs iterations including burn in and number of draws after burn in period. ## Dimensioning the matrix of beta/sigma/psi draws ## betav = matrix(0,repl,3) sigv = matrix(0,repl,1) psiv = matrix(0,repl,1) ## Generating underlying normal random draws to be used in Gibbs algorithm stnbeta = matrix(c(rnorm((repl*3),0,1)),repl,3) ## Specifying starting value for Gibbs chain sigv = true_sig ## Gibbs sampling for(j in 2:repl){ ## Generation of beta (vector) via its known multivariate normal conditional meanbeta = betahat varbeta = (sigv[j-1])^2*invxx betagen = chol(varbeta)%*%stnbeta[j,] + meanbeta betav[j,] = betagen ## Generation of sigma via its known inverted gamma conditional a = t(y-x%*%betav[j,])%*%(y-x%*%betav[j,]) zvec = matrix(c(rnorm(nobs,0,1)),nobs,1) chisq = t(zvec)%*%zvec sigv[j] = sqrt(a/chisq) ## Draws of psi from its marginal posterior pdf psiv[j] = (betav[j,2]+betav[j,3])/(sigv[j])^2 } ################################################################################# ## Estimation of the marginal posterior pdf of psi using kernel density smoothing ################################################################################# ## Remove burn in period of psiv which is from j=1 to j=101 to create the vector psi which conists of 1000 draws after the burn in period psi = psiv[102:repl] ## Histogram of draws of psi hist(psi,30,ylab="Percentage frequencies",xlab=expression(psi), main=expression(paste("Histogram of draws of ",psi))) ## kernel density plot of the draws of psi d=density(psi) plot(d,main=expression(paste("Kernel density estimate of marginal of ",psi)), xlab=expression(psi),ylab=expression(paste("p(",psi,"|y)"))) </code></pre> <p>My question is now how can I use the 1000 draws of psi from the above code to estimate the 5% quantile of$p(\psi|\mathbf{y})\$?</p>