Writing the Lebesgue–Stieltjes integral as a sum of equidistributed Dirac delta measures Problem set up:
Let $f: [0, 1] \to \mathbb R$ be an absolutely continuous function (thus a fortiori of bounded variation) such that its total variation on any open interval $(a, b)$ is $b-a$.
We say a sequence $x_k$ taking values in $[0, 1]$ is equidistributed if for any open interval $(a, b)$ it holds that
$$\lim_{n \to \infty} \frac{1}{n} \#\{j \mathrel| j < n, x_j \in (a, b)\} = b - a,$$
where $\#$ denotes the cardinality of a finite set.
Question:
Let $x_j \in [0, 1]$ be an equidistributed sequence. Does it hold that there exists a sequence $\{\varepsilon_n\}_{n \in \mathbb N}$ taking values in $\{-1, 1\}$ such that for all Riemann integrable functions $g: [0, 1] \to \mathbb R$ such that the Lebesgue–Stieltjes integral $\int_{[0, 1]} g(x) df(x)$ exists, we have
$$\lim_{n \to \infty} \frac{1}{n} \sum_{k = 0}^{n-1} \varepsilon_k g(x_k) = \int_{[0, 1]} g(x) \ df(x)?$$
 A: $\newcommand{\ep}{\varepsilon}$This is at best a partial answer, with a simple result in the positive direction, which also highlights certain difficulties which may lead to a counterexample.
Let $p_+$ and $p_-$ denote, respectively, densities of the positive and negative parts of the Lebesgue--Stieltjes measure $df$, so that $p_+p_-=0$ almost everywhere (a.e.). The condition that the total variation of $f$ on any open interval $(a,b)$ with $0\le a<b\le1$ is $b-a$ means that $p_++p_-=1$ a.e. and hence a.e.
\begin{equation*}
    p_+=1_S,\quad p_-=1_{S^c}
\end{equation*}
for some Legesgue-measurable set $S\subseteq[0,1]$, where $S^c:=[0,1]\setminus S$.
So,
\begin{equation*}
    I:=\int_{[0,1]}g(x)\,df(x)=\int_0^1 h(x)\,dx,
\end{equation*}
where $g$ is any Riemann-integrable function on $[0,1]$ and
\begin{equation*}
    h:=g1_S-g1_{S^c}; 
\end{equation*}
then, of course, $g$ is Lebesgue integrable and hence so is $h$, so that the integral $I$ exists.
Suppose now for a moment that, possibly up to the a.e.-equivalence, the density $p_+=1_S$ is Riemann integrable. Then, without loss of generality, $h$ is Riemann integrable as well. So, letting
\begin{equation*}
    \ep_k:=
    \begin{cases}
    1&\text{ if }x_k\in S,\\
    -1&\text{ if }x_k\in S^c,
    \end{cases}\tag{1}
\end{equation*}
we get $\ep_k g(x_k)=h(x_k)$ and hence
\begin{equation*}
\frac1n\sum_{k=0}^{n-1}\ep_k g(x_k)
=\frac1n\sum_{k=0}^{n-1}h(x_k) \to\int_0^1 h(x)\,dx=\int_{[0,1]}g(x)\,df(x),
\end{equation*}
as desired.
This reasoning falls apart if $I_S$ is not a.e.-equivalent to a Riemann-integrable function -- which happens, for instance, when $S$ is the Smith--Volterra--Cantor (fat) set. Letting now $(x_j)$ be a specific, convenient equidistributed sequence with values in $S^c$ and probing a large enough collection of Riemann-integrable functions $g$, one might find a counterexample to the conjecture in general. Indeed, (1) strongly suggests that the $\ep_k$'s should mimic the set $S$, taking of course the membership of the equidistributed $x_k$'s in $S$ into account. This mimicking task may be too hard to accomplish if $I_S$ is not a.e.-equivalent to a Riemann-integrable function. Therefore, I think in general there is a counterexample -- but it would be an instructive surprise otherwise!
