Estimation of multivariate probability density function (PDF) of soil properties based on limited site investigation data plays an important role in data driven-based geotechnical engineering analysis and design. However, this estimation is often intractable because the site-specific investigation data generally are multivariate, uncertain and unique, sparse, incomplete, sometimes corrupt and spatially variable. The estimation is more challenging when the uncertainties in model parameters of the distribution (i.e., mean, standard deviations and scale of fluctuation of soil properties) are taken into consideration. To address these difficulties, this study proposed a mixed sampling technique by combing the Gibbs sampler and transitional Markov Chain Monte Carlo sampler to effectively estimate the multivariate probability density function of soil properties. The proposed method is demonstrated through a virtual site example while the effect of statistical uncertainty in the scale of fluctuation on reliability analyses is illustrated in an infinite slope example. It is shown that the proposed method can not only handle the complex features of site investigation data but also effectively quantify the statistical uncertainties of the mean values, standard deviations and scale of fluctuation of the investigated soil properties. Ignoring the statistical uncertainties of scale of fluctuation can affect the estimation of mean and standard deviations of the soil parameters. When considering the statistical uncertainty of scale of fluctuation, the 95% confidence intervals for the profiles of the soil properties can rationally cover the actual soil data at the unobserved points. Ignoring the statistical uncertainty of scale of fluctuation can lead to the underestimation of the probability of failure of the infinite slope.