Characterizing probability distributions of geotechnical data plays an important role in data-centric geotechnics. On the one hand, geotechnical data are Multivariate, Uncertain, and Irregular (MUI), where the irregular characteristic implies that asymmetry and/or multimodality are often observed in the histograms of geotechnical data, so the corresponding probability distribution is Multivariate, Asymmetric, and Multimodal (MAM). On the other hand, many geotechnical datasets are unavoidably subjected to the issue of modelling and prediction stages missing values (called “dual-stage missing values”), so characterising the MAM distribution of geotechnical data with dual-stage missing values becomes an essential task. There are three fundamental difficulties for this purpose. The first is on joint Probability Density Function (PDF) modelling for a MAM distribution given data with modelling stage missing values. Many traditional and advanced approaches collapse in the presence of MAM distributions and modelling stages missing values, respectively. The second is on joint PDF prediction for a MAM distribution given data with prediction stage missing values. The third is on Credible Region (CR) construction of a MAM distribution as there is no unique CR of a MAM distribution given an exceedance probability only. We propose the three-stage BAyeSIan Copula-based Highest density region/contour (BASIC-H). Stage-1 constructs the posterior distribution of data with modelling stage missing values based on Copula theory and Bayesian inference. Stage-2 derives the posterior predictive distribution of data with prediction stage missing values based on marginalisation and conditionalisation of the posterior distribution. Stage-3 constructs the CRs for the posterior and predictive distributions adopting the reasonable constraint imposed by the Highest Density Region (HDR). Examples using simulated data, CLAY/10/7490 and CLAY/5/345 are presented to illustrate the capability of the proposed BASIC-H.