I have found that sampling only becomes intuitive once it stops looking like a mysterious numerical trick. At its core, it is a way of replacing a calculation we cannot do with a calculation that we can: draw a lot of sensible random values, and let their average tell us something useful. That sounds almost too simple, but it gets us surprisingly far.
The interesting part is that this simple idea is enough to handle distributions for which even the normalising constant is out of reach. Let us get there gradually.
Probability and Inference
Before getting to sampling, I find it useful to remember where the distributions come from. The two probability rules I return to most often are the sum rule and the product rule. For two events $A$ and $B$,
\[P(A \cup B) = P(A) + P(B) - P(A \cap B)\]and
\[P(A \cap B) = P(A \mid B)P(B).\]The product rule is especially useful because we can write it in either direction. Equating the two ways of writing the joint probability gives Bayes’ rule:
\[p(\theta \mid D) = \frac{p(D \mid \theta)p(\theta)}{p(D)}, \qquad p(D) = \int p(D \mid \theta)p(\theta)\,d\theta.\]Here $D$ is the data and $\theta$ is whatever the model is trying to learn. If we only want one best answer, maximum likelihood chooses one value,
\[\hat{\theta}_{\mathrm{ML}} = \arg\max_{\theta} p(D \mid \theta),\]which is often a perfectly reasonable thing to do. The Bayesian view instead keeps the whole posterior $p(\theta \mid D)$. I like thinking of it as refusing to throw uncertainty away too early: we can say not only what seems most plausible, but also how unsure we are about it.
The slightly annoying term is $p(D)$, sometimes called the evidence. It asks us to add up the likelihood over every possible value of $\theta$. For a one-dimensional toy problem that may be manageable; for a useful model, it is usually not. Sampling is one of the ways around that obstacle.
A Small Bayesian Update
Before worrying about difficult posteriors, it helps to see what an update feels like. Imagine an unknown coin bias $\theta$. A Beta prior, $\theta \sim \operatorname{Beta}(\alpha, \beta)$, is convenient here: it lets us treat $\alpha$ and $\beta$ roughly like prior heads and prior tails. After observing $h$ heads and $t$ tails, the posterior is simply
\[\theta \mid D \sim \operatorname{Beta}(\alpha + h, \beta + t).\]The two boxes on the left are your assumptions before data. Try a prior of 8 heads and 2 tails, for example: you start out rather convinced the coin favours heads. Then change the observed heads and tails. The dotted curve is what you started with; the rust curve is where the evidence moves you. This is the basic pattern we want to retain when the update is no longer analytically friendly.
Interactive
Update a coin-bias belief
First choose your prior belief, then add the coin flips you observed. You can drag the sliders or type a number.
Monte Carlo Estimation
Now forget coins for a moment. Suppose we want the area under a function $f(x)$ between $a$ and $b$. The familiar approach is integration. Monte Carlo does something more indirect: sample points $x_i$ uniformly in that interval, look at their heights $f(x_i)$, and average them.
If $X \sim \operatorname{Uniform}(a,b)$, then
\[\int_a^b f(x)\,dx = (b-a)\,\mathbb{E}[f(X)] \approx \frac{b-a}{N}\sum_{i=1}^N f(x_i).\]There is nothing magical hidden in the approximation. The average height of enough randomly chosen points approaches the average height of the function, and average height times width is area. The law of large numbers does the rest.
Try changing the function and the number of samples. Each dot shows one random place where we looked at the function. The estimate changes with every resample, but it settles down as $N$ grows. That noisy convergence is a useful intuition for most Monte Carlo algorithms.
Interactive
Estimate an area with random points
Change the sample count, then press resample a few times. Watch the estimate become less jumpy as you add more dots.
More generally, if $X \sim p(x)$, then
\[\mathbb{E}_{p}[f(X)] = \int f(x)p(x)\,dx \approx \frac{1}{N}\sum_{i=1}^N f(x_i).\]This is the whole promise and the whole catch: the approximation is easy once we have samples from $p$. What if $p$ itself is a complicated posterior?
For easy distributions, inverse-transform sampling gives one answer. Draw $U \sim \operatorname{Uniform}(0,1)$ and set $X = F^{-1}(U)$, where $F$ is the cumulative distribution function. But many useful $F^{-1}$ functions have no neat closed form. That is why we need other ways of turning easy randomness into harder randomness.
Sampling from a Posterior
Returning to Bayes’ rule, the posterior is proportional to the product of a likelihood and a prior:
\[p(\theta \mid D) \propto p(D \mid \theta)p(\theta).\]Even when we cannot calculate the evidence $p(D)$, we can often calculate the right-hand side. It is useful to call it an unnormalised score:
\[\tilde{p}(\theta \mid D) = p(D \mid \theta)p(\theta).\]It does not integrate to one, so it is not quite a probability density yet. Still, it tells us which values of $\theta$ look better relative to one another. Many sampling algorithms only need those relative scores, so the troublesome normalising constant never has to be computed explicitly.
Rejection Sampling
Rejection sampling is the most visual example. Start with a proposal distribution $q(\theta)$ that is easy to draw from. Then place an envelope above the unnormalised target, choosing $M$ large enough that
\[\tilde{p}(\theta \mid D) \leq Mq(\theta) \quad \text{for every } \theta.\]Draw a candidate $\theta \sim q(\theta)$ and a uniform number $u \sim \operatorname{Uniform}(0,1)$. Keep the candidate when
\[u \leq \frac{\tilde{p}(\theta \mid D)}{Mq(\theta)}.\]Geometrically, we are throwing points into the envelope and keeping the ones under the target curve. The accepted points have the target distribution. The demo uses a uniform proposal for a Beta-shaped target: rust points survive, grey points do not. Adjust the two shape parameters to see why a loose envelope makes the method inefficient.
Interactive
Keep points under a target curve
Change the shape of the target, then draw a new batch. Every point begins equally likely in the rectangle; only the points below the solid curve are kept.
The method is wonderfully direct, but it has an obvious weakness. If the proposal puts little mass where the posterior is large, the envelope must be tall and most candidates will be rejected. Better proposals, importance sampling, and Markov chain Monte Carlo are all different ways of dealing with that same problem.