Modeling and forecasting the behavior of nonlinear systems remains an active area of research. Yet, the governing laws of these systems are often unknown or prohibitively expensive to simulate. With the proliferation of measurement sensors and the advent of efficient computational hardware, our modern world is abundant with multi-fidelity data that enables us to construct data-driven models without relying solely on first principles. Transfer operator theory describes linear representations of nonlinear dynamical systems and provides a pathway to produce analytically informative algorithms. Recently, transfer operator theory has been combined with another concept that is popular in data science: reproducing kernel Hilbert spaces (RKHS). Both transfer operator and RKHS techniques frame their corresponding problems in a functional analytic setting. Therefore, it stands to reason that by combining methods from both of these two fields, we can design algorithms that provide compact and interpretable models that can be readily extended to control of nonlinear dynamical systems. We follow this thread into Gaussian process (GP) methods, and illustrate how the Bayesian framework can alleviate two pervasive problems with kernel-based Koopman algorithms. The first being sparsity: most kernel methods do not scale well and require an approximation to become practical. Using numerical experiments on simulated nonlinear and stochastic dynamical systems we show that not only can the computational demands of kernel-based algorithms be reduced, but also demonstrate improved resilience against sensor noise. The second problem involves hyperparameter optimization and dictionary learning to fine-tune the model. The main contribution of the work is the unification of GP regression and dynamic mode decomposition (DMD). We propose a Bayesian DMD algorithm that treats the Perron-Frobenius operator as a random variable, and enables propagating uncertainties in the eigenfunctions, and provides a criterion for reprojections.