A typical introduction to the Metropolis--Hastings algorithm, and hence to Markov chain Monte Carlo techniques in general, starts with the following assumptions on some probability distribution $P(x)$ which we'd like to sample from.  We are perfectly able to evaluate $P(x)$, or perhaps more accurately some multiple $C \cdot P(x)$ for some fixed $C$, though I don't quite see why you wouldn't just normalise and put $C = 1$. However, it is very hard to sample from $P(x)$; so hard, that we'd rather sample from some other distribution thousands of times for the sake of getting a sample of $P(x)$ via Markov chain methods.  I'd like to see an actual, real-life example in which such a scenario comes up. How can it be easy to evaluate $P(x)$ but hard to sample from it? I'd be interested in seeing a real-life situation involving say weather, diseases, neurons, air flow, whatever it is, that leads to a situation in which MCMC methods are a natural thing to ask for.
69973
My understanding is that when using a Bayesian approach to estimate parameter values:   The posterior distribution is the combination of the prior distribution and the likelihood distribution. We simulate this by generating a sample from the posterior distribution (e.g., using a Metropolis-Hasting algorithm to generate values, and accept them if they are above a certain threshold of probability to belong to the posterior distribution). Once we have generated this sample, we use it to approximate the posterior distribution, and things like its mean.   But, I feel like I must be misunderstanding something. It sounds like we have a posterior distribution and then sample from it, and then use that sample as an approximation of the posterior distribution. But if we have the posterior distribution to begin with why do we need to sample from it to approximate it?
72273
Given a random variable $X$ which follows a normal distribution $X \sim (\mu,\sigma^2)$ with unknown mean and known variance. Say we have $n$ observations $y_i$ where $1\le i \le n$ and each $y_i$ is sampled from a normal distribution $Y_i \sim N(x,\sigma_i^2)$ (not i.i.d.) where $x$ is sampled from $X$ and it is the same for all $\{Y_i\}$, but $\sigma_i$ is different for different $Y_i$.  The question is: how to compute the posterior distribution $X$ given the observations p(X|{y_i})? The question is different from the answer listed in  because $y_i$ is not sampled from i.i.d. distribution.
77854
I'm studying 'Introduction to Statistical Learning' by James, Witten, Hastie, Tibshirani.   In page 139, of their book, they began by introducing Bayes' Theorem $p_k(X)=P(Y=k|X=x) = \dfrac{\pi_kf_k(x)}{\sum_{l=1}^k \pi_l f_l(x)}$. $\pi$ is not mathematical constant, but denotes the prior probability. Nothing is strange in this equation.   The book claims that it wants to obtain an estimate for $f_k(x)$ that it can plug into the equation given above. To estimate $f_k(x)$, it assumes that is normal. In the one dimensional setting, $f_k(x)=\dfrac{1}{\sqrt{2\pi}\sigma}\exp(-\dfrac{1}{2\sigma^2}(x-\mu_k)^2)$, where $\mu_k$ and $\sigma^2_k$ are the mean and variance for the $k$th class. It is assumed that $\sigma^2_1 = \sigma^2_2 = \cdots = \sigma^2_K$. (I began to get confused from the last statement.)   Plugging $f_k$ into $p_x$, you have this quite messy equation (1) :   $$p_x(k)=\dfrac{\pi_k \frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{1}{2\sigma^2}(x-\mu_k)^2)}{\sum_{l=1}^K \pi_l \frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{1}{2\sigma^2}(x-\mu_l)^2)}. $$  Again, no surprises here since it is just substitution.     The Bayes' Classifier involves assigning an observation to the class for which equation (1) is the largest. Taking the log of Equation (1) and rearranging terms, it is not hard to show that this is equivalent to assigning the observation to the class for which the following is the largest:  $$\delta_k(x)=x \cdot \dfrac{\mu_k}{\sigma^2} - \dfrac{\mu_k^2}{2\sigma^2} + \log(\pi_k) $$  Question: I do not understand where this came from, and what it means. I tried doing the log of equation and it doesn't become this. Are we taking the derivative somewhere here, since this is the largest observation?
110924
I was wondering if someone could explain the following passage from my textbook.    In the third paragraph from the bottom it is stated that the simulated values come from the marginal distributions. How is this so? We simulate the value x2 from the condition density f(x2|x1)..  Also, there is mention of independence and again because of the dependence of the value of x2 on x1 and so forth, I do not see how this occurs.   Would greatly appreciate your help
40102
I have been trying to understand the idea of conjugate priors in Bayesian statistics for a while but I simply don't get it. Can anyone explain the idea in the simplest possible terms, perhaps using the "Gaussian prior" as an example?
84580
I am getting familiar with Bayesian statistics by reading the book Doing Bayesian Data Analysis, by John K. Kruschke also known as the "puppy book". In chapter 9, hierarchical models are introduced with this simple example: \begin{align} y_{ji}   &amp;\sim {\rm Bernoulli}(\theta_j)  \\ \theta_j &amp;\sim {\rm Beta}(\mu\kappa, (1-\mu)\kappa)  \\ \mu      &amp;\sim {\rm Beta}(A_\mu, B_\mu)   \\ \kappa   &amp;\sim {\rm Gamma}(S_\kappa, R_\kappa) \end{align} and the Bernoulli observations are 3 coins, each 10 flips. One shows 9 heads, the other 5 heads and the other 1 head.  I have used pymc to infer the hyperparamteres.  with pm.Model() as model: # define the          mu = pm.Beta('mu', 2, 2)     kappa = pm.Gamma('kappa', 1, 0.1)     # define the prior     theta = pm.Beta('theta', mu * kappa, (1 - mu) * kappa, shape=len(N))     # define the likelihood     y = pm.Bernoulli('y', p=theta[coin], observed=y)      # Generate a MCMC chain     step = pm.Metropolis()     trace = pm.sample(5000, step, progressbar=True)     trace = pm.sample(5000, step, progressbar=True)   burnin = 2000  # posterior samples to discard thin = 10  # thinning  pm.autocorrplot(trace[burnin::thin], vars =[mu, kappa])   My question is regarding the autocorrelation. How shall I interpret autocorrelation? Could you please help me to interpret the autocorrelation plot?    It says as the samples gets further from each other the correlation between them reduces. right? Can we use this to plot to find the optimum thinning? Does the thinning affect the posterior samples? after all, what is the use of this plot?
66664
The second paragraph of the  section of the Wikipedia article on resampling statistics, the values of a confidence interval for a p-value from a MC sampling is given:     After $N$ random permutations, it is possible to obtain a confidence interval for the p-value based on the Binomial distribution. For example, if after $N = 10000$ random permutations the p-value is estimated to be $\hat{p}=0.05$, then a 99% confidence interval for the true $p$ (the one that would result from trying all possible permutations) is $[0.044, 0.056]$.   What's the formula to calculate it?
66536
Background  I am designing a Monte Carlo simulation that combines the outputs of series of models, and I want to be sure that the simulation will allow me to make reasonable claims about the probability of the simulated outcome and the precision of that probability estimate.  The simulation will find the probability that a jury drawn from a specified community will convict a certain defendant. These are the steps of the simulation:   Using existing data, generate a logistic probability model (M) by regressing “juror first ballot vote” on demographic predictors. Use Monte Carlo methods to simulate 1,000 versions of M (i.e., 1000 versions of the coefficients for the model parameters). Select one of the 1,000 versions of the model (Mi). Empanel 1,000 juries by randomly selecting 1,000 sets of 12 “jurors” from a “community” (C) of individuals with specified demographic characteristic distributions. Deterministically calculate the probability of a first ballot guilty vote for each juror using Mi. Render each "juror’s" probable vote into a determinate vote (based on whether it is greater or less than randomly selected value between 0-1). Determine each "jury’s" “final vote” by using a model (derived from empirical data) of the probability a jury will convict, conditional on the proportion of jurors voting for conviction on the first ballot. Store the proportion of guilty verdicts for the 1000 juries (PGi). Repeat steps 3-8 for each of the 1,000 simulated versions of M. Calculate the mean value of PG and report that as the point estimate of the probability of conviction in C. Identify the 2.5 &amp; 97.5 percentile values for PG and report that as 0.95 confidence interval.   I am currently using 1,000 jurors and 1,000 juries on the theory that 1,000 random draws from a probability distribution&mdash;demographic characteristics of C or versions of M&mdash;will fill out that distribution.   Questions  Will this allow me to accurately determine the precision of my estimate? If so, how many juries do I need to empanel for each PGi calculation to cover C's probability distribution (so I avoid selection bias); may I use fewer than 1,000?  Thank you so much for any help!
91987
I am taking a course on Monte Carlo methods and we learned the Rejection Sampling (or Accept-Reject Sampling) method in the last lecture. There are a lot of resources on the web which shows the proof of this method but somehow I am not convinced with them.  So, in the Rejection Sampling, we have a distribution $f(x)$ which is hard to sample from. We choose a easy-to-sample distribution $g(x)$ and find a coefficient $c$ such that $f(x) \leq cg(x)$. Then we sample from $g(x)$ and for each draw, $x_i$, we also sample a $u$ from a standard uniform distribution $U(u|0,1)$.  The sample $x_i$ is accepted if it is $cg(x_i)u \leq f(x_i)$ and rejected otherwise.  The proofs which I came across usually just show that $p(x|Accept) = f(x)$ and stop there.   What I think about this process is that we have a sequence of variables $x_1,Accept_1,x_2,Accept_2,...,x_n,Accept_n$ and a $x_i,Accept_i$ pair corresponds to our i.th sample ($x_i$) and whether it is accepted ($Accept_i$). We know that each $x_i,Accept_i$ pair is independent of each other, such that:  $P(x_1,Accept_1,x_2,Accept_2,...,x_n,Accept_n) = \prod\limits_{i=1}^n P(x_i,Accept_i)$  For a $(x_i,Accept_i)$ pair we know that $P(x_i) = g(x_i)$ and $P(Accept_i|x_i) = \frac{f(x_i)}{cg(x_i)}$. We can readily calculate $p(x_i|Accept_i)$ but I don't understand how it suffices as a proof. We need to show that the algorithm works, so I think a proof should show that the empricial distribution of the accepted samples converge to $f(x)$ as $n\rightarrow\infty$. I mean, with $n$ being the number of all accepted and rejected samples:  $\frac{Number \hspace{1mm} of \hspace{1mm} samples \hspace{1mm} with \hspace{1mm} (A \leq x_i \leq B)}{Number \hspace{1mm} of \hspace{1mm} accepted \hspace{1mm} samples} \rightarrow \int_A^B f(x)dx$    as    $n\rightarrow\infty$.  Am I wrong with this thought pattern? Or is there a connection between the common proof of the algorithm and this?   Thanks in advance
95334
Assume we have $N$ random variables $X_1, \ldots, X_N$ with known (posterior) distributions that are easy to sample from. For simplicity, assume that I am interested in the expected value of the ten largest of these random variables. There will be some uncertainty surrounding this expected value as there is uncertainty with respect to the rank of a given $X_i$ in case the distributions of $X_1, \ldots, X_N$ overlap. I want to use simulation based methods. My approach is to simulate $M$ samples from the distributions of $X_1, \ldots, X_N$, where in each iteration, I construct a binary indicator of length $N$ indicating whether random variable $i$ is within the set of the ten largest samples or not. Averaging over these $M$ draws gives me the probability that random variable $i$ is member of the top ten group. Denote this probability as $\pi_i$. Define a binary random variable that indicates membership in the top ten group as $S_i$. $S_i$ follows a Bernoulli distribution s.t. $S_i \sim Ber(\pi_i)$. In a second step, I draw $L$ samples from the distribution of $S_i$ and from the distributions of $X_1, \ldots, X_N$. Given a draw of $S_i$, I can select the top ten out of the samples of $X_1, \ldots, X_N$ in each iteration and average the respective samples within the top ten group. I save this average, so the result after $L$ iterations will be the distribution of the average outcome within the top ten group, taking into account all relevant sources of uncertainty. My questions are: Is this simulation approach valid in general? From a Bayesian perspective, is this simulation valid or do I need to treat additional quantities (e.g. $S_i$ or $\pi_i$) as random variables that need a prior?
8970
In several different contexts we invoke the  to justify whatever statistical method we want to adopt (e.g., approximate the binomial distribution by a normal distribution). I understand the technical details as to why the theorem is true but it just now occurred to me that I do not really understand the intuition behind the central limit theorem.  So, what is the intuition behind the central limit theorem?   Layman explanations would be ideal. If some technical detail is needed please assume that I understand the concepts of a pdf, cdf, random variable etc but have no knowledge of convergence concepts, characteristic functions or anything to do with measure theory.
39755
I am studying different point estimate methods and read that when using MAP vs ML estimates, when we use a "uniform prior", the estimates are identical. Can somebody explain what a "uniform" prior is and give some (simple) examples of when the MAP and ML estimators would be the same?
123626
In a MCMC implementation of hierarchical models, with normal random effects and a Wishart prior for their covariance matrix, Gibbs sampling is typically used.  However, if we change the distribution of the random effects (e.g., to Student's-t or another one), the conjugacy is lost. In this case, what would be a suitable (i.e., easily tunable) proposal distribution for the covariance matrix of the random effects in a Metropolis-Hastings algorithm, and what should be the target acceptance rate, again 0.234?  Thanks in advance for any pointers.
41984
I'm reading  to work on the foundations of my understanding of Bayesian analysis. In the first chapter, about normal models, I'm having some trouble working the posterior distribution by myself.  For the model with known variance and conjugate prior for the mean $N(0, \sigma^2)$, the derivation seems straightforward (pages 29 and 30):  $$\pi(\mu\mid D_n) \propto \mathrm{exp}\{-\mu^2/2\sigma^2\}\mathrm{exp}\{-[n(\bar{x}-\mu)^2 + s^2]/2\sigma^2\}$$ $$\propto \mathrm{exp}\left\{-\frac{\mu^2 + n\mu^2 - 2n\mu\bar{x} +n\bar{x}^2 + s^2}{2\sigma^2}\right\}$$  $s^2$ drops off because it does not depend on the $\mu$ parameter; $n\bar{x}^2$ is almost what we need to complete the square, but not quite so, so I guess we also drop it off and substitute it for the proper term:  $$\propto \mathrm{exp}\left\{-\frac{(n+1)\mu^2 - 2n\mu\bar{x}}{2\sigma^2}\right\}$$ $$\propto \mathrm{exp}\left\{-\frac{(n+1)\left(\mu^2 - \frac{2n\mu\bar{x}}{n+1} + \frac{n^2\bar{x}^2}{(n+1)^2}\right)}{2\sigma^2}\right\}$$ $$\propto \mathrm{exp}\left\{-\frac{(n+1)\left(\mu - \frac{n\bar{x}}{n+1}\right)^2}{2\sigma^2}\right\}$$  So far, so good. But when it comes to the unknown mean, unknown variance model with prior $N(0, \sigma^2)$ and $\frac{1}{\sigma^2} \sim Exp(1)$, we can't drop off the $n\bar{x}^2$ term anymore because now the denominator inside the exponential function depends on the parameter $\sigma^2$.  So I have come to a similar point deriving the joint posterior for $\mu$ and $\sigma^2$:  $$\propto \mathrm{exp}\left\{-\frac{2 + \mu^2 + n\mu^2 - 2n\mu\bar{x} +n\bar{x}^2 + s^2}{2\sigma^2}\right\}\sigma^{-(n+5)}$$ $$\propto \sigma^{-1/2}\mathrm{exp}\left\{-\frac{(n+1)\mu^2 - 2n\mu\bar{x} +n\bar{x}^2}{2\sigma^2}\right\}\mathrm{exp}\left\{-\frac{2+s^2}{2\sigma^2}\right\}(\sigma^2)^{-(\frac{n+2}{2})-1}$$  The part to the right is already the core of the Inverse-Gamma function derived in the book, but I can't see how to complete the square for the normal density function without adding another term to the second Inverse-Gamma parameter. Indeed, in the general case, this parameter value should also have a $\frac{n\bar{x}^2}{n+1}$ term (as, e.g., ).   So I decided to compare both solutions by simulation using Metropolis-Hastings to obtain posterior draws of both parameters. Surprisingly (or not), the book solution uses the correct formula for the scale parameter for the Inverse-Gamma, $\frac{2 + s^2}{2}$, and if we add the term based on the sample mean, the density curve does not quite fit (Edit: not really, see below) . What am I missing here?    Edit: I'm starting to think that the book might indeed be mistaken. Completing the square in the first exponential makes it necessary to add a $\frac{n\bar{x}^2}{n+1}$ term that can only go in the second exponential term. So I tried my hand again with simulations, adding a JAGS model with the same definition to compare to my homemade solution (which had a coding error when I posted the question), and now it seems that the correct formula for the second Inverse-Gamma parameter is indeed $\frac{2 + s^2 + \frac{n\bar{x}^2}{n+1}}{2}$. Here's a plot with the simulated and theoretical densities for the posterior distribution (data = rnorm(100, 50, 10)):
44641
Assume I have a distribution governing the possible outcome from a single random variable X.  This is something like [0.1, 0.4, 0.2, 0.3] for X being a value of either 1, 2, 3, 4.  Is it possible to sample from this distribution, i.e. generate pseudo random numbers upon each of the possible outcomes given the probability of that outcome. So if I wanted to know what the probability of getting a 2 is, the sampling operation may return 0.34 or something like that.  The reason I ask is that I'm trying to implement an action selection policy for a reinforcement learning method based on a research paper. From what I gather from the paper, the author is able to sample the distribution by "mapping the uniform distribution U[0,1] through cumulative probability density functions obtained by adaptive numerical integration". From this he then samples the transition probabilities for each trial...  I would be grateful for any info on this...  Thanks in advance
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.
12675
Given a coin with unknown bias $p$, how can I generate variates — as efficiently as possible — that are Bernoulli-distributed with probability 0.5?  That is, using the minimum number of flips per generated variate.