This work leverages recent advances in probabilistic machine learning to discover governing equations expressed by parametric linear operators. Such equations involve, but are not limited to, ordinary and partial differential, integro-differential, and fractional order operators. Here, Gaussian process priors are modified according to the particular form of such operators and are employed to infer parameters of the linear equations from scarce and possibly noisy observations. Such observations may come from experiments or black-box computer simulations, as demonstrated in several synthetic examples and a realistic application in functional genomics. Employ probabilistic machine learning to discover governing equations expressed by parametric linear operators.Proper placement of Gaussian process priors allows one to efficiently infer model parameters via maximum likelihood estimation.A general treatment of inverse problems governed by linear operators, leading to model discovery from just a handful of noisy measurements.