This question seems to be rather complicated. I would guess the anwer depends for example on how two particles interact ("collide"), the cross section etc. I will have to ignore such subtleties here. My answer is more an educated guess than a qualified derivation, so take it with a large grain of salt. I would actually also be interested to know if there is a more qualified approach.
Let us first recall some things about statistical mechanics:
Assuming a canonical ensemble the probability to find the system (i.e. all particles) in the phase space volume $dpdq$ (where $q$ is the vector containing the coordinates of all particles, and $p$ the vector containing the momentum components of all particles, so they have a lot of components and the dimension of our phase space is very large) is proportional to $dqdp e^{-\beta H(q,p)}$ where $\beta=\frac{1}{k_b T}$ and $H$ is the Hamiltonian of the whole system. The integral over all of that is $Z$ (with some constants). If one wants the average of any observable $O(q,p)$ one integrates $dqdp e^{-\beta H(q,p)}O(q,p)$ and devides by $Z$ for normalization.
Let us now assume an ideal gas, i.e. $H(q,p)=\frac{p^2}{2m}$. This might or might not be a good assumption. How can it ever be valid if particles can collide, i.e. interact? Well maybe the interaction has a very short range. If we don't take the Hamiltonian of the ideal gas it should still be possible to go through the calculation but it will be a lot less clear. So I will just use the ideal gas as a hopefully instructive example.
The probability that particle $i$ collides with particle $j$ now depends on what we actually mean with "collide". One could imagine both particles to be something like little balls with a particular radius and if the positions of their centers are closer than twice that radius they collide. Because all particles are indistinguishable (or at least I also assume that) to get the average collision velocity per particle pair I only have to look at one pair. Instead of the "normal" canonical partition function (integral of $dqdp e^{-\beta H(q,p)}$ with some constants in front of it) we should only integrate over states where our collision condition from before is met. But that will turn out to have no effect at all because we assumed an ideal gas, i.e. $H$ has no position dependence and all position integrals will be trivial anyways. For the momentum integrals everything will cancel except the integral over the momentum $p_i$ of the $i$-th particle (attention: $p_i$ is a vector, not just the $i$-th component of $p$) and $p_j$. Here you will have the following:
$$ \frac{\int dp_idp_j (p_i-p_j)^2 e^{-\beta\frac{p_i^2+p_j^2}{2m}} }{\int dp_idp_j e^{-\beta\frac{p_i^2+p_j^2}{2m}}}$$
Here $p_i$ and $p_j$ both get integrated over the full $\mathbb{R}^3$. The integral in the denominator is just for normalization (it is the parti of the canonical partition function $Z$ that didn't cancel because the corresponding term in the numerator is the only one where the integral is different from the one in $Z$). It is just some Gaussian integral and is easy to calculate to be something like $\left(\frac{2m\pi}{\beta}\right)^3$. The integral in the numerator is also easy to calculate, the result should then be something like the average square of the "collision momentum" if all assumptions are true. I know that isn't yet what you wanted but you might want to calculate it or something like it to see if it makes sense at all.
To now identify the distribution for collision velocities we first remove the integral sign and the $(p_i-p_j)^2$ in the numerator (the denominator will stay exactly the same and we will ignore it for now). Then we do a change of variables for $p_i$: The new variable is $p_\Delta=p_i-p_j$, i.e. the relative momentum. This transforms the numerator into $dp_\Delta dp_j e^{-\beta\frac{(p_\Delta+p_j)^2+p_j^2}{2m}}=dp_\Delta dp_j e^{-\beta\frac{p_\Delta^2+2p_\Delta p_j+2p_j^2}{2m}}$. Now we take the $p_j$-Integral because we want only the relative momentum and have to average over the total momentum of one particle. Because of rotational invariance we can put the $z$-axis of $p_j$ in the direction of $p_\Delta$, so $p_\Delta^2+2p_\Delta p_j+2p_j^2$ becomes $p_\Delta^2+2\vert p_\Delta\vert p_{j_z}+2p_{j_x}^2+2p_{j_y}^2+2p_{j_z}^2$ where $p_{j_x}$ is the $x$-component of $p_j$ and so on. The whole numerator now is
$$ dp_\Delta \int dp_{j_x}dp_{j_y}dp_{j_z} e^{-\beta\frac{p_\Delta^2+2\vert p_\Delta\vert p_{j_z}+2p_{j_x}^2+2p_{j_y}^2+2p_{j_z}^2}{2m}}=dp_\Delta \left(\frac{m \pi}{\beta} \right) ^{\frac{3}{2}} e^{-\beta\frac{p_\Delta^2}{4m}} $$
If one now puts back the denominator and and compares with the probability $dp_\Delta f(p_\Delta)$ for the distribution $f$ one gets the normalized distribution
$$ f(p_\Delta)=\left(\frac{\beta}{4m \pi} \right) ^{\frac{3}{2}} e^{-\beta\frac{p_\Delta^2}{4m}} $$
This should now be our result for the "ideal gas" (plus possible collisions, so not completely ideal) with several assumptions which seem more or less plausible.
I am unsure whether this actually is the thing you are looking for but maybe it at least gets you one step closer to your goal. Remember that the integrals will be significantly less easy to solve once you introduce interactions between the particles or an external potential. You should still be able to write down the integrals (which become a lot more cumbersome because less stuff cancels) and take numerical solutions.
I hope my explanation, if you can call it that, was sufficient and you can replicate the calculation on your own, maybe for a more general Hamiltonian.
