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$."



Algorithm to sample from unknown probability distribution, given its projections?
This is too long for a comment. As noted above, this is [Radon inversion problem](https://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform), 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.


