AbstractVariability, stochastic or otherwise, is a central feature of neural circuits. Yet the means by which variation and uncertainty are derived from noisy observations of neural activity is often unprincipled, with too much weight placed on numerical convenience at the cost of statistical rigour. For two-photon imaging data, composed of fundamentally probabilistic streams of photon detections, the problem is particularly acute. Here, we present a complete statistical pipeline for the inference and analysis of neural activity using Gaussian Process Regression, applied to two-photon recordings of light-driven activity in ex vivo mouse retina. We demonstrate the flexibility and extensibility of these models, considering cases with non-stationary statistics, driven by complex parametric stimuli, in signal discrimination, hierarchical clustering and inference tasks. Sparse approximation methods allow these models to be fitted rapidly, permitting them to actively guiding the design of light stimulation in the midst of ongoing two-photon experiments.