It has been recognized that many complex dynamical systems in the real world require a description in terms of multiplex networks, where a set of common, mutually connected nodes belong to distinct network layers and play a different role in each layer. In spite of recent progress towards data based inference of single-layer networks, to reconstruct complex systems with a multiplex structure remains largely open. We articulate a mean-field based maximum likelihood estimation framework to solve this outstanding and challenging problem. We demonstrate the power of the reconstruction framework and characterize its performance using binary time series from a class of prototypical duplex network systems that host two distinct types of spreading dynamics. In addition to validating the framework using synthetic and real-world multiplex networks, we carry out a detailed analysis to elucidate the impacts of structural and dynamical parameters as well as noise on the reconstruction accuracy and robustness.