Summary: Let \(\mu\) be a Borel probability measure on a compact path-connected metric space \((X, \rho)\) for which there exist constants \(c, \beta \geq 1\) such that \(\mu(B) \geq c r^\beta\) for every open ball \(B \subset X\) of radius \(r > 0\). For a class of Lipschitz functions \(\Phi : [0, \infty) \to \mathbb{R}\) that are piecewise within a finite-dimensional subspace of continuous functions, we prove under certain mild conditions on the metric \(\rho\) and the measure \(\mu\) that for each positive integer \(N \geq 2\), and each \(g \in L^\infty (X, d \mu)\) with \(\| g \|_\infty = 1\), there exist points \(y_1, \ldots, y_N \in X\) and real numbers \(\lambda_1, \ldots, \lambda_N\) such that for any \(x \in X\), \[ \left| \int\limits_X \Phi (\rho (x, y)) g(y)\, \mathrm{d}\mu (y) - \sum\limits_{j = 1}^N \lambda_j \Phi (\rho (x, y_j)) \right| \leqslant CN^{-\frac{1}{2} - \frac{3}{2\beta}} \sqrt{\log N}, \] where the constant \(C > 0\) is independent of \(N\) and \(g\). In the case when \(X\) is the unit sphere \(\mathbb{S}^d\) of \(\mathbb{R}^{d + 1}\) with the usual geodesic distance, we also prove that the constant \(C\) here is independent of the dimension \(d\). Our estimates are better than those obtained from the standard Monte Carlo methods, which typically yield a weaker upper bound \(N^{-\frac{1}{2}} \sqrt{\log N}\).