You have to go back to see where your equation $(1)$ came from.
For one dimension the probability of particles having a velocity between $\vec v_{\rm x}$ and $\vec v_{\rm x}+ d\vec v_{\rm x}$ is given by
$$f(\vec v_{\rm x})\;d\vec v_{\rm x} =\sqrt{\frac{m}{ 2\pi kT}}\exp \left[\frac{-mv_{\rm x}^2}{2kT}\right]\,d\vec v_{\rm x}$$
The speed distribution is given by
$$f(v_{\rm x})\;dv_{\rm x} =2\,\sqrt{\frac{m}{ 2\pi kT}}\exp \left[\frac{-mv_{\rm x}^2}{2kT}\right]\,dv_{\rm x}=\sqrt{\frac{2m}{ \pi kT}}\exp \left[\frac{-mv_{\rm x}^2}{2kT}\right]\,dv_{\rm x}$$
the factor $2$ being there because the speed (magnitude) of $\vec v_{\rm x}$ is the same as that of $-\vec v_{\rm x}$.
This is in agreement with your equation $(3)$.
Equation $(1)$ came from the idea that in three dimensions there is no preferred direction so
$f(\vec v_{\rm x},\vec v_{\rm xy},\vec v_{\rm z})\,d\vec v_{\rm x}d\vec v_{\rm y}d\vec v_{\rm z} = f(\vec v_{\rm x})d\vec v_{\rm x}\,f(\vec v_{\rm y})d\vec v_{\rm y}\,f(\vec v_{\rm z})d\vec v_{\rm z}$
You now have to count all the speeds which are the same ie the magnitude of the velocity $v$ is the same where $v^2 = v^2_{\rm x}+v^2_{\rm y}+v^2_{\rm z}$.
In this three dimensional case the volume of a shell of radius $v$ and thickness $dv$ is $d\vec v_{\rm x}d\vec v_{\rm y}d\vec v_{\rm z}= 4 \pi v^2 dv$ is being considered which results in your equation $(1)$.
$$f(v) \,dv = \sqrt{\left(\frac{m}{2 \pi kT}\right)^3}\, 4\pi v^2 \exp \left[\frac{-mv^2}{2kT}\right] \,dv$$
The equivalent distribution for two dimensions with an area of a ring of radius $v$ and thickness $dv$ is $d\vec v_{\rm x}d\vec v_{\rm y}= 2 \pi v\, dv$ and $v^2 = v^2_{\rm x}+v^2_{\rm y}$ is
$$f(v) \,dv = \left(\frac{m}{2 \pi kT}\right)\, 2\pi v\, \exp \left[\frac{-mv^2}{2kT}\right] \,dv$$
Best Answer
According to your formulas, $$\langle v_x \rangle = \left( \frac{m}{2\pi k_B T} \right)^{3/2} \int_{-\infty}^\infty v_x e^{-\frac{m[(v_x - a)^2 + v_y^2 +v_z^2]}{2k_B T}} d v_x d v_y dv_z $$
Shifting $v_x$ by defining $v'_x = v_x - a$ we have
$$\langle v_x \rangle = \left(\frac{m}{2\pi k_B T}\right)^{3/2} \int_{-\infty}^\infty (v'_x +a) e^{-\frac{m[v_x'^2 + v_y^2 +v_z^2]}{2k_B T}} d v'_x d v_y dv_z = a \left( \left(\frac{m}{2\pi k_B T}\right)^{3/2} \int_{-\infty}^\infty e^{-\frac{m[v_x'^2 + v_y^2 +v_z^2]}{2k_B T}} d v'_x d v_y dv_z\right) = a$$ I used that the $v'_x e^{-(...)}$ integral vanishes because it is odd in $v'_x$, and the remaining $a e^{-(...)}$ integral simplifies with the prefactor to give simply $a$ as expected.
This result just follows from identifying this distribution with a gaussian distribution of mean $(a,0,0)$ and standard deviation $\sigma = \sqrt{\frac{k_B T}{m}}$