Methods from stochastic dynamical systems theory have been instrumental in understanding the behaviours of chemical reaction networks (CRNs) arising in natural systems. However, considerably less attention has been given to the inverse problem of synthesizing CRNs with a specified behaviour, which is important for the forward engineering of biological systems. Here, we present a method for generating discrete-state stochastic CRNs from functional specifications, which combines synthesis of reactions using satisfiability modulo theories and parameter optimization using Markov chain Monte Carlo. First, we identify candidate CRNs that have the possibility to produce correct computations for a given finite set of inputs. We then optimize the parameters of each CRN, using a combination of stochastic search techniques applied to the chemical master equation, to improve the probability of correct behaviour and rule out spurious solutions. In addition, we use techniques from continuous-time Markov chain theory to analyse the expected termination time for each CRN. We illustrate our approach by synthesizing CRNs for probabilistically computing majority, maximum and division, producing both known and previously unknown networks, including a novel CRN for probabilistically computing the maximum of two species. In future, synthesis techniques such as these could be used to automate the design of engineered biological circuits and chemical systems.