A probabilistic data-driven approach that models the filtered reaction rate in large-eddy simulation (LES) is investigated. We propose a novel framework that incorporates a conditional generative adversarial network and a Gaussian mixture model to take into account the statistical fluctuations that are present in LES of turbulent reacting flows due to non-resolved subgrid structures, which cannot be predicted by purely deterministic models and machine learning algorithms. The data from a direct numerical simulation of turbulent premixed combustion are spatially filtered using a wide range of filter widths and employed for the training. We extract physically relevant parameters from the database and reduce the input features to the network to the most influential ones based on the result of feature importance analysis. The trained model is then tested on unseen timesteps and untrained LES filter widths, where it is able to accurately predict the distribution of the filtered reaction rate.