Sampling from a multivariate probability density function in python
There a few different paths one can follow here.
(1) If P(x,y,z) factors as P(x,y,z) = P(x) P(y) P(z) (i.e., x, y, and z are independent) then you can sample each one separately.
(2) If P(x,y,z) has a more general factorization, you can reduce the number of variables that need to be sampled to whatever's conditional on the others. E.g. if P(x,y,z) = P(z|x, y) P(y | x) P(x), then you can sample x, y given x, and z given y and x in turn.
(3) For some particular distributions, there are known ways to sample. E.g. for multivariate Gaussian, if x is a sample from a mean 0 and identity covariance Gaussian (i.e. just sample each x_i as N(0, 1)), then y = L x + m has mean m and covariance S = L L' where L is the lower-triangular Cholesky decomposition of S, which must be positive definite.
(4) For many multivariate distributions, none of the above apply, and a more complicated scheme such as Markov chain Monte Carlo is applied. 
Maybe if you say more about the problem, more specific advice can be given.
I have a multivariate probability density function P(x,y,z), and I want to sample from it. Normally, I would use numpy.random.choice() for this sort of task, but this function only works for 1-dimensional probability densities. Is there an equivalent function for multivariate PDFs?