Eliciting Dirichlet and Gaussian copula prior distributions for multinomial models

In this paper,we propose novelmethods of quantifying expert opinion about prior distributions formultinomial models. Two different multivariate priors are elicited using median andquartile assessments of themultinomial probabilities. First, we start by eliciting a univariate beta distribution for the probability of each category. Then we elicit the hyperparameters of the Dirichlet distribution, as a tractable conjugate prior, from those of the univariate betas through various forms of reconciliation using least-squares techniques. However, a multivariate copula function will give a more flexible correlation structure between multinomial parameters if it is used as their multivariate prior distribution. So, second, we use beta marginal distributions to construct a Gaussian copula as a multivariate normal distribution function that binds these marginals and expresses the dependence structure between them. The proposed method elicits a positive-definite correlation matrix of this Gaussian copula. The two proposed methods are designed to be used through interactive graphical software written in Java.


Introduction
Bayesian statistical methods provide a formal mechanism for taking into account prior knowledge. In many circumstances, prior knowledge is based on historical data that are only recorded in the form of the personal experience of experts. Then expert opinion must be quantified as a prior distribution if the information is to be used. This can be done through an elicitation process which is defined as asking the expert for information that is then used as a basis for encoding the parameters of a subjective prior distribution. The method of quantifying opinion must be designed for the sampling model of interest and the form of prior distribution that is to be used. One important sampling model that has attracted limited attention in the elicitation context is the multinomial model. This model is used in many scientific and industrial applications. For example, it is frequently applied to voting behavior in political science (e.g. Chamberlin and Featherston 1986), to the compositions of rocks in geology (e.g. Weltje 2002) and to patterns of consumer selection preferences in microeconomics (e.g. Manski 2007). The model is so common that good elicitation methods for quantifying opinion about its parameters are clearly needed. For example, in their recent applied research, Fu et al. (2012) and Conn et al. (2013) stated that an elicited informative prior for the multinomial model would be a possible alternative.
Here we are concerned with the most common case of the multinomial model, where each observation has the same probability of falling into any specified category and observations are independent of each other. Then observations can be modelled to follow a multinomial distribution with, say, probability p i that an observation falls in the ith category. To quantify an expert's opinions about the p i , we assume some distributional form to represent her opinion and uncertainty. Instead of directly asking the expert to assess the parame-ters of the distribution, we ask her to make assessments that determine appropriate values for the parameters of that distribution, which are referred to as its hyperparameters.
The Dirichlet distribution, which is the conjugate prior distribution for the parameters of a multinomial distribution, is one natural choice for representing an expert's opinion. A number of elicitation methods that follow this approach have been proposed in the literature. For example, Bunn (1978) encoded a Dirichlet prior from assessed beta marginal distributions through hypothetical future samples. Chaloner and Duncan (1987) based their elicitation method on assessing a sample size and modal values of the Dirichlet variates. Dorp and Mazzuchi (2004) introduced a numerical algorithm for encoding the Dirichlet parameters from a number of quantile assessments that is equal to the number of parameters. Recently, two elicitation methods have been proposed for the Dirichlet distribution by reconciling the assessed parameters of its univariate beta conditional distributions  and of its univariate marginal beta distributions (Zapata-Vázquez et al. 2014).
The first method proposed in the current paper elicits the hyperparameters of the Dirichlet distribution from those of its marginal beta distributions through forms of reconciliation that use least-squares techniques. The usage of least-squares techniques for reconciliation has been discussed in Albert et al. (2012). Our proposed method differs from those of Bunn (1978) and Chaloner and Duncan (1987) as we require only medians and quartiles to be assessed. We assess more quartiles than Dorp and Mazzuchi (2004) and then reconcile them to encode hyperparameters. Our method encodes a Dirichlet distribution using its beta marginal distributions rather than its conditional distributions as in Elfadaly and Garthwaite (2013). Our proposed reconciliation methods are different from those that have been independently developed in Zapata-Vázquez et al. (2014) for reconciling beta marginal distributions into a Dirichlet distribution for a set of proportions.
However, while the standard Dirichlet distribution offers tractability and mathematical simplicity, it has been criticized as too inflexible to represent a broad range of prior information about the parameters of multinomial models (Aitchison 1986;O'Hagan and Forster 2004). Specifically, the Dirichlet variates are always negatively correlated and the Dirichlet distribution contains the same number of parameters as the number of categories. With k categories, k −1 parameters are required to describe the proportion in each category, leaving only one parameter to describe their dependence structure. This may sometimes be insufficient to give a reasonable representation of an expert's opinions.
In a Dirichlet prior distribution, the dependence structure that is imposed by only one parameter is a pattern of negative correlations between the multinomial probabilities of all categories. This negative correlation pattern is not appropriate in many practical situations. For example, if the categories represent the political parties in some country, with the corresponding multinomial probabilities reflecting the voting behavior in the next elections, a high probability of voting for a specific radical party might be accompanied by rather high probabilities of voting for other radical parties. A suitable prior in this case should imply a positive correlation pattern between parties of the same political views.
In the case of multinomial models, a particular difficulty is to elicit assessments that satisfy all the constraints of mathematical coherence. For example, the probabilities of each category must be non-negative and sum to one, which we refer to as the unit-sum constraint. Other constraints are implicit and less intuitive. We believe that a well-designed elicitation method should lead to coherent assessments without the expert having to be conscious of coherence constraints. This is not straightforward, especially with dependent variables in multivariate distributions. The expert has no built-in prior distribution to assess directly. Instead, she must be asked to quantify her opinions through questions she can both comprehend and answer. We believe it is important that the expert should be able to focus solely on her opinions when giving assessments.
To date, perhaps due to the difficulty discussed above, elicitation methods for multinomial sampling have mostly been proposed for modelling opinion by a Dirichlet distribution. An exception is the work by Elfadaly and Garthwaite (2013), who model opinion by a Connor-Mosimann distribution. However, this form of generalized Dirichlet distribution provides a structure that is only slightly richer than the standard Dirichlet distribution. Methods that yield a more flexible distribution than the Dirichlet distribution are long overdue. The present paper focuses on elicitation methods that yield a much more flexible distribution than the Dirichlet prior distribution.
Motivated by the deficiencies of the standard Dirichlet distribution, several authors have constructed new families of distributions for proportions that allow more general types of dependence structures. These distributions include, for example, the generalized Dirichlet distribution (Connor and Mosimann 1969), the multivariate normal distribution of log contrasts (Leonard 1975), and the multivariate logistic normal distribution (Aitchison 1986;Elfadaly and Garthwaite 2016).
A possible general multivariate distribution, that can serve as a prior distribution for multinomial models, is constructed through using a multivariate copula function. A copula is a multivariate function that represents a joint multivariate cumulative distribution function (CDF) in terms of one-dimensional marginal CDFs. Hence, it joins marginal distributions into a multivariate distribution that has those marginals. The theoretical foundation of copula functions was introduced in Sklar (1959), where he gave Sklar's Theo-rem which states that any joint distribution can be written in a copula form. The marginal distributions can thus be chosen independently from the dependence structure that is represented by the copula function. For an introduction to copulas, see for example Joe (1997), Frees and Valdez (1998) or Nelsen (1999).
The use of copula functions to elicit multivariate distributions has been considered in the literature; see Jouini and Clemen (1996) or Clemen and Reilly (1999), among others. In general, the joint distribution can be elicited by first assessing each marginal distribution. Then the dependence structure is elicited through the copula function. Different families and classes of copula functions have been defined for both bivariate and multivariate distributions. Jouini and Clemen (1996) used bivariate and multivariate Archimedean families and Frank's families of copulae to aggregate multiple experts' opinions about a random quantity. The book by Kurowicka and Cooke (2006) lists and compares different copula functions and discusses many other methods of modelling high-dimensional dependence. However, the simplest and most intuitive family of copulae is the inversion copula (Nelsen 1999). This has the form where G i are the known marginal CDFs of x 1 , . . . , x k , F (1,...,k) is their assumed multivariate CDF, and its marginals are F i . Hence, the marginal functions G 1 , . . . , G k are coupled through F (1,...,k) , via the quantile functions evaluated for each marginal, into a new multivariate distribution given by the copula function C. It has to be mentioned here that the arguments of the inversion copula function in (1) are the quantile functions F −1 i [.] (i = 1, . . . , k). It is therefore appealing to exploit the inversion copula in our elicitation method in which we ask the expert to assess only medians and quartiles. Experimental work in the literature suggests that the assessment of quantile functions is often the best way to elicit information from experts. See for example Garthwaite et al. (2005) and the references therein. Specifically, we believe that medians and quartiles are easier for an expert to assess than other quantile functions as they can be obtained by the first two steps of equally likely subdivisions using the bisection method (Garthwaite et al. 2005).
The distribution F (1,...,k) in (1) has sometimes been taken as a multivariate t-distribution (Demarta and McNeil 2005), or even as a Dirichlet distribution (Lewandowski 2008). Most commonly though, it is selected as a multivariate normal distribution with a joint CDF, F (1,...,k) = Φ k,R , and standard normal marginal CDFs, F i = Φ. This gives a Gaussian copula (Clemen and Reilly 1999) that is parameterized by the correlation matrix R of a multivariate normal distribution. To elicit R, Clemen and Reilly (1999) suggest that a pairwise rank-order correlation between each X i and X j , such as Spearman's ρ i, j or Kendall's τ i, j , should be assessed. Then properties of the multivariate normal distribution are used to transform them into the product-moment Pearson correlation r i, j as follows: The product-moment correlation matrix R is then formed from the elements r i, j . In practice, however, it may be preferable to elicit correlations on a numerical scale rather than rank-order correlations. In this paper, we use quartile assessments to obtain the product-moment correlation matrix R. Clemen and Reilly (1999) favour eliciting rank-order correlations, not product-moment Pearson correlation, on the grounds that the latter cannot necessarily be transformed through the function Φ −1 [G i (.)]. In contrast, rank-order correlations are transformation respecting regardless of the choice of the marginal CDF G i (.), i.e. they are invariant under strictly monotone increasing transformations of the form Φ −1 [G i (.)], for any G i (.). To elicit these correlations, Clemen and Reilly (1999) mention three methods that can be used either separately or together. However, none of the methods is guaranteed to yield a positive-definite correlation matrix. On the other hand, the elicitation method of Kadane et al. (1980) has been designed to encode a positive-definite covariance matrix of a multivariate t-distribution as a conjugate prior for the hyperparameters of a normal multiple linear regression model. Their method can be useful in a variety of multivariate elicitation problems that require the assessment of positive-definite matrices (Garthwaite et al. 2005. These clearly include the problem addressed here, so that we modify their method to form a Gaussian copula. Interestingly, no other elicitation method for Gaussian copulas seems to encode a positive-definite correlation matrix using quartile assessments. To fill this gap, the second proposed method in this paper elicits a Gaussian copula function with beta marginal distributions as a prior distribution for multinomial models. Our approach simultaneously overcomes two problems of the method of Clemen and Reilly (1999). First, we transform the assessed conditional quartiles of X i and X j , through Φ −1 [G i (.)], which enables productmoment correlations to be computed on the normal scale without the need of rank-order correlations. Second, the conditional quartiles are assessed according to the structural elicitation procedure of Kadane et al. (1980), which guarantees that the elicited correlation matrix is positivedefinite.
The paper is organized as follows. In Sect. 2, we discuss our choice of distribution, the Dirichlet prior and the Gaussian copula prior, for representing expert uncertainty relevant to proportions of a categorical variable. We show how this can be used to specify a prior for the probabilities of a multinomial distribution. In Sect. 3 we describe the assessment tasks that an expert performs, and in Sect. 4 we derive hyperparameters of both forms of prior distributions from the elicited assessments. Sections 3 and 4 also describe an interactive graphical interface that has been developed to implement the two proposed elicitation methods. In Sect. 5, these methods and the implementing software are used by an environmental ecologist to quantify his opinion about the future proportions of different phenotypes of a specific snail species in the UK. Concluding comments are given in Sect. 6.
2 Two multivariate prior distributions for a multinomial model

The multinomial likelihood
Let the random vector W = (W 1 , . . . , W k ) be multinomially distributed with k categories, n trials and a vector of probabilities p = ( p 1 , . . . , p k ), so that where w i = n, 0 ≤ w i ≤ n, 0 ≤ p i ≤ 1 and p i = 1.

The Dirichlet prior distribution
A conjugate prior for the parameter vector p is the Dirichlet distribution of the form where Γ (.) is the gamma function defined by Γ (a) = ∞ t=0 t a−1 e −t dt, N = a i and a i are the hyperparameters representing the relative size of each category, so a i > 0.
The expectation and variance of p i are a i /N and a i (N − a i )/[N 2 (N + 1)] and the covariance between p i and p j is −a i a j /[N 2 (N + 1)], where i = 1, . . . , k; j = 1, . . . , k and i = j. The marginal distribution of each p i (i = 1, . . . , k) is a beta distribution (e.g. Wilks 1962), where To elicit the vector of hyperparameters a = (a 1 , . . . , a k ), we will exploit this direct relationship between the Dirichlet distribution and its marginal beta distributions.

Constructing a Gaussian copula prior distribution
Rather than assuming that the marginal beta distributions in (5) stem from a Dirichlet distribution, we would like to allow a more flexible dependence structure via another joint distribution, with the aim of better representing the expert's opinion. A flexible tool for this task is given by the copula function, as discussed in Sect. 1. This allows us to choose the marginal distributions independently from the dependence structure between them. The best-known example of the inversion copula in (1) is the Gaussian copula (Clemen and Reilly 1999). For any random vector X = (X 1 , . . . , X k ), the Gaussian copula is defined at the point (x 1 , . . . , x k ) as Here Φ k,R is the CDF of a k-variate normal distribution with zero means, unit variances, and a correlation matrix R that reflects the desired dependence structure. The function Φ is the marginal standard univariate normal CDF, and G i (i = 1, . . . , k) are any chosen marginal CDFs. The marginal mean and variance of each x i are expressed through its marginal CDF, G i (x i ), hence, no other parameters are required for Φ k,R beyond the correlation matrix R. To be used as a prior distribution, the density function of the Gaussian copula is needed. Fortunately, since Φ k,R and Φ are differentiable, the Gaussian copula density function is easily obtained by differentiating (7) with respect to x i (i = 1, . . . , k); see for example Clemen and Reilly (1999). This gives where is the density function corresponding to the choice of G i (.) (i = 1, . . . , k), and I k is the identity matrix of order k.
To construct a Gaussian copula function in the case of a multinomial model, it might seem natural to assume that each marginal distribution of p i (i = 1, . . . , k) is a beta distribution, as in (5), and try to construct a Gaussian copula function for the multivariate distribution of p. However, applying the normalizing transformations is the assessed beta CDF, would raise three problems. First, encoding marginal beta distributions directly on the elements of p might result in mathematically incoherent marginal distributions with mean values that do not sum to one. Second, the inverse normalizing transformations , while preserving the marginal domain of each probability, 0 ≤ p i ≤ 1, do not necessarily fulfil the unit-sum constraint k i=1 p i = 1. Third, under the unit-sum constraint on p, the normally distributed vector y k = (Y 1 , . . . , Y k ) contains a redundant variable and hence leads to a singular multivariate normal distribution.
To overcome these problems, we define new variables Z 1 , . . . , Z k , where The inverse transformations are given by Instead of assuming a beta marginal distribution of the form (5) for each p i (i = 1, . . . , k), we assume that each Z i (i = 1, . . . , k − 1) has a marginal beta distribution with parameters α # i and β # i to be encoded as in Sect. 4.1. Then we construct a Gaussian copula function for the multivariate distribution of z k−1 = (Z 1 , . . . , Z k−1 ). The beta marginal distributions preserve the domain of each Z i (i = 1, . . . , k − 1) on [0, 1], while the Gaussian copula accounts only for their dependence structure. Here, no other constraints need to be fulfilled on z k−1 , while the unit-sum constraint on p is always attained by virtue of the transformations in (9) and (10).
Using the Gaussian copula function, the dependence structure of the multivariate distribution of Z i , and consequently of p i , will have high flexibility rather than the limited dependence structure imposed by the Dirichlet distribution. In the underlying multinomial model, each multinomial parameter p i represents the probability that an observation falls in the ith category. Hence, Z i (i = 1, . . . , k − 1) can be defined as the probability that an observation falls in category i, given that the first i − 1 categories have been eliminated and the available categories are only the last k − i + 1 categories. Connor and Mosimann (1969) . Thus the generalized Dirichlet distribution can be regarded as a special case of the Gaussian copula in which Z 1 , . . . , Z k−1 are independent. If, moreover, , so the standard Dirichlet distribution is also a special case of the Gaussian copula prior. This confirms that the Gaussian copula function imposes a prior distribution on p with a dependence structure that is more flexible than that of the standard and generalized Dirichlet distributions.
The prior distribution for p is defined in terms of the prior distribution of z k−1 . Let G i (Z i ) be the CDF of the beta distribution of Z i with hyperparameters α # i and β # i (i = 1, . . . , k − 1) to be elicited in Sect. 4.1. We assume that the joint density of z k−1 is given by a Gaussian copula density, . The correlation matrix R and the parameters of the beta densities are the hyperparameters of the Gaussian copula prior distribution.
Note that the marginal distributions of this joint density are still the same beta marginal distributions. After eliciting the hyperparameters of the beta distribution for each Z i (i = 1, . . . , k − 1) as detailed in Sect. 4.1, the prior distribution is totally known except for the matrix R. Here, R is the correlation matrix of the multivariate normally distributed vector y k−1 ; it needs to be encoded effectively and must be a positive-definite matrix. The aim is to choose R so as to model the expert's opinion about the dependence between the p i s. In Sect. 4.3, we introduce a method, inspired by Kadane et al. (1980), to encode the correlation matrix R. Although the density in (11) is neither a multivariate normal for p k−1 = ( p 1 , . . . , p k−1 ) nor z k−1 = (Z 1 , . . . , Z k−1 ), and the correlation matrix of z k−1 is not an explicit function in R, the dependence structure between the elements of p k−1 , and consequently that of Z k−1 , is indirectly reflected by R. Moreover, we can still use the multivariate normal properties to form a positive-definite matrix R by considering the following normalizing transformations, Under the main assumption of the Gaussian copula construction, and from (12), the vector y k−1 = (Y 1 , . . . , Y k−1 ) has a multivariate normal distribution with zero means, unit variances and a correlation matrix R, so y k−1 ∼ MVN(0, R). In our proposed approach, the matrix R is encoded as a covariance and a correlation matrix of the multivariate normal random vector y k−1 . This is achieved by utilizing the monotone increasing property of the transformations in (12). We assess conditional quartiles of p k−1 , then transform them into those of z k−1 and y k−1 using (9) and (12), respectively.
Correlation coefficients between the elements of y k−1 can then be estimated using their conditional quartiles and utilizing properties of the multivariate normal distribution. This is described in Sects. 3.2.2, 3.2.3 and 4.3.

Assessment tasks
With both the Dirichlet and Gaussian copula prior distributions, marginal distributions are each of the beta type. In both cases, the first set of assessments determines the hyperparameters of the marginal beta distributions. For the standard Dirichlet distribution, as detailed in Sect. 3.1 the expert is asked to assess marginal medians and quartiles of the probability of each category. These quantiles are the median and quartiles of the expert's subjective probability distribution for each p i ; the distance between the quartiles reflects her uncertainty about the true value of p i . The marginal median and quartile assessments are obtained in Sect. 3.1.1. They are used to encode the marginal beta distribution of each p i (i = 1, . . . , k). In the second part of the elicitation process, detailed in Sect. 3.2, we obtain the required assessments to construct a Gaussian copula prior distribution. Specifically, in Sect. 3.2.1, we obtain the marginal median and quartile assessments that are required to encode the marginal beta distributions of Z i (i = 1, . . . , k − 1). Further sets of conditional medians and quartiles are obtained in Sects. 3.2.2 and 3.2.3; they will be used to encode the correlation matrix R of the Gaussian copula prior.
To ensure that the expert understands the request to assess her median estimate of a proportion, she is told: "You should consider it equally likely that the true proportion is above the median or below the median. For example, suppose you assess the median as 0.6, you should think it equally likely that the proportion will be above 0.6 as that it will be below 0.6. Hence, if you were asked 'Would I rather bet on the proportion being above 0.6 or below 0.6?', you should find it hard to decide." For the lower (upper) quartile assessment, the expert is told: "The lower (upper) quartile divides the probability below (above) the median in half. Suppose your assessed median is 0.3 and your assessed lower quartile is 0.2. Then your probability that the proportion is below 0.2 should equal your probability that it is between 0.2 and 0.3."

Obtaining marginal assessments for p
For each category in turn, the expert is asked to assess a lower quartile, a median and an upper quartile for the probability p i , i = 1, . . . , k, say L i , m i and U i , respectively. In the practical evolutional example detailed in Sect. 5, the expert assessed his median, lower and upper quartiles for the future proportion of each of the 6 different phenotypes of the most common snail species in the UK, see Table 1. Here we adopt a marginal approach, so that the categories are interchangeable and their order does not matter. In Sect. 4.1, the marginal medians and quartiles assessed here are used to encode the parameters α i and β i of each marginal beta distribution of p i for the standard Dirichlet distribution.

Obtaining marginal assessments for z k−1
For category i (i = 2, . . . , k − 1), the expert is asked to assume that the first i − 1 categories have been eliminated from the analysis and that she should only consider the remaining k − i + 1 categories, namely categories i, i + 1, . . . , k. That is, the expert should assume that the observation does not fall in one of the first i − 1 categories, but which of the remaining category it falls in is unknown. Given this situation, the expert is asked to assess the median of the probability of each of the categories i, i + 1, . . . , k.
Only the assessment for the first of these categories is used and we denote it m # i,0 . The other assessments are elicited to improve internal consistency; assessments of means should sum to one so the assessments of the medians should be close to one.
In fact, these assessments represent the marginal unconditional medians of each Z i = p i /(1 − i−1 j=1 p j ), for i = 2, . . . , k −1. After these medians have been assessed, the expert assesses the marginal unconditional lower and upper quartile of each Z i . We denote these assessments as L # i,0 and U # i,0 (i = 2, . . . , k − 1). Clearly, since Z 1 = p 1 , we have L # 1,0 = L 1 , m # 1,0 = m 1 and U # 1,0 = U 1 and these are assessed If a Gaussian copula prior distribution is to be encoded, we define m * 1,0 as the median of p 1 , and m * i,0 as the median of ( . Instead, we compute their values from the corresponding marginal medians, m # i,0 , of Z i (i = 1, . . . , k − 1) assessed above. This can be done by utilizing Lemma 1 (below) to obtain m * i,0 in terms of m # i,0 . The former will be used as the conditioning values in Sects. 3.2.2 and 3.2.3.
Lemma 1 Under the unit-sum constraint on p, and the multivariate normality of y k−1 , we have See the Appendix for a proof of this lemma. Assessments obtained in this section are used in Sect. 4.1 to encode the hyperparameters α # i and β # i of the marginal beta distribution of Z i (i = 1, . . . , k − 1) for the Gaussian copula prior distribution.

Assessing conditional quartiles
After finishing the marginal assessment tasks in Sect. 3.2.1, the expert is asked to give conditional quartile assessments for specified categories under the assumption that the initial assessments are the correct values for the other categories. Specifically, the expert is asked to assume that p 1 equals its initially assessed median value, i.e. p 1 = m * 1,0 , and then gives a lower quartile L * 2 and an upper quartile U * 2 for p 2 . For each remaining p j , j = 3, . . . , k − 1, the expert assesses the two quartiles L * j and U * j conditional on p 1 = m * 1,0 , . . ., p j−1 = m * j−1,0 . For the last category, the lower (upper) quartile L * k (U * k ) of p k is automatically shown to the expert once she assesses the upper (lower) quartile U * k−1 (L * k−1 ) of p k−1 . The two quartiles L * k and U * k are shown to the expert as a guide to help her choose L * k−1 and U * k−1 . When p 1 = m * 1,0 , . . . , p k−2 = m * k−2,0 , from the unit-sum constraint it follows that L * k−1 + U * k = U * k−1 + L * k = 1 − m * 1,0 − · · · − m * k−2,0 . In practice, these conditional assessments constitute an easy task for the experts to perform. They normally think of successive category assessments in a conditional way given their already assessed values for the preceding categories. For conditional assessments, the order of the categories is important. We recommend ordering them from the most probable category to the least probable, unless there is a natural order to the categories. This ordering leads to assessment tasks that involve less skewed distributions. See Elfadaly and Garthwaite (2013) for more detail and examples on this recommended category ordering.
L # i and U # i are used to estimate the correlation matrix R of the Gaussian copula prior, as will be discussed in Sect. 4.3.

Assessing conditional medians
First the expert is asked to assume that the true value of p 1 is m * 1,1 , while her median assessment of p 1 had been m * 1,0 . The value specified for m * 1,1 does not equal m * 1,0 and we let η * 1 = m * 1,1 − m * 1,0 . Given the information that p 1 = m * 1,1 , the expert is asked to change her previous medians m * j,0 of each p j , and her new (conditional) median is denoted by m * j,1 ( j = 2, . . . , k). In practice, this question can be phrased to the expert as: "Suppose that the proportion in the first category were actually changed from your original assessment to the value given on the interactive software graph. If that were the case, what would be your median estimates of the proportions in the other categories?" Here, the expert is asked to sequentially assume that each conditional median assessment she gives at this stage is the true value and that she should condition on each of them to be true while she is assessing the remaining conditional medians. For example, when the expert is asked to assess m * 4,1 , she does not only assume that p 1 = m * 1,1 , but she also assumes that p 2 = m * 2,1 and p 3 = m * 3,1 . In each successive step i, for i = 2, . . . , k −2, the expert is asked to suppose that the true value of p i is m (while her initially calculated median of p i had been m * i,0 ) and that the true values of p 1 , . . ., p i−1 are equal to their initially calculated medians, m * 1,0 , . . . , m * i−1,0 , respectively. Given this information, the expert is asked to revise her initial median values m * j,0 of p j , and the new assessments are denoted by m * j,i ( j = i + 1, . . . , k). Again, as explained in the case of i = 1, these assessments are also sequentially conditional assessments for i = 2, . . . , k − 2, i.e. m * j,i is assessed under the condition that p 1 = m * 1,0 , . . . , The assumed conditioning true value m * i,i should be (a) reasonably different from the initial value m * i,0 , so that the expert revises the previously calculated assessments at the other categories by a measurable change, but (b) not too different from m * i,0 , so that it is still a plausible value according to the expert's belief. Moreover, as will be shown later, for mathematical coherence the arbitrary values η * i s have to be chosen such that η * i = 0 (i = 1, . . . , k − 2). To fulfil these constraints, the lower quartile of each category is taken as the value for m * i,i . Hence, we put m * 1,1 = L * 1,0 and m * . The expert uses the software to make her assessments on an interactive graph. An example is illustrated in Fig. 1, where the categories represent 6 phenotypes of snails: yellowbanded/ yellow-unbanded/ pink-banded/ pink-unbanded/ brown-banded/ brown-unbanded. The conditioning set of median values are the three right hand side (red) bars at the first three categories. The expert is asked to assess how her new median values, shown as the three right hand side (one blue and two orange) bars at the last three categories, will change based on the new conditioning set. The initially calculated medians of all categories are shown as the left hand side (grey) bars. For mathematical coherence, the expert's assessments must fulfill the first statement of Lemma 2 below. The second statement of Lemma 2 gives the conditional medians of Z j , say m # j,i , in terms of the corresponding conditional medians of p j , m * j,i ( j = i + 1, . . . , k). The former are used in Sect. 4.3 to estimate the correlation matrix R of the Gaussian copula prior distribution.
Lemma 2 Under the unit-sum constraint on p, and the multivariate normality of y K −1 , for i = 1, . . . , k − 2, we have See the Appendix for a proof of this lemma.
By scaling their values, the software normalizes the values of m * i+1,i , . . ., m * k,i so that they satisfy Lemma 2. Specifically, the normalized value of m The expert is told that the conditional medians must satisfy a mathematical equation and that a good representation of her assessments has been calculated. The software displays the normalized conditional medians as the thinner (yellow) bars in Fig. 1. The expert has the option of changing these medians until she feels that the suggested normalized set m * i+1,i , . . ., m * k,i gives an acceptable representation of her opinion. In Fig. 1, the expert's conditional assessments are the three right-most (blue and orange) bars for the last three categories. These are reasonably coherent, so the suggested values in the thinner (yellow) bars are almost equal to these assessments. The current assessment task stops at step k − 2, as we do not ask for any conditional assessments for the last remaining category p k . Since the condition of summing to one should always be fulfilled, conditioning on specific values of all p 1 , . . . , p k−1 gives a fixed value for p k .

Encoding the hyperparameters 4.1 Encoding the hyperparameters of the beta marginal distributions
To encode a standard Dirichlet prior distribution, the expert has already assessed a lower quartile, L i , a median, m i , and an upper quartile, U i , for p i (i = 1, . . . , k), as discussed in Sect. 3.1.1. But, to encode a Gaussian copula prior distribution, the expert has already assessed a lower quartile, L # i,0 , a median, m # i,0 , and an upper quartile, U # i,0 , for Z i (i = 1, . . . , k − 1), as discussed in Sect. 3.2.1. The quartile assessments L i , m i and U i are used here to encode the hyperparameters α i and β i of the marginal beta distribution of each p i (i = 1, . . . , k). Similarly, L # i,0 , m # i,0 and U # i,0 are used to encode the hyperparameters α # i and β # i of the marginal beta distribution of each Z i (i = 1, . . . , k − 1). The same method of encoding beta hyperparameters is used in both cases, so we denote a general set of quartile assessments byL i,0 ,m i,0 andŨ i,0 and show how we use them to encode the beta hyperparameters denoted byα i andβ i .
The method proposed in Elfadaly and Garthwaite (2013) is used to encode the two parametersα i andβ i of each marginal beta distribution. In their approach, the three assessed quartilesL i,0 ,m i,0 andŨ i,0 are substituted in a normal approximation of the beta distribution to obtain initial values of the beta parameters. These are then used as the initial values in a numerical least-squares optimization which encodes the beta parameters as non-negative valuesα i andβ i that minimize where G(q;α i ,β i ) is the CDF of a beta distribution with parametersα i andβ i at the point q. The least-squares optimization enables fast and accurate computation. See, for example, Albert et al. (2012) and the references therein. If a Dirichlet prior is to be encoded, the parameters α i = α i and β i =β i will be directly reconciled to obtain the Dirichlet hyperparameters in Sect. 4.2. But for encoding a Gaussian copula prior, we take α # i =α i and β # i =β i ; these encoded beta marginal distributions of Z 1 , . . . , Z k−1 have no constraint to satisfy-under the transformations in (9) and (10), the unit-sum constraint on p is always fulfilled for any encoded beta marginal distributions of Z 1 , . . . , Z k−1 .

Encoding the Dirichlet hyperparameters
After beta marginal distributions have been assessed for each category's probability p i (c.f. Sect. 4.1), the beta parameters are then adjusted to estimate the Dirichlet hyperparameter vector as follows. First, note that the system of equations in (6) does not guarantee a consistent solution for a. This is because the expert will seldom be fully consistent and give marginal assessments that are exactly mathematically coherent. From (6), each marginal step of the elicitation process provide two estimates for a i and N i , namely, for i = 1, . . . , k, we get a i = α i and N i = α i + β i = k j=1 a j . Moreover, the estimated hyperparameters must fulfill the unit-sum constraint on the expected values of the probabilities. Specifically, if the expected values are denoted by μ i (where μ i = a i /N i , for i = 1, . . . , k) they must satisfy k i=1 μ i = 1. If the assessments are not coherent and do not satisfy these constraints, which is usually the case, initial assessments will need to be reconciled. Even if the expert were very consistent and gave reasonably coherent assessments, reconciliation would be needed for correcting small rounding errors. Lindley et al. (1979) investigated the reconciliation of assessments that are inconsistent with the laws of probabilities (incoherent). They developed least-squares procedures as reconciliation tools that may be used for any expert's incoherent assessments. In the spirit of their approach, we propose the following three options to reconcile incoherent estimates of the μ i and N into mathematically coherent estimates μ * i , N * , respectively.
(I) Normalize μ i 's into μ * i 's that add up to one, i.e. put μ * i = μ i / k j=1 μ j , for i = 1, . . . , k. (II) Minimize the sum of squares of differences between μ * i and μ i , subject to the constraints 0 < μ * i < 1 and k i=1 μ * i = 1. So, the μ * i are taken as the values that minimize where λ is a Lagrange multiplier. (III) The precision of each p i , i.e. the inverse of its variance, can be used as a weight to reflect the expert's confidence about each of her assessments (Lindley et al. 1979). These weights are used in a weighted least-squares procedure. Hence we determine the μ * i s that minimize With (I) and (II), N * is set equal to the average of the N i 's, while with (III) it is set equal to their weighted average as Estimating μ * i and N * , using any of the options listed above, makes it easy to estimate a i by a * i , where a * i = μ * i N * for i = 1, . . . , k. For options (II) and (III), the software computationally estimates λ as zero up to a tolerance of 10 −4 .
The software elicits three hyperparameter vectors of the Dirichlet distribution, one vector for each of the above options. Each vector is then used to compute the corresponding pairs of marginal beta parameters as given in (6). Three quartiles of each beta marginal distribution are computed numerically for each Dirichlet hyperparameter vector. The three sets of quartiles are then displayed to the expert and she is asked to select the set of quartiles that best represents her opinion. See Fig. 2: the middle three button groups at the bottom of the screen offer the options that the expert can choose for her quantile assessments. (The other buttons are for navigation to the previous screen or the next screen.) The vector with the selected set of quartiles will be taken Fig. 2 A feedback screen for Dirichlet elicitation. The following message is shown to the expert: The probabilities of the different categories must add to one. The quartile assessments for the different categories must also meet certain requirements in order to be internally consistent. This program gives three options for reconciling your assessments to meet these requirements. Select the option that best matches your opinion and then click 'Next' as the final elicited hyperparameter vector of the Dirichlet prior.
The expert is still able, however, to modify any or all of the selected set of quartiles, in which case beta parameters are computed again and the final Dirichlet hyperparameter vector is computed through averaging, as in option (I) above.

Encoding the Gaussian copula hyperparameters
The normalizing one-to-one functions in (12) are used to transform the assessed conditional quartiles of z k−1 into conditional quartiles of y k−1 , yielding conditional expectations, variances and covariances of the multivariate normal variables. To form a positive-definite correlation matrix R, let y i = (Y 1 , . . . , Y i ) and R i = Var(y i ), for i = 1, . . . , k − 1, where R 1 = Var(Y 1 ) = 1 and the final matrix R = R k−1 . Supposing that R i−1 has been estimated as a positive-definite matrix, an encoding method, based on that of Kadane et al. (1980), is used to form a positive-definite correlation matrix R i (i = 2, . . . , k − 1). The mathematical details of the method are shown in the Appendix. Since R 1 = 1 > 0, by mathematical induction the full correlation matrix R = R k−1 is guaranteed to be positive-definite.
We have to note that, according to this method of elicitation, the variances on the main diagonal of R, say r i,i , i = 1, . . . , k − 1, cannot be guaranteed to each equal one, except for the first element r 1,1 . It is easy, however, to transform R into R * = ARA, where A is a diagonal matrix with a 1,1 = 1 and a i,i = 1/ √ r i,i , for i = 2, . . . , k − 1. Then R * can be used as the correlation matrix in the Gaussian copula function and satisfies both the unit variances and positive-definiteness requirements. The unit variances in the correlation matrix R * guarantee that each marginal distribution G i (Z i ) is still a beta distribution with the same hyperparameters α # i and β # i that were encoded, for i = 1, . . . , k − 1, in Sect. 4.1.

Evaluating the encoded priors
A good elicitation method encodes a prior distribution that captures the expert's knowledge and beliefs faithfully and as closely as possible. Thus, a reasonable way to compare the performance of the two encoded prior distributions is to check which of them produces values that are closer to the expert's original assessments. The aim of the Gaussian copula prior is to capture the correlation structure better than the Dirichlet prior, so the conditional medians computed from each encoded prior distribution can be compared to the expert's original assessments of the conditional medians to check whether this aim has been achieved. We implement this approach at the end of the evolution example in Sect. 5 to evaluate the encoded priors.
A different quality of probability assessments is calibration, which measures the agreement between an expert's assessed probabilities and reality. Events that are given a prior probability of occurring of, say, 0.3 should ideally occur about 30 % of the time. Similarly, a set of subjective 50 % prediction intervals should contain the target event half the time. A body of psychological research has examined the calibration of probability assessments in varied contexts; see (O'Hagan et al. 2006, chapter 4) for a review. How-ever, as noted by (O'Hagan et al. 2006, p. 161), "comparing the elicited probabilities or distributions with the true values of the events or variables in question does not evaluate the success of elicitation; it does not compare the elicited probabilities with the expert's underlying opinions."

Example: Evolution of the snails in the UK
Environmental scientists are interested in the effect of climate change on the evolutionary responses of different species. The brown-lipped banded snail (Cepaea nemoralis) is known to have a short generation time and tends to sensitively adapt its shell colour and banding pattern according to its thermal environment. Silvertown et al. (2011) tested the evolutionary changes in this type of snail by comparing a recent dataset with historical datasets from the last century. Data on the frequency of each phenotype of this species were collected by volunteers in 15 European countries through the citizen science project Evolution MegaLab (Worthington et al. 2012). Thousands of public volunteers contributed data to the project by recording the shell colour and banding pattern of the snails they had observed in each country.
An environmental expert who has worked with these data is interested in the modelling of ecosystems under future climate change. In this example, the expert (who is male) used the elicitation software to quantify his opinion about the expected proportions of each phenotype of Cepaea nemoralis snails in the UK in 2050 under the assumption that only a little is done to reduce climate change (a scenario called Representative Climate Forcing 6, under which greenhouse gas emissions continue to rise without strong mitigation). He advised that three main phenotypes can be considered: yellow, pink or brown; each of them can be banded or unbanded. Hence, any single snail can be categorized in one of six different phenotypes, so the problem can be formulated as a multinomial model with six categories. The probabilities p 1 , . . . , p 6 are the proportions of each phenotype that might be observed in the UK. Our method and software were used to quantify the expert's opinion about his uncertainty regarding the parameters of this multinomial model, representing his opinion by both a Dirichlet prior distribution and a Gaussian copula prior.
After initializing the software and defining the model, the expert assessed his medians, m i (i = 1, . . . , 6), of the proportion of each of the following 6 phenotypes of Cepaea nemoralis snails: yellow-banded/ yellow-unbanded/ pinkbanded/ pink-unbanded/ brown-banded/ brown-unbanded. Then the expert assessed lower and upper quartiles (L i and U i ) for the proportion of each phenotype. His assessed medians and quartiles are given in Table 1.
These assessments were used to encode a marginal beta prior distribution, with hyperparameters α i and β i , for the proportion of phenotypes in each category using equation (13) in Sect. 4.1. These can then be used to determine the hyperparameters of a Dirichlet distribution using any of the three least-squares options detailed in Sect. 4.2. The expert chose option (II) where the values of the Dirichlet hyperparameters were encoded as given in Table 2.
The median values and quartiles of the coherent marginal beta distributions of the Dirichlet variates were computed and presented to the expert as feedback in Fig. 3. These are shown as the left hand side (grey) bars for each category. During this feedback stage he was invited to accept or revise these quantities by clicking on the right hand side (blue) bars of each category in Fig. 3. The initial median values, m i (i = 1, . . . , 6), as assessed by the expert have a sum that is nearly equal to one, so the coherent medians suggested by the software in Fig. 3 were close to his median assessments as shown in Table 1. However, the expert felt that some of the suggested quartiles were too far apart to be a reasonable representation of his opinion, and he chose to revise the suggested quartile values for the last three categories. These were then used to encode another Dirichlet distribution with hyperparameters given in Table 3 and this is taken as the expert's Dirichlet prior distribution. Table 4 shows the final coherent set of medians and quartiles after the expert's revisions computed from the Dirichlet hyperparameters in Table 3.
The expert's initial assessments in Table 1 give hyperparameter values α i and β i of the marginal beta distributions that result in different estimates of the Dirichlet hyperparameter N . Specifically, they give N i = α i + β i (i = 1, . . . , 6) that range between 11.75 and 73.27. The variability of these values can be used as a diagnostic check to determine whether the Dirichlet prior distribution might be a reasonable representation of the expert's knowledge. As concluded Fig. 3 The coherent marginal quartiles for the Dirichlet distribution given as the left (grey) boxes, with the expert's revised values, given as the right (blue) boxes. The following message is shown to the expert: This screen shows the medians and quartiles that the software has determined for your prior distribution. You should change them if they do not form a reasonable representation of your opinion. The grey (left) boxes are there to show the coherent quartiles in Elfadaly and Garthwaite (2013), a substantial variation in the values of N i (i = 1, . . . , k) indicates that a Dirichlet prior does not capture the expert's opinions well. This suggests that here the Dirichlet prior may not be an adequate representation of the expert's opinion, in which case the Gaussian copula prior is to be preferred. The expert gave further assessments to encode a Gaussian copula prior as well, which gains much greater flexibility in the correlation structure of the parameters of his multinomial model. To encode a Gaussian copula prior, the expert was asked to assess a median value for the proportion of the yellow-banded phenotype of snails in the UK in 2050. He assessed this as m # 1,0 = 0.446, the median of Z 1 = p 1 .
As mentioned before in Sect. 3.2.1, m # 1,0 m 1 = 0.449. To improve the internal consistency between this assessment and the forthcoming ones, the expert was also asked at this stage to assess the median values of the proportions of all other phenotypes. These were assessed as 0.198, 0.147, 0.076, 0.091, 0.041. Clearly, these assessments are quite consistent with the experts's assessments m i in Table 1.
The expert was then asked to assume that a randomly selected snail is not yellow-banded. For this snail he was asked to assess the median probability, m # 2,0 , that it belongs to the yellow-unbanded category. To help the expert, the software uses the median values assessed at the previous stage to suggest a value for the median of Z 2 = p 2 /(1 − p 1 ) as m # 2,0 = 0.198/(1 − 0.446) = 0.358. Suggested median values for the remaining categories were also computed and presented to the expert (See Fig. 4). He accepted the value of m # 2,0 = 0.358 and all the other suggested median values as a reasonable representation of his opinion.
Then at each step i (i = 3, 4, 5), assuming that a selected snail does not belong to the first i − 1 phenotypes, the expert accepted the suggested value of m # i,0 as the median of , together with the remaining median values at each step.
Similarly, the expert was asked to assess a lower and an upper quartile (L # i,0 and U # i,0 ) for the probability Z 1 = p 1 and for Z i (i = 2, . . . , 5) as the probability of belonging to the ith phenotype given that the snail could not be categorized in the first i − 1 categories. The full set of median and quartile assessments obtained at this stage are shown in Table 5. These assessments were used in the optimization equation (13) . . . , 5). The values of these hyperparameters are given in Table 6.
To encode a correlation matrix for the Gaussian copula prior, the expert gave conditional assessments that quantified his opinion about the dependence structure between the marginal beta distributions. Specifically, he assessed conditional quartile values, L * i and U * i , for each ith (i = 2, . . . , 5) category, under the condition that the assessed medians for the previous j ( j = 1, . . . , i − 1) categories were actually their true values. The expert's five pairs of assessments for the conditional lower and upper quartiles are given in Table 7. The quartiles for the last category in Table 7 were automatically computed by the software when the expert assessed two quartiles for the fifth category. This is also illustrated in Fig. 5. Next, conditional on the proportion for the first category being m * 1,1 = 0.349 (his assessment of its lower quartile, L # 1,0 ), the expert gave conditional median assessments, m * j,1 ( j = 2, . . . , 6), for the proportions of the five remaining categories. The number of conditions was then increased in stages. Table 8 gives all the conditional median assessments, m * j,i (i = 1, . . . , 4; j = i + 1, . . . , 6), where the underlined values constitute the conditioning set at each stage. The assessments reported at each row of Table 8 were sequentially obtained under the assumption that the previously assessed Assessing conditional quartiles for the last two categories (L * 5 = 0.081 and U * 5 = 0.100). The following message is shown to the expert: Give your lower and upper quartile assessments on category five, assuming that your median assessments for the first four categories with (red) boxes are correct. The implied quartiles for the last category-the right-most (yellow) box-are also shown. To be logical, your assessment must be below the upper-most short horizontal (red) dotted line and should be between the other two short (orange) horizontal dotted lines  Fig. 1 shows the expert's assessments of conditional medians m * 1,4 , m * 1,5 and m * 1,6 for the last three categories, given that p 1 = m * 1,1 = 0.349, where the conditional median assessments m * 1,2 = 0.279 and m * 1,3 = 0.162 are also assumed to be the true values of p 2 and p 3 , respectively. This is one of the sub-stages required to obtain the assessments in the first row of Table 8. However, the expert chose not to change his conditional median assessments at different sub-stages. He only revised his assessments when the underlined conditioning values were changed, i.e. from one row to the next in Table 8.
This was the last assessment task, after which the software output the encoded hyperparameters, α # i and β # i , of the marginal beta prior distributions of Z i (i = 1, . . . , 5) as in Table 6. These were computed as detailed in Sect. 4.1 using Eq. (13). The dependence structure between these beta Table 9 The encoded correlation matrix, R * , of the Gaussian copula prior marginals was quantified as a multivariate Gaussian copula function with the encoded correlation matrix, R * , given in Table 9. This has been encoded as detailed in Sect. 4.3. The encoded matrix in Table 9 does not give covariances between the beta distributed variables, Z 1 , . . . , Z 5 . Instead, it gives the covariances between the transformed normal variates, Y 1 , . . . , Y 5 . The Gaussian copula multivariate distribution is parameterized by both the marginal beta parameters in Table 6 and the correlation matrix in Table 9. These determine the experts's Gaussian copula prior distribution. The software produces a file containing the WinBUGS model (Spiegelhalter et al. 2007) that uses the Gaussian copula prior distribution. This is given at the end of the Appendix.
To examine the assessed prior distributions, we compute the conditional medians from each of the two encoded prior distributions and compare them to the expert's assessments in Table 8. The aim is to check whether the Gaussian copula has captured the correlation structure of the multinomial probabilities better than the Dirichlet distribution. Table 10 gives the conditional medians computed from the Dirichlet distribution, with the encoded hyperparameters in Table 3, together The underlined values represent the conditioning set in each row. The discrepancies between the conditional medians of each prior and the expert's assessments are given as bracketed percentages with the corresponding conditional medians computed from the Gaussian copula prior, with the encoded hyperparameters in Tables 6 and 9. The original conditional medians assessed by the expert are given again in Table 10 to facilitate the comparison. For both prior distributions, the discrepancies of the computed conditional medians from the expert's assessments are also presented in the same table as percentages of the assessments. Comparing the conditional median values with the corresponding expert's original assessments it can be seen that the conditional medians computed from the proposed Gaussian copula prior are much closer to the expert's assessments than are those computed from the Dirichlet distribution. This illustrates the potential advantage from using the Gaussian copula prior rather than the Dirichlet prior distribution.
The elicitation process took about an hour to complete. The expert stressed the importance of the convenient order of categories when giving conditional assessments, commenting that ordering the categories in a suitable sequence made it easier for him to think about the conditions.

Concluding comments
Our aim in this paper was to propose methods for eliciting expert opinion about the hyperparameters of both the Dirichlet distribution and the Gaussian copula distribution. So as to simplify the separate tasks that an assessor must perform, the elicitation process for a multivariate prior has been decomposed into a number of processes for eliciting univariate beta distributions. The hyperparameters of these beta distributions are estimated from marginal median and quartile assessments obtained from a field expert. Two novel methods have been proposed for reconciling the beta mar-ginal distributions into a multivariate prior distribution. The first method reconciles the beta distributions as marginals of the well-known standard conjugate Dirichlet distribution, giving a prior distribution that is easily applied in practice.
In the second method, instead of assuming a Dirichlet prior, we use a Gaussian copula function with beta marginals to model the joint prior distribution of multinomial probabilities. The proposed Gaussian copula prior assumes that the dependence structure between the transformed probabilities can be represented by a multivariate normal distribution, where the marginal prior distribution of each of them is still expressed as a beta distribution. This requires further conditional quartile assessments to describe the correlation structure between these probabilities. The monotonicity of the Gaussian copula transformation allows conditional quartiles of the multinomial probabilities to be transformed into normal quartiles. The latter are used to obtain productmoment correlations for normal variates. This powerful technique of transforming quartiles avoids the difficulties encountered when transforming product-moment correlations. Structural assessments of the conditional quartiles are used to ensure that the encoded correlation matrix is positivedefinite. Although the encoding of a Gaussian copula requires more assessments than the Dirichlet distribution, it gives a more flexible multivariate prior distribution. A practical example where the methods are used was given and the two encoded prior distributions compared.
The two methods have been implemented in a single userfriendly software package that is freely available at http:// statistics.open.ac.uk/elicitation. The user has the option of eliciting either of the two priors or both of them simultaneously. For the Dirichlet prior distribution, the available software outputs the encoded hyperparameter vector a * . The expectations and variances of each p i are also given. For the Gaussian Copula prior, the output includes the encoded pairs of beta parameters α # i and β # i (i = 1, . . . , k − 1), together with the encoded correlation matrix R * . then, from Eq. (17), replacing the expected values by the median function, M(.), we get which, utilizing monotonicity, gives From (20), statement (21) is also valid for j = k. This coincides with the definition of Z k = 1, and completes the proof of the second statement of Lemma 2.
Encoding the correlation matrix R i (i = 2, . . . , k − 1) Let R i be partitioned as follows where R i−1 r i = Cov(y i−1 , Y i ), and σ 2 i = Var(Y i ). Although the Gaussian copula function imposes the constraint that Var(Y i ) = 1, we will find another estimate for σ 2 i using the conditional variance of Y i elicited in (19). The reason for that, as will be shown later, is to follow the approach of Kadane et al. (1980), which will ensure the positive-definiteness of the matrix R i . In what follows, we use the conditional median assessments to estimate r i .
This holds for j = 1, . . . , i − 1, so we have a system of i − 1 equations of the form where u i = (m i,1 , . . . , m i,i−1 ) and Provided that η j = 0 ( j = 1, . . . , i − 1), the upper diagonal matrix Q i−1 is non-singular, and hence r i = Q −1 i−1 u i . The existence of this solution for r i is guaranteed by choosing a non-zero value for η * j ( j = 1, . . . , i − 1) when eliciting conditional medians in Sect. 3.2.3. It can be seen from the relationship in (18) that η j = 0 if and only if η * j = 0 ( j = 1, . . . , i − 1).
Since Var(Y i |y i−1 ) = σ 2 i − r i R i−1 r i , we can use the assessed conditional variance given by V i,i−1 in (19) to estimate the unconditional variance as σ 2 i = V i,i−1 + r i R i−1 r i . Then, using the Schurr complement, the matrix R i is positivedefinite if and only if σ 2 i −r i R i−1 r i > 0, which is guaranteed from (19), since V i,i−1 > 0.

The WinBUGS model
The following is the WinBUGS model for the example in Sect. 5.