Average over all possible states
We assume that the system randomly samples states and that the average behavior coincides with the average over all possible states. We don't know the exact trajectory between the states and we don't care. We just hope that the system rattles around enough that the average over a long time or many particles will be similar to the average over all possible states.
The key to statistical mechanics is assuming that the state of the system is randomized over time and it can visit all possible states, subject to whatever constraints we have. We then can calculate macroscopic quantities (pressure, density, magnetization, temperature, even things like color—see the blackbody spectrum, etc.) by averaging over all possible states:
$$\left< A \right> = \sum_{\mathrm{all~possible~states}~i}{A(\mathbf{X}_i) P(\mathbf{X}_i)},$$
where $A$ is the macroscopic quanity, $\mathbf{X}_i$ is a state of the system (usually specifying the 3D positions and momenta for each particle), and $P(\mathbf{X}_i)$ is the probability of the state $\mathbf{X}$.
Here, I'm assuming the states are discrete (countable), if they are continuous, then the sum becomes an integral. For particles, this is usually an integral over momentum ($p$) and position ($q$) states. If we work with systems in contact with an environment of fixed temperature $T$, then the probability is $\exp (-\frac{E}{kT})$, where the $E$ is the energy of the state and $k$ is Boltzmann's constant. This relation comes from assuming that energy is randomly distributed between the system and environment (in other words, maximization of the total entropy). Let's also assume we have $N$ particles and that these particles are conserved. Under these conditions, we can calculate the average of some quantity $A$ by:
$$ \left< A \right> = \frac{\int \mathrm{d}^{3N}\mathbf{q} \int \mathrm{d}^{3N}\mathbf{p} \, A(\mathbf{q},\mathbf{p}) \exp \left[ -E(\mathbf{q},\mathbf{p})/(kT) \right] }{ \int \mathrm{d}^{3N}\mathbf{q} \int \mathrm{d}^{3N}\mathbf{p} \, \exp \left[ -E(\mathbf{q},\mathbf{p})/(kT) \right] }$$
To make it simple, let's say we have only one particle moving in one dimension that bounces elastically off the walls of the 1D box of length L and has only kinetic energy $E = \frac{1}{2}mv^2 = p^2/(2m)$. We want to calculate the mean energy:
$$ \left< E \right> = \frac{\int_0^L \mathrm{d}q \int_{-\infty}^\infty \mathrm{d}p \, p^2/(2m) \exp \left[ -p^2/(2mkT) \right] }{ \int_0^L \mathrm{d}q \int_{-\infty}^{\infty} \mathrm{d}p \, \exp \left[ -p^2/(2mkT) \right] }$$
The $q$ integrals give $L$ in both numerator and denominator, so they cancel.
I'll use sympy to do the $p$ integrals:
from sympy import *
kT, m = symbols('kT m', real=True, positive=True)
p = symbols('p', real=True)
E = p**2/(2*m)
numer = integrate(E*exp(-E/kT), (p, -oo, oo))
denom = integrate(exp(-E/kT), (p, -oo, oo))
mean_energy = simplify(numer/denom)
print(mean_energy)
So the average energy of this particle is $\left<E\right> = kT/2$. From this you can get a typical speed (root-mean-square speed). Since $\frac{1}{2}m\left<v^2\right> = \left<E\right>$, so $v_\mathrm{RMS} = \sqrt{kT/m}$
The point is, we average over all possible momentum states of the system (from $-\infty$ to $\infty$) and can use that to determine what the average energy and typical velocity are. We don't need to know the detailed trajectory of the particle in time. We don't need to think about how the particle bounces off the walls or how it exchanges energy with the environment through the walls. We just assume that on long times it effectively randomly samples states and that the average behavior will be the average over all states it can access.
So the field of statistical mechanics is about clever ways to average over states and come up with probability distributions for different quantities in model systems.
