construction of a random measure with a given mean Let me first pose a trivial question. 

Given a Borel probability measure $\mu$ on the real line, is it possible to construct a purely atomic random measure $M$ whose mean is $\mu$?

The answer is obviously yes: take $M = \delta_{X}$, where $\delta_x$ is the Diract measure sitting at $x$ and $X$ is a random variable distributed according to $\mu$. In fact take $M$ to be the empirical measure $\frac{1}{n} \sum_{i=1}^n \delta_{X_i}$ where each $X_i \sim \mu$ suffices.
Now here is my question:

Given a Borel probability measure $\mu$ on the real line and $\delta > 0$, is it possible to construct a purely atomic random measure $M$ whose mean is $\mu$ and $d(M, \mu) < \delta$ almost surely, where $d$ is some metric on the space of probability measures (e.g. the Wasserstein distance, the Lévy–Prokhorov metric or the Kolmogorov distance, i.e., the sup-distance between distribution functions)?

Of course for an arbitrary distance this is not always possible. For example, if $d$ is the total variation distance and $\mu$ is atomless, then $d(M,\mu) = 2$ a.s. But for a weaker distance, will this be possible? The intuition is the following: consider the empirical measure $\frac{1}{n} \sum_{i=1}^n \delta_{X_i}$ where $X_i$ are iid generated from $\mu$. Then as $n\to\infty$, it will converge to the mean $\mu$ a.s. under those distances. Therefore within any $\delta$-ball centered at $\mu$, there are lots of atomic measures, i.e., $\mathbb{P} \{d(M,\mu) < \delta\}$ is very close to 1. But can we achieve exactly 1, i.e., can we construct an $M$ supported on those measures which are close to the desired mean?
 A: Yes -- the idea is to use the random delta-measure described in the first part of your question. However, in order to obtain a good approximation one has to subdivide the real line into small intervals (instead of taking the empirical averages over the whole line). 
For simplicity let us consider just the transportation metric $d$. Given a probability measure $\mu$ denote by $\tilde\mu$ the image of $\mu$ under the map $x\mapsto\delta_x$. Obviously, if the support of $\mu$ is contained in a length $\delta$ interval, then $d(\mu,\lambda)\le \delta$ for $\tilde\mu$-a.e. measure $\lambda$. Now subdivide $\mathbb R$ into intervals $I_i$ of length $\delta$, denote by $\mu_i$ the normalized restriction of $\mu$ onto $I_i$, take independent (actually, independence is not needed here) random measures $\lambda_i$ with distributions $\tilde\mu_i$, and put $\lambda= \sum \mu(I_i)\lambda_i$. Then $\mathbf E\lambda=\mu$ and almost surely $d(\mu,\lambda)\le\delta$.
A: Let $\{ X_n \}$ be an i.i.d. sequence with common law $\mu$. Denote the empirical measure by $L_n = \frac{1}{n} \sum_{i=1}^n \delta_{X_i}$. Define the stopping time $T = \min \{n: d(L_n, \mu) \leq \delta \}$. Then $T < \infty$ a.s. (which follows from Glivenko-Cantelli and LLN if we are considering transport metric like 2-Wasserstein). Therefore $L_T$ is a well-defined random measure, which satisfies $d(L_T, \mu) \leq \delta$ automatically. Moreover, $L_T$ is purely atomic with finitely many atoms almost surely. Observe that $L_T$ has mean $\mu$. To see this, note that for any positive Borel function $f$, 
$\mathbb{E}[\int f d L_T]$
$ = \sum_{n \geq 1} \frac{1}{n} \mathbb{E}[\sum_{i=1}^n f(X_i) | T = n] \mathbb{P}\{T = n\}$
$ = \sum_{n \geq 1} \mathbb{E}[f(X_1) | T = n] \mathbb{P}\{T = n\}$
$ = \mathbb{E}[f(X_1)]$
$ = \int f d \mu$
where the second equality is by symmetry.
