Mutational processes, such as the molecular effects of carcinogenic agents or defective DNA repair mechanisms, are known to produce different mutation types with characteristic frequency profiles, referred to as mutational signatures. Non-negative matrix factorization (NMF) has successfully been used to discover many mutational signatures, yielding novel insights into cancer etiology and targeted therapies. However, the NMF model is only a rough approximation to reality, and even small departures from this assumed model can have large negative effects on the accuracy and reliability of the results. We propose a new approach to mutational signatures analysis that improves robustness to misspecification by using a power posterior for a fully Bayesian NMF model, while employing a sparsity-inducing prior to automatically infer the number of active signatures. In extensive simulation studies, we find that our proposed approach recovers more true signatures with greater accuracy than current leading methods. On whole-genome sequencing data for six cancer types from the ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium, we find that our method is able to accurately recover more signatures than the current state-of-the-art.