What is the best way to sample from joint distributions with independent marginals? Suppose we have an $n$-dimensional joint distribution where all its marginals are independent. That is, if the joint density function is $p(x_1,\ldots,x_n)$, then $p(x_1,\ldots,x_n)=p_1(x_1)\cdots p_n(x_n)$, where $p_1,\ldots,p_n$ are marginal densities, and all these marginals are known and quite simple.
Now we want to get samples ${\bf x}_k=(x_{k,1},\ldots,x_{k,n}),k=1,\ldots,m$ from $p$ where each sample is an $n$-dimensional vector, and the objective is to use these samples to estimate the expectation of $E(h(X))$ where $h$ is a real-valued function, X is a random variable distributed according to $p$, by computing the mean of $h({\bf x}_1),\ldots,h({\bf x}_m)$.
For this purpose $m$ is better to be a very large number. Anyone knows what is the most efficient way to do so besides MCMC?

The brutal naive way is to sample $x_{k,i}$ from $p_i$ for every $i=1,\ldots,n$ for every $k$. This is not desirable when $m$ is large.
We exclude MCMC because the cost is too high for our application. All marginals are independent, known, and simple. We don't want to involve this heavy machinery.
A possible alternative may be that, after we sample ${\bf x}$ from $p$, we then take turns to re-sample each dimension of ${\bf x}$ from corresponding marginal, but we are not sure if this is correct (i.e. if the estimate is unbiased, will the convergence will be much slower). Anyone can help prove or show counterexample of this?

 A: The efficiency you can reach will depend on the type of function you have for $h({\bf X})$ as well as on the implementation and CPU-cost of the various computations, and of course the problem at hand.
The resampling method that you suggest is absolutely fine. The additional vectors thus obtained, are however not fully independent. So you get a trade of between a slower convergence due to that, versus a faster computation because you only need to generate and evaluate the probabilities of the few dimensions you sample, while most of the components of the vector ${\bf x}_k$ remain the same.
How well this works and how much better the estimated average will be cannot be told in advance without knowledge of the function $h({\bf X})$ and the various probabilities.
The ultimate generalisation of this idea is to take every possible combination of the components of the $k$ vectors you have. So in stead of
$$
E(h({\bf X})) = \sum_{k=1}^m p({\bf x}_k) h({\bf x}_k) = \sum_{k=1}^m \left( \prod_{i=1}^n p_i(x_{k,i}) \right) h({\bf x}_k)
$$
you could take $m^n$ different (but not independent) vectors:
$$
E(h({\bf X})) = \sum_{i_1=1}^m p_1(x_{i_1,1})\dots  \sum_{i_n=1}^m p_n(x_{i_n,n}) h(\{x_{i_1,1},\dots,x_{i_n,n}\}) 
$$
which allows for very efficient calculations and is the most information you could get out of the $m$ independent vectors that were generated innitially. This of course only works if the function $h({\bf X})$ consists of additive and multiplicative terms of subsets of the various components.
Note that the estimated error in $E(h({\bf X}))$ should be determined taking into account that there are only $m$ independent vectors.
