Algorithm to sample from unknown probability distribution, given its projections? The Cramér–Wold theorem states that, if $X$ is a random variable living in $\mathbb{R^d}$, then the distributions of all one-dimensional projections of $X$ uniquely determine the distribution of $X$.
But can this be made algorithmic? I have an algorithm that approximately finds the distribution of any one-dimensional projection of a certain random variable $X$. Can I use this algorithm to efficiently approximately sample from $X$?
If yes, then -- hopefully this is not asking too much -- are there error bounds on how good the sampling is? Ideally I would want a statement similar to "If one can approximate the distribution of any projection of $X$ to within $\ell_1$ distance at most $\varepsilon$ from its true distribution, then one can efficiently sample from a distribution with at most $\ell_1$ distance [some function of $d$ and $\varepsilon$] from $X$."
 A: This is too long for a comment. As noted above, this is Radon inversion problem, which can be reduced to the Fourier inversion problem. As is noted in Wikipedia, the inverse problem is ill-posed, which roughly speaking means that no error bounds of the type you want are possible. Let me try and give a simple explanation why.
Suppose we know the density of $X\cdot v$ for any unit-norm vector norm $v$, with $L_2(\mathbb{R})$ error $\varepsilon$. This means that we know its Fourier transform $t\mapsto\mathbb{E}(\exp(it(X\cdot v)))$, also with $L_2(\mathbb{R})$ error $\varepsilon$ (maybe up to $2\pi$ factor). This means that we know the d-dimensional Fourier transform $\hat{X}:v\mapsto \mathbb{E}(\exp(i(X\cdot v)))$ with an error of order $\varepsilon$ in the weighted $L^2(\mathbb{R}^d,\omega)$ space, with weight $\omega(v)=1/|v|^{d-1}$ (the Jacobian factor coming from passing to spherical coordinates). But this can correspond to arbitrarily large error in the usual $L_2(\mathbb{R}^d)$ sense, if the error is on high frequencies. When we recover $X$ from its Fourier transform, this arbitrarily large error is preserved.
Update: The following argument gives an easy error bound involving Wasserstein distance. Let $\hat{Y}$ be the Fourier transform recovered from the data, and $\hat{X}$ the actual Fourier transform, so that $||\hat{X}-\hat{Y}||_{L^2(\mathbb{R}^d,\omega)}\leq c_d \varepsilon$. Multiply $\hat{Y}$ by $\exp(-\delta|v|^2/2)$ for a small $\delta$ and take the inverse Fourier transform. What we will get is $f_{X+\sqrt{\delta}Z}+f_{\text{err}}$, where $f_{X+\sqrt{\delta}Z}$ is the density of $X+\sqrt{\delta}Z$, with $Z$ standard Gaussian, and $f_{\text{err}}$ is the inverse Fourier transform of $\exp(-\delta|v|^2/2)(\hat{X}-\hat{Y})$. Now, $X+\sqrt{\delta}Z$ is within $C\sqrt{\delta}$ in the $L_1$ Wasserstein distance from $X$. On the other hand, the maximum of $\exp(-\delta|v|^2)|v|^{d-1}$ is of order $\delta^{-(d-1)/2}$. Therefore,
$$||e^{-\frac{\delta|v|^2}{2}}(\hat{X}-\hat{Y})||_{L_2(\mathbb{R}^d)}\leq C\delta^{(1-d)/4}||(\hat{X}-\hat{Y})||_{L_2(\mathbb{R}^d,\omega)}\leq C\delta^{(1-d)/4}\varepsilon.$$
So, if we know the densities of all the projections up to $L_2$ error $\varepsilon$, then we can recover a distribution that is within $C\delta^{(1-d)/4}\varepsilon$ in $L_2$ distance from a distribution that is within $C\cdot \sqrt{\delta}$ in the $L_1$ Wasserstein distance from the actual distribution, where $\delta$ is in our disposal.
