This paper presents Monte Carlo simulations of tunneling recombination in random distributions of defects. Simulations are carried out for four common luminescence phenomena in solids exhibiting tunneling recombination, namely continuous wave infrared stimulated luminescence (CW-IRSL), thermoluminescence (TL), isothermal thermoluminescence (iso-TL) and linearly modulated infrared stimulated luminescence (LM-IRSL). Previous modeling work has shown that these phenomena can be described by the same partial differential equation, which must be integrated numerically over two variables, the elapsed time and the donor-acceptor distance. We here present a simple and fast Monte Carlo approach which can be applied to these four phenomena, and which reproduces the solution of the partial differential equation, without the need for numerical integrations. We show that the method is also applicable to cases of truncated distributions of nearest neighbor distances, which characterize samples that underwent multiple optical or thermal pretreatments. The accuracy and precision of the Monte Carlo method are tested by comparing with experimental data from several feldspar samples.