The effectiveness of non-parametric methods for regression comes at the price of high computational complexity. In fact, these methods scale as O(N3), where N is the number of available data points. One possible option to address this issue consists in introducing a set of fictitious (pseudo-) inputs of size M ≪ N such that the computational effort is reduced to O(M2N). To estimate hyper-parameters and pseudo-inputs, a non-convex optimization problem needs to be solved. As opposed to the conventional gradient-based approach used in the literature, this paper proposes a stochastic scheme leveraging Markov Chain Monte Carlo methods. Numerical comparisons show that the latter returns a more efficient set of pseudo-inputs, leading to a superior performance in terms of mean squared error.