Basically, make a cumulative probability distribution (CDF) array. Basically, the value of the CDF for a given index is equal to the sum of all values in P equal to or less than that index. Then you generate a random number between 0 and 1 and do a binary search (or linear search if you want). Here's some simple code for it.
from bisect import bisect
from random import random
P = [0.10,0.25,0.60,0.05]
cdf = [P[0]]
for i in xrange(1, len(P)):
cdf.append(cdf[-1] + P[i])
random_ind = bisect(cdf,random())
of course you can generate a bunch of random indices with something like
rs = [bisect(cdf, random()) for i in xrange(20)]
yielding
[2, 2, 3, 2, 2, 1, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2, 2, 2, 2]
(results will, and should vary). Of course, binary search is rather unnecessary for so few of possible indices, but definitely recommended for distributions with more possible indices.
random.choices
(note the 's' at the end) which allows submitting relative weights.random.choices
information into the top answer, since the interface is substantially the same.