Why does Metropolis Hastings algorithm need/make a sequence of interdependent samples, i.e., a Markov chain in order to generate a posterior distribution?
17655
Maybe the concept, why it's used, and an example.
31732
My background is computer science. I am fairly new to monte carlo sampling methods and, although I understand the math, I have hard time coming up with intuitive examples for importance sampling. More precisely, could someone provide examples of:    an original distribution one cannot sample from but one can estimate an importance distribution which can be sampled from and adequate for this original distribution.
32455
I am trying to understand how to derive the posterior distribution of a parameter $\mu$ given data vector $z$, $P(\mu|z)$, where  $$ \mu \sim N(0,A) $$  and  $$  z|\mu \sim N(\mu,1). $$  Obviously from Bayes theorem  $$ P(\mu|z) = g(\mu) f(z|\mu) /f(z). $$  where $$ f(z)= \int g(\mu) f(z|\mu) d \mu . $$  I can write $g(\mu) f(z|\mu)$ as a pointwise product of normal densities:  $$  \frac{1}{2 \pi \sqrt{A}} \exp \bigg(-\frac{\mu^2+A(z-\mu)^2}{2A}\bigg) $$  the solution is given in Efron's book on Large Scale Inference, p. 2. :   $$ \mu|z \sim N(zB,B)$$ where $$B=\frac{A}{A+1} $$  I would appreciate advice on how to approach the problem (and the answer to this question is a proof). In particular I do not understand what to do with the integral in the numerator of Bayes theorem.   EDIT Following answer by @Neil_G I took a next approach:  We have:  \begin{align} g(\mu)   &amp;\propto \exp\bigg(- \frac{\mu^2}{2A}\bigg)  \\[8pt] f(z|\mu) &amp;\propto \exp\bigg(z \mu - \frac{1}{2}(z^2+\mu^2)\bigg) \end{align}  so \begin{align} P(\mu|z) &amp;\propto \exp\bigg(- \frac{1}{2A} \mu^2 + z\mu - \frac{1}{2} z^2 - \frac{1}{2} \mu^2\bigg)  \\[8pt]  &amp;= \exp \bigg(-\frac{1}{2B} (\mu^2+Bz^2-2Bz \mu)\bigg)  \\[8pt] &amp;\propto \exp \bigg(-\frac{1}{2B} (\mu^2-2Bz \mu)\bigg)  \\[8pt] &amp;\propto \exp \bigg(-\frac{ (\mu - Bz)^2}{2B}\bigg)  \end{align}  which completes the proof.
107354
I have been trying to get a sense of the different problems in frequentist settings where MCMC is used. I am familiar that MCMC (or Monte Carlo) is used in fitting GLMMs and in maybe Monte Carlo EM algorithms. Are there more frequentist problems where MCMC is used?
316
In the textbook "Causal Inference for Statistics" by Rubin and Imbens, the following argument is made on pg. 39:     "In part of this text we view our sample of size N as a random sample   from an infinite super-population. In that case we employ slightly   different formulations of the restric- tions on the assignment   mechanism. Sampling from the super-population generates a joint   sampling distribution on the quadruple of unit-level variables (Yi(0),   Yi(1), Wi, Xi), i = 1, . . . , N. More explicitly, we assume the   (Yi(0), Yi(1), Wi, Xi) are independently and identically distributed   draws from a distribution indexed by a global parameter."   I am wondering why the bolded portion is needed to be assumed. In the context of observational studies what happens if 1) we have dependence, and 2) if the tuple is not jointly identically distributed.   Is this a general assumption that can be relaxed? Where do we really need it? Specifically:  1) What happens if we have dependence?  2) What happens if they are NOT identically distributed?
8152
I am trying to compute this posterior distribution:  $$ (\theta|-)=\frac{\prod_{i=1}^{n}p_i^{y_i}(1-p_i)^{1-y_i}}{\sum_{\text{all}\,\theta,p_i|\theta}\prod_{i=1}^{n}p_i^{y_i}(1-p_i)^{1-y_i}} $$  The problem is that the numerator, which is the product of a bunch of $\text{Bernoulli}(p_i,y_i)$ probabilities is too small. (My $n$ is large, about 1500).   Hence, the posterior values for all $\theta$ all get calculated to be 0 (I am doing calculations in R).  To clarify, each $y_i$ has its own $p_i$, together these $p_i$'s make a vector of $n$ elements for $n$ $y$'s. Each $\theta$ has its own $n$-element vector of $p_i$.  EDIT: Adding a reproducing example (for the numerator)  p &lt;- sample(seq(0,1,by=0.01), 1500, replace=T) y &lt;- sample(c(0,1), 1500, replace=T) dbern(y, p) # 1500-element vector, each element is &lt; 1 prod(dbern(y, p)) # produce 0 exp(sum(log(dbern(y, p)))) # produce 0 since the sum is very negative
92523
Consider an irreducible but possibly periodic Markov chain on a finite state space with transition matrix $P$. We know there exists a unique stationary distribution $\pi$. If the Markov chain were aperiodic, we would have $P^n_{ij} \to \pi(j)$ as $n \to \infty$. This fails if the chain is periodic, but we have convergence of the Cesaro averages: $$\frac{1}{n}\sum_{k=1}^n P^k_{ij} \to \pi(j) \text{ as } n \to \infty.$$ Can anyone point me to a reference that states this fact? Every reference I've seen only considers convergence for aperiodic chains, or "fixes" periodicity by considering a lazy version of the chain. Alternatively, is there a simple way to obtain this result using the result for aperiodic chains?
42774
I know that in the case of Gaussian mixture, the "intuition" is that you're drawing from a PDF which itself is just a sum of weighted Gaussian PDFs.  I don't understand the intuition behind how the sum of two iid Gaussian RVs leads to a Gaussian RV. I know the maths behind it, but I'm having a hard time visualizing how the sum of two Gaussian RVs leads to a Gaussian RV (and not a mixture).
116483
I am new to MCMC and reading a  regarding Gibbs sampling. However, there are two parts in the paper I cannot understand and get stuck. The first part is equation 2.3 in page 168. It says that if $X^{'}_{j}$ is drawn from the process of Gibbs sampling, $X^{'}_{j}\sim f(x|Y^{'}_{j}=y^{'}_{j})$. It further claims under &quot;general conditions,&quot; the distribution of $X^{'}_{j}$ converges to $f(x)$ as $n\rightarrow\infty$. Can anyone help explain what these general conditions are and why/how this statement is true. The second part is at the beginning of page 169, expression 2.9 says $\hat{f(x)}=\frac{1}{m}\sum^{m}_{i=1}f(x|y_{i})$, and the theory behind this equation is $E[f(x|Y)]=\int f(x|y)f(y)dy$. This confuses me because my understanding about monte carlo integration is that for example, $E[X]=\frac{1}{m}\sum^{m}_{i=1}X_{i}$ is true due to the law of large number that $\frac{X_{i}}{m}$ converges to $X_{i}\times p(X_{i})$ when $m\rightarrow\infty$. However, I don't see how equation 2.9 can be explained by the concept of monote carlo integration. Can anyone explain the intuition or math behind this equation? Thanks! EDIT 1 (My Answer/Thought for Part 2) I have thought about this question, and since no one responds yet, I will try my answer first and see whether someone has a better explanation. For the second part, instead of using $f(x|y_{i})$, let's focus on the density of a specific value $x^{a}$, which is $f(x^{a}|y_i)$. For simplicity, now assume $y_i$ is randomly drawn from a density $f(y)$ instead of $f(y|x_i)$ but the result could be generalized to $f(y|x_i)$. Both $f(x|y_i)$ and $f(y)$ are know distribution. Now we first randomly draw a $y_i$ from $f(y)$ and once we have $y_i$, we further randomly draw a $x_i$ from $f(x|y_i)$. Suppose we do this process $n$ time, and $n\rightarrow\infty$. Given a specific $y^{a}$ and $x^{a}$, we do know the value of $f(x^{a}|y^{a})$, which is just the density of $x^{a}$ cnditional on $y^{a}$. In other words, by Monte Carlo integration, for a given set of $y_i$'s and $x^{a}$, we can dervive that $\frac{1}{n}\sum^{n}_{i=1}f(x^{a}|y_i)\approx\int f(x^{a}|y_i)f(y_i)dy_i=\int f(x^{a}|y)f(y)dy=\int \dfrac{f(x^{a},y)}{f(y)}f(y)dy=\int f(x^{a},y)dy=f(x^{a})$ We can do the same process to dervie $f(x^{b})$, $f(x^{c})$, and etc. In the end, we will have $f(x^{a})$, $f(x^{b})$, $f(x^{c})$, ..., which all together constrcut $f(x)$. Therefore, $\hat{f(x)}=\frac{1}{m}\sum^{m}_{i=1}f(x|y_{i})$.
71759
I have this confusion related to the predictive distribution of gaussian process. I was reading this      I didn't get how the integration gave that result. What is P(u*|x*,u). Also how come the covariance of the posterior distribution is $\sigma^2(\sigma^2I+K)^{-1}K$
84532
I'm taking an Intro to Bayes course and I'm having some difficulty understanding predictive distributions. I understand why they are useful and I'm familiar with the definition, but there are some things I don't quite understand.   1) How to get the right predictive distribution for a vector of new observations  Suppose that we have built a sampling model $p(y_i | \theta)$ for the data and a prior $p(\theta)$. Assume that the observations $y_i$ are conditionally independent given $\theta$.   We have observed some data $\mathcal{D} = \{y_1, y_2, \, ... \, , y_k\}$, and we update our prior $p(\theta)$ to the posterior $p(\theta | \mathcal{D})$.   If we wanted to predict a vector of new observations $\mathcal{N} = \{\tilde{y}_1, \tilde{y}_2, \, ... \, , \tilde{y}_n\}$, I think we should try to get the posterior predictive using this formula $$ p(\mathcal{N} | \mathcal{D}) = \int p(\theta | \mathcal{D}) p ( \mathcal{N} | \theta) \, \mathrm{d} \theta = \int p(\theta | \mathcal{D}) \prod_{i=1}^n p(\tilde{y}_i | \theta) \, \mathrm{d} \theta, $$ which is not equal to $$ \prod_{i=1}^n  \int p(\theta | \mathcal{D}) p(\tilde{y}_i | \theta) \, \mathrm{d} \theta, $$ so the predicted observations are not independent, right?  Say that $\theta | \mathcal{D} \sim$ Beta($a,b$) and $p(y_i | \theta) \sim$ Binomial($n, \theta$) for a fixed $n$. In this case, if I wanted to simulate 6 new $\tilde{y}$, if I understand this correctly, it would be wrong to simulate 6 draws independently from the Beta-Binomial distribution that corresponds to the posterior predictive for a single observation. Is this correct? I don't know how to interpret that the observations are not independent marginally, and I'm not sure I understand this correctly.  Simulating from posterior predictives  Many times when we simulate data from the posterior predictive we follow this scheme:  For $b$ from 1 to $B$:  1) Sample $\theta^{(b)}$ from $p(\theta | \mathcal{D})$.  2) Then simulate new data $\mathcal{N}^{(b)}$ from $p(\mathcal{N} | \theta^{(b)})$.   I don't quite know how to prove this scheme works, although it looks intuitive. Also, does this have a name? I tried to look up a justification and I tried different names, but I had no luck.  Thanks!
84330
Suppose that I have two univariate marginal distributions, say $F$ and $G$, which I can simulate from. Now, construct their joint distribution using a , denoted $C(F,G;\Sigma)$. All the parameters are known.  Is there a non-MCMC method for simulating from this copula?
109602
This post relates to my original question , but this time focusing on a more fundamental misunderstanding of how MCMC actually works.  When using Bayesian MCMC for parameter inference with a single model, it is usual to ignore the "evidence" integral and write:  $$P(\theta | D) \propto P(D | \theta) P(\theta)$$  For numerical reasons, it is also common to work with logarithms:  $$ln[P(\theta | D)] \propto ln[P(D | \theta)] + ln[P(\theta)]$$  I've used this successfully on a number of occasions, typically using Python and . However, for reasons explained in my , I've recently been wondering about "normalising" my likelihood by dividing it by the number of observations (I realise this probably sounds odd, but please see my other question for why I thought it was worth a try).  My likelihood assumes i.i.d. Gaussian errors, so my log-likelihood is just the sum of the individual data likelihoods. My prior is uniform and I'm currently using an improper log-prior function, which simply returns zero if the parameter values are "within range" and '-inf' otherwise.   I initially assumed that dividing by the number of data points wouldn't make any difference, because the RHS is still proportional to the log-posterior:  $$ln[P(\theta | D)] \propto \frac{1}{n}ln[P(D | \theta)] + ln[P(\theta)]$$  (My prior only ever returns 0 or '-inf', so multiplying that by $1/n$ makes no difference).  However, when I re-run my MCMC, I get very different results. As a sense check, I tried a simple OLS linear regression using the same code.   With my original (not normalised) likelihood, I get sensible results: the MCMC recovers the "true" parameters and the 95% CI broadly agrees with the Frequentist results from 'statsmodels':    However, with the "normalised" likelihood, I get very different (and obviously wrong) results:    I think I understand why this is, at least from the point of view of a simple Metropolis algorithm. The probability of accepting the proposed jump, $x_{i+1}$, compared to the current location, $x_i$, is:  $$\alpha = \frac{P(x_{i+1})}{P(x_i)}$$  Which, re-writing in terms of my log-posterior, $L$, becomes:  $$\alpha = \frac{exp(L(x_{i+1}))}{exp(L(x_i))} = exp(L(x_{i+1}) - L(x_i))$$  On the other hand, for my "normalised" version, all the likelihoods are divided by $n$, so this becomes:  $$\alpha = [exp(L(x_{i+1}) - L(x_i))]^\frac{1}{n}$$  which is different, and might explain why my uncertainty bounds suddenly become much larger. (I understand that emcee is more sophisticated than just basic Metropolis, but I guess a similar principle applies?).  My question:  What am I missing/misunderstanding to reconcile these two results, please? On the one hand, it's OK to ignore the "normalising constant" in Bayes' equation, because anything proportional to the log-posterior is sufficient. On the other hand, multiplying the log-posterior by a constant $\frac{1}{n}$ seems to have a dramatic effect on my results?  Thank you!
47005
The Bayesian bootstrap was introduced by  as a Bayesian analog of the original bootstrap. Given dataset $X=\{x_1, \dots, x_N\}$, instead of drawing weights $\pi_{n}$ from the discrete set $\left\{0, \frac{1}{N}, \ldots, \frac{N}{N}\right\},$ the Bayesian approach treats the vector of weights $\pi$ as unknown parameters and derives a posterior distribution for $\pi$. Rubin (1981) used an improper, non-informative prior, $\prod_{i=1}^{n} \pi_{i}^{-1},$ which when combined with the multinomial likelihood, leads to a Dirichlet(1,...,1) distribution for the posterior distribution of $\pi$. In other words, our prior is \begin{equation}     p(\boldsymbol{\pi}) =Dirichlet(\boldsymbol{\alpha}), \quad \text{with}\ \boldsymbol{\alpha} = [0,\dots,0]. \end{equation} and posterior is \begin{equation}     p(\boldsymbol{\pi}|\boldsymbol{x}) =Dirichlet(\boldsymbol{\alpha}), \quad \text{with}\ \boldsymbol{\alpha} = [1,\dots,1]. \end{equation} Now my questions are:  I was asked the following questions which I was not able to answer: How can you have a posterior distribution that a/ does not depend on data and b/ is a uniform distribution?  Are both the prior and the posterior non-informative? I understand that the posterior is a uniform distribution, which is non-informative. Also, I see that the prior is referred to as non-informative prior. Does that mean it's flat?   I believe that section 5 in  addresses these questions, but I do not comprehend that discussion. Any clarification or pointing out what I may be misunderstanding would be appreciated.  EDIT: I just noticed one more issue when computing the posterior. Let $d=\left(d_{1}, \ldots, d_{K}\right)$ be the vector of all possible distinct values of $X$, and let $\pi=\left(\pi_{1}, \cdots, \pi_{K}\right)$ be the associated vector of probabilities $$ P\left(X=d_{i} \mid \pi\right)=\pi_{i}, \quad \sum \pi_{i}=1 $$ Let $x_{1}, \ldots, x_{n}$ be an i.i.d. sample from the equation above and let $n_{i}$ be the number of $x_{j}$ equal to $d_{i}$. If we use improper prior above over the sampling weights $\pi$, we can compute the posterior over $\pi$ \begin{align*}     p(\boldsymbol{\pi}|X) &amp;\propto p(X|\boldsymbol{\pi})p(\boldsymbol{\pi})\\     &amp; \propto \prod_{i}\pi_i^{n_i}\prod_{i}\pi_{i}^{\alpha_i-1}\\     &amp; \propto \prod_{i}\pi_i^{n_i}\prod_{i}\pi_{i}^{-1}\\     &amp; \propto \prod_i\pi_i^{n_i-1}. \end{align*} How does this yield a flat Dirichlet posterior? Are we assuming $n_i=1$ for $i=1,\dots,K$? In that case, is the vector of all possible observations $d=\left(d_{1}, \ldots, d_{K}\right)$ (the original sample that we resample from) our observation?
96082
Prove that if $\mu$ is a stationary distribution for an irreducible Markov chain on $S$, then $\mu(j) &gt; 0 \ \forall j \in S$. (only using the fact that Markov chain is irreducible and that $\mu P = \mu$)   How can I do that? I have used this fact while proving other limit theorems but can't think about how to do it.
84149
I work in the field of data mining and have had very little formal schooling in statistics. Lately I have been reading a lot of work  that focuses on Bayesian paradigms for learning and mining, which I find very interesting.   My question is (in several parts), given a problem is there a general framework by which it is possible to construct a statistical model? What are the first things you do when given a dataset of which you'd like to model the underlying process? Are there good books/tutorials out there that explain this process or is it a matter of experience? Is inference in the forefront of your mind when constructing your model or do you first aim to describe the data before you worry about how to use it to compute?  Any insight would be greatly appreciated! Thanks.
43314
Robby McKilliam says in a comment to  post:     It should be pointed out that, from the frequentists point of view, there is no reason that you can't incorporate the prior knowledge into the model. In this sense, the frequentist view is simpler, you only have a model and some data. There is no need to separate the prior information from the model   Also, , @jbowman says that frequentists use regularization by a cost/penalty function, while bayesians can make this a prior:     Frequentists realized that regularization was good, and use it quite commonly these days - and Bayesian priors can be easily interpreted as regularization.   So, my question is, can frequentists in general incorporate into their models what Bayesians specify as priors? Taking the regularization as an example, is the cost/penalty function really integrated into the model, or is this a purely artificial means of adjusting the solution (as well as making it unique)?