What is the best way of getting random numbers in NumPy?

Your approach is fine. An alternative is to use the function numpy.random.uniform():

>>> numpy.random.uniform(-1, 1, size=10)array([-0.92592953, -0.6045348 , -0.52860837,  0.00321798,  0.16050848,       -0.50421058,  0.06754615,  0.46329675, -0.40952318,  0.49804386])

Regarding the probability for the extremes: If it would be idealised, continuous random numbers, the probability to get one of the extremes would be 0. Since floating point numbers are a discretisation of the continuous real numbers, in realitiy there is some positive probability to get some of the extremes. This is some form of discretisation error, and it is almost certain that this error will be dwarved by other errors in your simulation. Stop worrying!

Note that numpy.random.rand allows to generate multiple samples from a uniform distribution at one call:

>>> np.random.rand(5)array([ 0.69093485,  0.24590705,  0.02013208,  0.06921124,  0.73329277])

It also allows to generate samples in a given shape:

>>> np.random.rand(3,2)array([[ 0.14022471,  0.96360618],        [ 0.37601032,  0.25528411],        [ 0.49313049,  0.94909878]])

As You said, uniformly distributed random numbers between [-1, 1) can be generated with:

>>> 2 * np.random.rand(5) - 1array([ 0.86704088, -0.65406928, -0.02814943,  0.74080741, -0.14416581])

From the documentation for numpy.random.random_sample:

Results are from the “continuous uniform” distribution over the stated interval. To sample Unif[A, b), b > a multiply the output of random_sample by (b-a) and add a:

 (b - a) * random_sample() + a

Per Sven Marnach's answer, the documentation probably needs updating to reference numpy.random.uniform.