Let $H = [-1, 1]^d$ be the $d$-dimensional hypercube, and let $\mu \in \text{int} H$.
Under these conditions, I can explicitly construct a tractable probability measure $P$, supported on on $H$, which has $\mu$ as its mean. For my purposes, tractability means:
(a) $P$ has a density, which can be written down and cheaply evaluated.
(b) $P$ can be exactly sampled from, efficiently.
This is a relatively simple task, since the problem basically decouples across dimensions; if you can find a probability measure $P_i$ on $[-1,1]$ with mean $\mu_i$, then setting $P = P_1 (x_1) \cdots P_d (x_d)$ suffices.
I've been using the family of random variables given by $P(x) \propto \exp( \langle \theta, x \rangle)$ to do this; the mean of each coordinate is then given by $m_i(\theta) = \coth \theta_i - \theta_i^{-1}$ which is surjective onto $(-1,1)$ and so one can always find a $\theta$ which gives you the mean you want.
Now, suppose I want to extend this, and construct a probability measure on $H$ which also has the covariance matrix I want, i.e. let $\Sigma$ be a general $d \times d$ positive semi-definite matrix. Then, I want all the above, plus:
\begin{align} \textbf{Var}_{X \sim P} [X] &= \Sigma \end{align}
This is more or less possible when $\Sigma$ is diagonal, but otherwise I haven't figured out how to do it. Of course, in the general case, I will need $\Sigma$ to be sufficiently small, as I can't have arbitrarily high variance when constrained to a cube,
I've considered sampling from a truncated Gaussian, i.e. sampling from $\mathcal{N}(\mu, \Sigma)$ until you land in $H$. However, this won't generally have a tractable density, due to the normalising constant. It also won't quite have the right moments, though I can basically forgive this; I ultimately don't mind too much if I only have
\begin{align} \mathbf{E}_{X \sim P} [X] &\approx \mu \\ \textbf{Var}_{X \sim P} [X] &\approx \Sigma \end{align}
with some control of the bias. I should perhaps highlight that it's not sufficient for my purposes to approximate $\Sigma$ by a diagonal matrix.
So, the concise version of my question is: is there a way of explicitly constructing probability measures on the hypercube which have the mean and covariance matrix which I want it to have?
There needn't be a 'canonical' / optimal' answer (whatever that might mean here); something functional will entirely suffice. It's also preferable to have a construction where it's easy to adapt things to sample from a range of different $(\mu, \Sigma)$.
