Ok, I didn't expect that to work. I tried some other distributions myself and didn't find any other ones that worked even when while attempting some really elementary ones (like a 2D normal distribution). It's good to know polynomials work, but the fact remains that you basically have to be creative to get it to work and there are no guarantees. It's not like there's a robust MCMC sampler that just works for arbitrary PDFs. Something like that is quite high on my wish list, frankly speaking, and I'm a bit surprised that WL seems to be behind the curve with that, given how many applications there are for functionality like that.

Many thanks, Sjoerd and Marco. So it seems, as Sjoerd suggested, Wolfram hasn't implemented multivariate random number generation for user-defined distributions. It's certainly possible to use conditional probability for simulation, but let's hope Wolfram provide handy, automatic functions like RandomVariate for this purpose very soon.

I'm guessing that this is because sampling from arbitrary multi-dimensional distributions is generally not a simple problem which they haven't got around to implement yet. There are some trivial exceptions like ProbabilityDistribution[1, {x, 0, 1}, {y, 0, 1}], but in general it seems like you can only use build-in multi-dimensional distributions with RandomVariate.

In your case, you could try to sample the unit disk by transforming the coordinates to r and theta. This will factorize the distribution so you can draw from two 1D distributions, which should be easy to do.