Geotechnical data are typically Multivariate, Uncertain, and Irregular (MUI), so the probability distribution of geotechnical data is Multivariate, Asymmetric, and Multimodal (MAM). Probability Density Function (PDF) modelling and Credible Region (CR) construction are two key issues for a MAM distribution. There are two fundamental difficulties in characterizing a MAM distribution. The first is on joint PDF modelling as many traditional approaches collapse for a MAM distribution. Copula theory has attracted special attention for this purpose but very few works attempted to tackle the critical problem of probabilistic prediction on target variables using available information of remaining variables based on the copula-based joint PDF. The second is on CR construction of a MAM distribution as it cannot find a unique CR of a MAM distribution given an exceedance probability only. There is still a lack of a unified approach for CR construction for a MAM distribution of geotechnical data. Aiming to resolve these two fundamental difficulties, we propose the BAyeSIan Copula-based Highest density region/contour (BASIC-H) for providing a systematic framework of PDF modelling and CR construction of a MAM distribution. This framework contains Stage-PDF and Stage-CR. Stage-PDF fuses the copula theory and Bayesian inference to develop optimal, robust, and hyper-robust predictions on the posterior distribution and posterior predictive distribution. Stage-CR adopts the constraint for the CR that the probability density of every point inside the CR is at least as large as the probability density of any point outside, which is the same as the idea of the HDR (Highest Density Region). The Monte Carlo Simulation (MCS), based on the developed optimal, robust, and hyper-robust posterior distributions and posterior predictive distributions, is performed for estimation of the probability density boundary, which is a key parameter for constructing the HDR. Examples using simulated data and Quaternary clay data are presented to illustrate the capabilities of the BASIC-H in PDF modelling and CR construction of MAM distributions of geotechnical data.