We introduce a novel procedure that, given sparse data generated from a stationary deterministic nonlinear dynamical system, can characterize specific local and/or global dynamic behavior with rigorous probability guarantees. More precisely, the sparse data is used to construct a statistical surrogate model based on a Gaussian process (GP). The dynamics of the surrogate model is interrogated using combinatorial methods and characterized using algebraic topological invariants (Conley index). The GP predictive distribution provides a lower bound on the confidence that these topological invariants, and hence the characterized dynamics, apply to the unknown dynamical system. The proposed method is applied to a simple one-dimensional system to capture the existence of fixed points, periodic orbits, connecting orbits, bistability, and chaotic dynamics.