Joint factor analysis (JFA) has been successfully applied to speaker verification tasks to tackle speaker and session variability. In the sense of Bayesian statistics, it is beneficial to take account of the uncertainties in JFA to better characterize its speaker enrollment and verification processes, e.g. representing target speaker model by posteriori distribution of latent speaker factors and evaluating model likelihood by integrating over all latent factors. However, in a JFA model which has a large number of latent factors, it is computationally demanding to carry out these things in their exact form. In this paper, an alternative approach based on variational Bayesian is developed to explore uncertainties in JFA in an approximate yet efficient way. In this method, fully correlated posteriori distribution is approximated by a variational distribution of factorial form to facilitate inference; and a tight lower bound on model likelihood is derived. Experimental results on the 10sec4w-10sec4w task of the 2006 NIST SRE show that variational Bayesian JFA could obtain better performance than JFA using point estimate.