Let $X$ be a metric space and $\mu$ be any probability measure on $X$. Suppose that $f:X\to X$ is a continuous $\mu$-measurable function and that $(S_n)_n$ is a sequence of subsets of the natural numbers $\mathbb{N}$ such that $\# S_n\to\infty,n\to\infty$ for the counting measure $\#$.
Define the map $\alpha_n(B) := \frac{1}{\# S_n}\sum_{l\in S_n}\mu(f^{-l}(B))$ for $B\subset X$. It is then not hard to see that $\alpha_n$ is a probability measure on $X$.
Suppose then that $(n_k)_k$ is a subsequence of indices such that the subsequence $(\alpha_{n_k})_k$ converges in weak* topology to a probability measure $\alpha$ over $X$, i.e. for any bounded continuous function $g:X\to\mathbb{R}$ we have
$$\lim_{k\to\infty}\left|\int_Xg(y)d\alpha_{n_k}(y) - \int_X g(y)d\alpha(y)\right| = 0$$
where
$$\int_Xg(y)d\alpha_{n_k}(y) = \frac{1}{\# S_{n_k}}\sum_{l\in S_{n_k}}\int_X(g\circ f^l)(y)d\mu(y)$$
as one can verify easily.
(Extra assumption:) Take it as granted that if $f:X\to X$ is a continuous map then its pushforward map $f_*\mu(B) := \mu(f^{-1}(B)), B\subset X$ is continuous in weak* topology. I.e. if $\mu_n\to\mu$ in weak* as $n\to\infty$ then $\lim_{n\to\infty}f_*(\mu_n) = f_*\mu$.
(Problem:) I want to show that the prior limit measure $\alpha$ is $f$-invariant, that for any $\alpha$ measurable subset $B\subset X$ we have $\alpha(B) = \alpha(f^{-1}(B))$.
(Solution attempt/thoughts:)
Let $g:X\to \mathbb{R}$ be a fixed bounded and continuous function. By continuity of $f_*$,
$$f_*\alpha = f_*\left(\lim_{k\to\infty}\alpha_{n_k}\right) = \lim_{k\to\infty}\frac{1}{\# S_{n_k}}\sum_{l\in S_{n_k}}\mu\circ f^{-l - 1}$$
Therefore I would want to conclude that the difference
$$I := \left|\frac{1}{\# S_{n_k}}\sum_{l\in S_{n_k}}\int_X(g\circ f^{l})(y) - (g\circ f^{l+1})(y)d\mu(y)\right|$$
can be made arbitrarily small after some sufficiently large $k$. But since nothing guarantees that $l\in S_n\implies l+1\in S_n$ I have no idea how to relate these two quantities of the integrand.
One could use triangle inequality and the fact that $\alpha_{n_k}$s are probability measures to show that
$$\left|\int_X(g\circ f^{l}) - (g\circ f^{l+1})(y)(y)d\mu(y)\right|\leq 2\sup_{y\in X}\left|g(y)\right|$$
for any $l\in S_n$. But without additional information you would get this same bound for every summand occurring in $I$ and so the averaging is lost.
What to do? I am not necessarily looking for a full-blown answer, but I'll appreciate any help/hints that you might be able to give. Thanks!
I feel there are some assumptions missing here. For one thing, if you take say $$X = \mathbb{R}, f(x) = -x$$ and all $S_n$ to be subsets of even number, then $$\alpha_{n_k} = \mu,$$ which converges to $\mu$, but $\mu$ is not necessarily $f$-invariant.
I will update this with a solution once a clarification is made.