<p>This study proposes a two-stage data-driven approach based on the Gaussian process (GP) for predicting the remaining useful life (RUL) of engineering systems using both historical data and current observational data. This method effectively fuses multi-source degradation data, efficiently estimates GP hyperparameters, and further predicts the RUL of the systems. In the first stage, the transitional Markov Chain Monte Carlo method is applied within the Bayesian modeling framework to obtain the posterior distribution of the GP parameters. An explicit mean function is incorporated into the GP, where its parameters and kernel hyperparameters are consolidated into a unified setting for joint estimation from historical datasets. In the second stage, historical information is utilized as the prior information for GP parameter estimation, and a Monte Carlo simulation strategy is used in conjunction with the updated parameter distribution. This parameter distribution is further used for efficient calculation of the probability distribution of the RUL. The effectiveness of the proposed method is validated through two real examples involving battery capacity degradation and crack propagation. The results demonstrate that the proposed two-stage GP method provides an effective and accurate solution for the posterior distribution of hyperparameters and RUL prediction. While traditional point estimation methods often yield similar prediction results, they typically fail to provide comprehensive uncertainty quantification information.</p>