Promising quantum technologies such as superconducting qubits suffer from destructive noise with a characteristic 1/f spectrum. It has been shown that this is due to the presence of strongly coupled two-level defects. Their exact distribution would usually not be known. The resulting non-Markovian dynamics of the qubit is hence difficult to model -- a signifacant hurdle in qubit control and correction. We develop an algorithm to find a Lindblad-type model given a set of qubit population measurements. It allows for a number of coherently coupled fictitious two-level systems standing in for the defects. Then it stochastically explores the space of possible models. It aims to keep the Hilbert-space dimension as low as possible to simplify evaluation, while identifying the systems and processes required to achieve the observed dynamics. Our goal is to create a tool that can produce effective models for qubits in noisy environments, as well as investigate more general questions on the reducibility of non-Markovian models.