The multinomial logistic Gaussian process is a flexible non-parametric model for multi-class classification tasks. These tasks are often involved in solving a pattern recognition problem in real life. In such contexts, the multinomial logistic function (or softmax function) is usually assumed to be the likelihood function. But, exact inferences for this model have proved challenging problem because it requires high-dimensional integration. In this paper, we propose approximate variational Bayesian inference for the multinomial logistic Gaussian process model. First, we compute the second-order approximation for the logarithm of the logistic likelihood function using Taylor series expansion, and derive the posterior distributions of all hidden variables and model parameters using the variational Bayesian inference method. Second, we derive the predictive distribution of the latent classification variable corresponding to the relevant test data point using the characteristics of the Cauchy product for a standard Gaussian process using a learning model parameter. We conducted experiments to verify the effectiveness of the proposed model using a number of synthetic and real datasets. The results show that the proposed model has superior classification capability to existing methods.