This paper focuses on variational inference with intractable likelihood functions that can be unbiasedly estimated. A flexible variational approximation based on Gaussian mixtures is developed, by adopting the mixture population Monte Carlo (MPMC) algorithm in \cite{cappe2008adaptive}. MPMC updates iteratively the parameters of mixture distributions with importance sampling computations, instead of the complicated gradient estimation of the optimization objective in usual variational Bayes. Noticing that MPMC uses a fixed number of mixture components, which is difficult to predict for real applications, we further propose an automatic component--updating procedure to derive an appropriate number of components. The derived adaptive MPMC algorithm is capable of finding good approximations of the multi-modal posterior distributions even with a standard Gaussian as the initial distribution, as demonstrated in our numerical experiments.