Accurate numerical modeling of laser shock processing, a typical complex physical process, is very difficult because several input parameters in the model are uncertain in a range. And numerical simulation of this high dynamic process is very computational expensive. The Bayesian Gaussian process method dealing with multivariate output is introduced to overcome these difficulties by constructing a predictive model. Experiments are performed to collect the physical data of shock indentation profiles by varying laser power densities and spot sizes. A two-dimensional finite element model combined with an analytical shock pressure model is constructed to obtain the data from numerical simulation. By combining observations from experiments and numerical simulation of laser shock process, Bayesian inference for the Gaussian model is completed by sampling from the posterior distribution using Morkov chain Monte Carlo. Sensitivities of input parameters are analyzed by the hyperparameters of Gaussian process model to understand their relative importance. The calibration of uncertain parameters is provided with posterior distributions to obtain concentration of values. The constructed predictive model can be computed efficiently to provide an accurate prediction with uncertainty quantification for indentation profile by comparing with experimental data.