Neural decoding for brain-computer interfaces (BCI) traditionally relies on spike sorting, which often suffers from inaccuracies and loss of information. Leveraging advances in high-density probes and computational methods, we propose a spike sorting-free decoding method using a mixture of Gaussians to model the distribution of spike features, encoding uncertainty without clustering spikes. We allow the mixing proportion of the MoG to change over time in response to the behavior and develop variational inference methods to fit the resulting model and to perform decoding. We benchmark our method with an extensive suite of recordings from different animals and probe geometries, demonstrating that our proposed decoder can consistently outperform current methods based on thresholding (i.e. multi-unit activity) and spike sorting.