How to create a discrete representation of a distribution function? I have a polynomial function representing a distribution function: $x^4 -6x^3 + 9x^2 ; 0 \le x \le 3$ but I need to get a discrete set of $ n $ data points representative of this distribution. I've seen there's many existing code libraries to do this for more standard distributions like binomial and normal, but is there a way I can do this for a given polynomial distribution function?
I had the idea of trying to divide the region into rectangles of equal area, then choosing one value from within each rectangle, but I would not know how to go about solving for each rectangle.

There are a lot of ways to do that. Here is one example that is easy to understand but is not efficient in Mathematica code for n=1000,
testxall = {}; testyall = {}; n = 1;
 While[n <= 1000,
  testx = 3 RandomReal[];
  testy = 5.0625 RandomReal[];
  If[testy <= testx^4 - 6 testx^3 + 9 testx^2, 
   testxall = AppendTo[testxall, testx]; 
   testyall = AppendTo[testyall, testy]; n++]]
lp1 = ListPlot[Transpose[{testxall, testyall}]];
pl1 = Plot[{x^4 - 6 x^3 + 9 x^2, 5.0625}, {x, 0, 3}];
Show[pl1, lp1]


The answer you asked for is the testxall array. What this does is generate rectangular {x,y} coordinates, and only count those that are less than or equal to the function. Note that the density is more compact in the x direction than in the y direction, but that is irrelevant.
Another, more efficient, way is to randomly generate probabilities, and that would be rather more commonly used.
