As far as I can tell, MCMC integration (e.g. VEGAS) is performed by sampling from a distribution proportional to $f(x)$ using MCMC, then building a density estimator $g(x)$ using these samples (for example, in case of VEGAS it is a factorized histogram $g(x_1, \dots, x_n) = g(x_1) \dots g(x_n)$), and then using this estimator to perform importance sampling and reweighting of values of $f(x)$, like $$I(f) = \mathbb E_{x \sim g(x)} \frac{f(x)}{g(x)}.$$ Correct me if I am wrong. It is not exactly clear to me whether one can use values $f_i$ at points $x_i \sim h(x) = \frac{f(x)}{\int f(x) dx}$ to directly approximate the final integral without introducing a density model $g(x)$ that unavoidably introduces an error into the estimate, especially in higher dimensions where density estimation is hard?
Let me give an example of an error I am talking about. Consider that $f(x)$ is just a standard 1D Gaussian pdf so it sums up to one. We receive three $(x_1, x_2, x_3)$ points near zero and corresponding values $(f_1, f_2, f_3)$ and then estimate $g_{K}(x)$ from these $x_i$ points using a kernel density estimator. Depending on the choice of the kernel width $K$ the density estimate in these points $g(x_i)$ would fluctuate significantly, i.e. for small $K$ values of $g(x_i)$ would be high and therefore $\sum_i \frac{f_i}{g(x_i)}$ would be low.
Here $K$ corresponds to our assumption about the smoothness of $f$ and if it is wrong, we get wrong results. The approximation used in VEGAS assumes that the function $f: R^n \to R$ factors well into the product of scalar functions.
Is there a way to avoid these possibly wrong assumptions? Assume that we already have plenty of values $(f_1, \dots , f_n)$ measured at points $(x_1, \dots , x_n)$ and we know that $x_i$ follow the distribution $x_i \sim h(x) \propto f(x)$?