The prediction of chemical reactivity and thus the design of new reaction systems are the key challenges in chemistry. Here, we develop an unbiased general-purpose reaction sampling method, the stochastic surface walking based reaction sampling (SSW-RS) method, and show that the new method is a promising solution for reactivity prediction of complex reaction systems. The SSW-RS method is capable of sampling both the configuration space of the reactant and the reaction space of pathways, owing to the combination of two recently developed theoretical methods, namely, the stochastic surface walking (SSW) method for potential energy surface (PES) exploration and the double-ended surface walking (DESW) method for building pathways. By integrating with first principles calculations, we show that the SSW-RS method can be applied to investigate the kinetics of complex organic reactions featuring many possible reaction channels and complex hydrogen-bonding networks, as demonstrated here using two examples, epoxypropane hydrolysis in aqueous solution and β-D-glucopyranose decomposition. Our results show that simultaneous sampling of the soft hydrogen-bonding conformations and the chemical reactions involving hard bond making/breaking can be achieved in the SSW-RS simulation, and the mechanism and kinetics can be predicted without a priori information on the system. Unexpected new chemistry for these reactions is revealed and discussed. In particular, despite many possible pathways for β-D-glucopyranose decomposition, the SSW-RS shows that only β-D-glucose and levoglucosan are kinetically preferred direct products and the 5- or 7-member ring products should be secondary products derived from β-D-glucose or levoglucosan. As a general tool for reactivity prediction, the SSW-RS opens a new route for the design of rational reactions.