A pseudo-marginal Markov chain Monte Carlo (PMCMC) method is proposed for nonnegative matrix factorization (NMF). The sampler jointly simulates the joint posterior distribution for the nonnegative matrices and the matrix dimensions which indicate the number of the nonnegative components in the NMF model. We show that the PMCMC sampler is a generalization of a version of the reversible jump Markov chain Monte Carlo. An illustrative synthetic data was used to demonstrate the ability of the proposed PMCMC sampler in inferring the nonnegative matrices and as well as the matrix dimensions. The proposed sampler was also applied to a nuclear magnetic resonance spectroscopy data to infer the number of nonnegative components.