In this paper we introduce a new class of multivariate unimodal distributions, motivated by Khintchine's representation. We start by proposing a univariate model, whose support covers all the unimodal distributions on the real line. The proposed class of unimodal distributions can be naturally extended to higher dimensions, by using the multivariate Gaussian copula. Under both univariate and multivariate settings, we provide MCMC algorithms to perform inference about the model parameters and predictive densities. The methodology is illustrated with univariate and bivariate examples, and with variables taken from a real data-set. example, [29].The traditional approach is to model the data via a mixture model, see [31] for finite mixture distributions, and, more recently, [13]. The idea, which carries through to infinite mixture models, is that each component in the mixture is represented by a parametric density function, typically from the same family, which is usually the normal distribution. This assumption, often made for convenience, supposes each cluster or group can be adequately modeled via a normal distribution. For if a group requires two such normals, then it is deemed there are in fact two groups. Yet two normals could be needed even if there is one group and it happens to be skewed, for example.On the other hand, if we start by thinking about what a cluster could be represented by in terms of probability density functions, then a unimodal density is the most obvious choice. Quite simply, with lack of further knowledge, i.e. only observing the data, a bimodal density would indicate two clusters. This was the motivation behind [25], which relies heavily on the unimodal density models being defined on the real line, for which there are no obvious extensions to the multivariate setting.There are representations of unimodal density functions on the real line (e.g. [18] and [10]) and it is not difficult to model such a density adequately. See, for example,[3], [23] and [25]. Clearly, the aim is to provide large classes of unimodal densities and hence infinite dimensional models are often used.However, while there is an abundance of literature on unimodal density function modeling on the real line, there is a noticeable lack of work in two and higher dimensions. The reasons are quite straightforward; first there is a lack of a representation for unimodal distributions outside of the real line and, second, the current approaches to modeling unimodal densities on the real line do not naturally extend to higher dimensions.The aim in this paper is to introduce and demonstrate a class of unimodal densities for a multivariate setting. While marginal density functions will be modeled using the [18] representation on the real line, the dependence structure will be developed using a Gaussian copula model. To elaborate, using an alternative form of representation of