Related to #3. If we add size(), indexing and iteration, the discrete_distribution class becomes to look much like std::vector. I'm tempted to add vector-like functions like push_back() and resize(). It is normal to incrementally add reactions in the setup phase of a Gillespie simulation, so the vector-like functions will be actually useful.
std::discrete_distribution has probabilities() member function to get normalized probability distribution of the events. Maybe I should add this feature for compatibility.