Let me first remark on the nature of your proposal. It's a good one, but we must understand that any realization of this random process will rarely have exactly the intended species proportion in each polygon: it will only attain those proportions on average (as statistical expectations). The actual proportions will vary from their expectations from one polygon to another.
If you want to achieve exactly the intended proportions, then you are need a random sample that is stratified by the polygons. For approaches to that (harder) problem, please see the thread at Stratified Random Sampling with large dataset in ArcGIS?.
Data preparation
A convenient way to represent the polygon-specific species distribution would be to create four grids, one for each species. Each cell gives the chance of that species being there. The grids may contain only non-negative values (and no NoData values inside any polygons) and they must sum to 1 at each polygon cell.
Convert these grids to cumulative distribution grids by taking their cumulative sum. The order doesn't really matter, so to illustrate let's suppose these species-probability grids are called "P1", "P2", "P3", and "P4". Then one way to create the cumulative sum grids is to compute
C0 = 0
C1 = "P1"
C2 = "C1" + "P2"
C3 = "C2" + "P3"
Simulation
With this preprocessing accomplished, the actual simulation is remarkably simple and efficient. As proposed in the question, create a grid of uniform random values in the interval [0,1]. Let's call this "U". The key to this solution is to use the Less Than Frequency
operation. When applied to "U" and the set {"C0", "C1", "C2", "C3"}, it will produce a grid with values in the set {0, 1, 2, 3, 4}--but it should have no zero values, because every value in "U" will exceed zero (the values in "C0"). These numbers encode the randomly selected species: 1 for the species represented in "P1", 2 for the species in "P2", and so on. That is the desired reclassification.
Analysis of the algorithm
To see why this method works, consider the chance that the Less Than Frequency
grid equals 2 (say). This will happen when exactly two of {"C0", ... "C3"} are less than "U". Because the values of the "C" grids are always in ascending order, this will happen precisely when "U" lies between "C1" and "C2". The chance of this occurring is "C2" - "C1" = "P2", as intended. Because the "C" grids are actually grids, and not just numbers, their values may vary from one cell to the next: that is how they accomplish this cell-specific reclassification operation.