Reliable estimation of equilibrium flow shear from reciprocating probe measurements is challenging since the quantity of interest corresponds to the second derivative of the observable plasma potential. In addition, a time series of the plasma potential obtained by plunging a probe is affected by both the probe head position and plasma fluctuations, complicating the estimation of equilibrium components and their errors. We tackle this problem by employing Gaussian process regression that is able to infer even the derivatives of a spatial or temporal profile in the form of a probability distribution function. The proposed inference framework is validated by using synthetic data generated by gyrofluid simulations. While the inference result based on a single plunge is unstable in certain spatial locations, we have obtained reasonable agreement between the inference result and the true flow shear profile by combining data sets taken from several plunges.