The most direct way to generate velocities for a set of N particles is to use the distribution function, which in principle provides the relative probability of a star to have a certain position and velocity.

There is a simple rejection method which can then be used to obtain a sampling consistent with the original distribution function. This method, known as the von Neumann rejection method, has been described so often, that here, we only provide a quick description of its principle using a simple illustration. Let’s imagine we wish to sample a positive distribution which depends on one variable (i.e. velocity) described by . The process is then : 0- First, define the range of allowed v’s (vmin - vmax). and then derive the maximum Pmax value for in this range.
 1 For one particle, select a velocity so that .
 2 For that particle, select P_i such that
 3 the rejection technique then tells you to reject that draw if , but accept otherwise.
 4 if it is accepted, proceed with the next particle, and otherwise go back to step 1 for that particle. This guarantees that the obtained sampled set of particles follows the functional form you originally choose.

The von Neumann rejection method is efficient only if is as close to Pmax as possible for all v’s (constant function). If this is not the case (e.g. strong peak at one v value, but very low otherwise), the number of rejection points will be high and the technique very inneficient. It is therefore advised to allow some change of variable as to uniformize the $P(v) function before the rejection technique is applied.

A description of this technique is available in the Numerical Recipes (Press et al. 1993), and has been used by many to generate a set of N particles consistently with some a priori distribution functions. One of the best know examples is the Plummer sphere as mentioned by Aarseth, Henon, Wielen 1974 (A&A 37, 183), and has been applied for axisymmetric multicomponents systems by e.g. Kuijken & Dubinski 1994 (MNRAS 2777, 1341).