The initial velocities are drawn from a Gaussian distribution with variance
$$\sigma_i^2=\frac{k_{\textrm{B}}T}{m_i},$$
where $k_{\textrm{B}}$ denotes Boltzmann's constant, $T$ is the temperature and $m_i$ is the mass of the $i^{\textrm{th}}$ particle.

Thus, the problem boils down to generate random numbers from a gaussian distribution using uniformly distributed random numbers. This is fortunately quite simple: the Wikipedia article <https://en.wikipedia.org/wiki/Box–Muller_transform> shows some very common algorithms how to transform uniform random numbers into gaussian random numbers.

Let's put everything together: the velocity of the $i^{\textrm{th}}$ particle is computed via
$$v_i=\sqrt{\frac{k_{\textrm{B}}T}{m_i}}\,\mathcal{N}(0,1),$$
where $\mathcal{N}(0,1)$ is a gaussian random number with variance 1 and mean 0.

In case of vectorial velocities, initialise all components of each particle's velocity vector with this scheme, and remember to use a new gaussian random number for each component.