Flexible Prior Elicitation via the
Prior Predictive Distribution

The process of performing Bayesian statistical inference usually starts by building a joint probability distribution of observable variables/measurements $\operatorname{\boldsymbol{Y}}$ and unobservable parameters $\operatorname{\boldsymbol{\theta}}$. The corresponding marginal distribution with respect to $\operatorname{\boldsymbol{\theta}}$ is referred to as the prior distribution and the marginal distribution with respect to $\operatorname{\boldsymbol{Y}}$ is referred to as the prior predictive distribution. According to the Bayesian paradigm, the prior distribution should be designed independently of the measurement outcomes, that is to say, it must reflect our prior knowledge about the parameters $\operatorname{\boldsymbol{\theta}}$ before seeing the actual independent measurements $\operatorname{\mathbf{y}}_{1},\operatorname{\mathbf{y}}_{2},\ldots$ (i.e., realizations of $\operatorname{\boldsymbol{Y}}$) obtained in the experiments (Berger, 1993; O’Hagan, 2004). After having obtained the measurements, the posterior distribution of $\operatorname{\boldsymbol{\theta}}$ arises from the joint distribution by conditioning on $\operatorname{\mathbf{y}}_{1}$, $\operatorname{\mathbf{y}}_{2}$, $\ldots$ (O’Hagan, 2004).

Let $\operatorname{\boldsymbol{Y}}=[Y_{1}\ldots Y_{S}]$ be a $S$-dimensional vector of observable variables and denote the sample space $\Omega$ as a subset of $\mathbb{R}^{S}$. Hereafter we denote by $\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}\sim\pi_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}$ our data probability distribution conditioned on the parameters. We also write $\operatorname{\boldsymbol{\theta}}\sim\pi_{\operatorname{\boldsymbol{\theta}}}$ where $\operatorname{\boldsymbol{\theta}}\in\Theta\subseteq\mathbb{R}^{D}$ and $\pi_{\operatorname{\boldsymbol{\theta}}}$ belongs to a given family of parametric distributions, say $\mathcal{F}_{\operatorname{\boldsymbol{\lambda}}}$ indexed by a hyperparameter vector $\operatorname{\boldsymbol{\lambda}}$. Then, by marginalizing out the parameters $\operatorname{\boldsymbol{\theta}}$, the prior predictive distribution is given by

$\displaystyle\pi_{\operatorname{\boldsymbol{Y}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\lambda}})=\displaystyle\int_{\Theta}\pi_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\theta}})\pi_{\operatorname{\boldsymbol{\theta}}}(\operatorname{\boldsymbol{\theta}}|\operatorname{\boldsymbol{\lambda}})\ \mathrm{d}\operatorname{\boldsymbol{\theta}}.$

(1)

The prior predictive distribution is not to be confused with the marginal likelihood of observed data, which is obtained by marginalization over $\operatorname{\boldsymbol{\theta}}$ of the observed data’s sampling distribution times the prior (e.g., Jeffreys and Zellner, 1980).

Given any subset $A\subseteq\Omega$, the prior predictive probability of $A$, denoted as $\mathbb{P}(\operatorname{\boldsymbol{Y}}\in A|\operatorname{\boldsymbol{\lambda}})$, can be obtained by exchanging the order of integration via the Fubini-Tonelli theorem (Folland, 2013) as

	$\displaystyle\mathbb{P}_{A|\operatorname{\boldsymbol{\lambda}}}$	$\displaystyle:=\displaystyle\int_{A}\pi_{\operatorname{\boldsymbol{Y}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\lambda}})\ \mathrm{d}\operatorname{\mathbf{y}}$
		$\displaystyle=\operatorname{\mathbb{E}}_{\operatorname{\boldsymbol{\theta}}}\big{(}\mathbb{P}_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}({\operatorname{\boldsymbol{Y}}\in A|\operatorname{\boldsymbol{\theta}}})\big{)}.$		(2)

See supplementary materials for details. The hyperparameter vector $\operatorname{\boldsymbol{\lambda}}$, which defines a particular prior from the set of all priors $\mathcal{F}_{\operatorname{\boldsymbol{\lambda}}}$, will be treated as constant. Hence, no prior needs to be assigned to it. Instead, the values of $\operatorname{\boldsymbol{\lambda}}$ will be obtained during the prior predictive elicitation method presented below.

Our assumption is that the output elicitation procedure provides information as probabilistic assignments regarding the data vector $\operatorname{\boldsymbol{Y}}$ falling within a fixed set of mutually exclusively and exhaustive events $\mathbf{A}$. Such collection of assignments can be considered as the data available for inferring the prior, and is not to be confused by actual measurement data following the generative model. Our focus here is in the mathematical machinery required for converting this information into prior distributions, not taking any stance on how the information is collected from the expert. However, we will briefly discuss the elicitation process itself in Section 3.4.

Let $\mathbf{A}$ $=$ $\{A_{1},\ldots,A_{n}\}$ be a partition of the sample space $\Omega$. Throughout the elicitation procedure, the expert supplies their expected opinions regarding the quantities $\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}$ for all $i=1,\ldots,n$. The expert’s judgements themselves are not fully deterministic and retain some uncertainty. Also, the expert may be more comfortable to make statements for certain partitions of $\Omega$ than for others.

To account for the uncertainty in the probability quantifications of $\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}$, we assume that the obtained judgements $\operatorname{\mathbf{p}}$ follow a Dirichlet distribution (Ferguson, 1973) with base measure given by the prior predictive probabilities $\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}$ and precision parameter $\alpha$. Hence, for any chosen partition $\mathbf{A}$ of size $n$, we denote the distribution of $\operatorname{\mathbf{p}}$ as

$\displaystyle\operatorname{\mathbf{p}}|\alpha,\operatorname{\boldsymbol{\lambda}}\sim\mathcal{D}(\alpha,[\mathbb{P}_{A_{1}|\operatorname{\boldsymbol{\lambda}}}\cdots\mathbb{P}_{A_{n}|\operatorname{\boldsymbol{\lambda}}}]),$

(3)

where $\mathcal{D}(\cdot)$ stands for Dirichlet distribution and whose multivariate density function reads

$\displaystyle\mathcal{D}(\operatorname{\mathbf{p}}|\alpha,\operatorname{\boldsymbol{\lambda}})$

$\displaystyle=\dfrac{\Gamma(\alpha)}{\prod_{i=1}^{n}\Gamma(\alpha\hskip 2.27626pt\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}})}\prod_{i=1}^{n}p_{i}^{\alpha\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}-1}.$

(4)

Naturally, we require $\sum_{i=1}^{n}\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}=1$. The Dirichlet density (4) accounts for the uncertainty inherent to the numerical quantification of the probability vector $\operatorname{\mathbf{p}}$ due to, for example, biases introduced through the mechanisms of elicitation processes (the way in which questions are made), practical imperfection (imprecision) of experts’ judgements in probabilistic terms or poor judgements on the effect of parameters in the output model. For details and in-depth discussion, see O’Hagan and Oakley (2004), O’Hagan (2019) and Sarma and Kay (2020) .

The hyperparameter $\alpha$ measures how well the prior predictive probability model is able to represent (or reproduce) the probability data provided in the elicitation process. The larger the values of $\alpha$, the less variance around the expected value $\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}$. For practical use of this principle, we can find the maximum likelihood estimate (MLE) $\hat{\alpha}$ of $\alpha$, which can be directly understood in terms of the deviance between the prior predictive probability and the experts opinion. More specifically, we have

$\displaystyle\hat{\alpha}\approx\dfrac{n/2-1/2}{{\rm KL}(\operatorname{\mathbb{P}_{\boldsymbol{A}|\operatorname{\boldsymbol{\lambda}}}}\ ||\ \operatorname{\mathbf{p}})}$

(5)

where $\operatorname{\mathbb{P}_{\boldsymbol{A}|\operatorname{\boldsymbol{\lambda}}}}=[\mathbb{P}_{A_{i}|\operatorname{\boldsymbol{\lambda}}}\cdots\mathbb{P}_{A_{n}|\operatorname{\boldsymbol{\lambda}}}]^{\top}$ and ${\rm KL}(\operatorname{\mathbb{P}_{\boldsymbol{A}|\operatorname{\boldsymbol{\lambda}}}}\ ||\ \operatorname{\mathbf{p}})$ is the Kullback-Leibler divergence between the two distributions. The practical interpretation is that for small KL values, we would not be able discriminate the prior predictive probability from the probability data provided by the expert. See supplementary materials for the proof of Equation (5).

Even though we work in a Bayesian context looking to recover a prior distribution, the core procedure of our method applies classical statistical inference. Given a numerical vector of probabilities from the elicitation process, the goal is to show that we are able to find the value of certain parameters (in this case the hyperparameters $\operatorname{\boldsymbol{\lambda}}$ and concentration $\alpha$ parameter) of the Dirichlet probabilistic model (3) which would have most likely generated this particular data (of user’s subjective knowledge). In other words, we are aiming to obtain the maximum likelihood estimator (MLE).

To study the MLE, we consider the limit where the partitioning is made increasingly more fine grained by increasing $n$ towards infinity. However, we still only obtain information from the user once (i.e., for a single partitioning). That is, the user is providing more and more information about the probabilities, but does not repeat the procedure multiple times. As we will show below, the MLE is consistent under these circumstances, providing the true $\operatorname{\boldsymbol{\lambda}}$ when $n\rightarrow\infty$, under reasonable assumptions.

Recall that equations (3) and (4) represent the probabilistic model of $\operatorname{\mathbf{p}}$ conditioned on the parameters $\operatorname{\boldsymbol{\eta}}=(\operatorname{\boldsymbol{\lambda}},\alpha)$. Suppose the implied true prior distribution of the expert has hyperparameter values $\operatorname{\boldsymbol{\lambda}}_{0}$ and denote $\operatorname{\boldsymbol{\eta}}_{0}=(\operatorname{\boldsymbol{\lambda}}_{0},\alpha_{0})$. Take the size of the partition $n$ to be large and denote the log-likelihood as $T_{\operatorname{\boldsymbol{\eta}}}(\operatorname{\mathbf{p}})=\log\mathcal{D}(\operatorname{\mathbf{p}}|\alpha,\operatorname{\boldsymbol{\lambda}})$ with expectation $Q_{\operatorname{\boldsymbol{\eta}}_{0}}(\operatorname{\boldsymbol{\eta}})=\operatorname{\mathbb{E}}_{\mathcal{D}}(T_{\operatorname{\boldsymbol{\eta}}}(\operatorname{\mathbf{p}}))$.

We show that the expected log-likelihood is maximized at $\operatorname{\boldsymbol{\eta}}_{0}$. By Jensen’s inequality, we know that

$\displaystyle\operatorname{\mathbb{E}}_{\mathcal{D}}$

$\displaystyle\left[-\log\dfrac{\mathcal{D}(\operatorname{\mathbf{p}}|\alpha,\operatorname{\boldsymbol{\lambda}})}{\mathcal{D}(\operatorname{\mathbf{p}}|\alpha_{0},\operatorname{\boldsymbol{\lambda}}_{0})}\right]>-\log\operatorname{\mathbb{E}}\left[\dfrac{\mathcal{D}(\operatorname{\mathbf{p}}|\alpha,\operatorname{\boldsymbol{\lambda}})}{\mathcal{D}(\operatorname{\mathbf{p}}|\alpha_{0},\operatorname{\boldsymbol{\lambda}}_{0})}\right]=0,$

(6)

yielding

$\displaystyle Q_{\operatorname{\boldsymbol{\eta}}_{0}}(\operatorname{\boldsymbol{\eta}}_{0})=\operatorname{\mathbb{E}}_{\mathcal{D}}(T_{\operatorname{\boldsymbol{\eta}}_{0}}(\operatorname{\mathbf{p}}))>\operatorname{\mathbb{E}}_{\mathcal{D}}(T_{\operatorname{\boldsymbol{\eta}}}(\operatorname{\mathbf{p}}))=Q_{\operatorname{\boldsymbol{\eta}}_{0}}(\operatorname{\boldsymbol{\eta}}),$

which holds for all $\operatorname{\boldsymbol{\eta}}$. The expectation $\operatorname{\mathbb{E}}_{\mathcal{D}}(\cdot)$ is taken with respect to the distribution (4). The technical condition to ensure uniqueness of the MLE is that the probabilistic model (4) must be identifiable¹¹1In practise, this may not be an issue when fitting the model. However, we believe it is important to understand the theoretical properties of the inference process so that we can avoid problems in the optimisation procedures.. That is, equality of likelihoods must imply equality of parameters: $\mathcal{D}(\operatorname{\mathbf{p}}|\alpha_{1},\operatorname{\boldsymbol{\lambda}}_{1})$ $=$ $\mathcal{D}(\operatorname{\mathbf{p}}|\alpha_{2},\operatorname{\boldsymbol{\lambda}}_{2})\Rightarrow\operatorname{\boldsymbol{\eta}}_{1}=\operatorname{\boldsymbol{\eta}}_{2}$ for all $\operatorname{\mathbf{p}}$. Otherwise we may encounter multiple maxima and thus the prior distribution in the set $\mathcal{F}_{\operatorname{\boldsymbol{\lambda}}}$ is not unique.

Next, we demonstrate how the proposed approach can be used for concrete modelling problems, by detailing the procedure for the widely-used family of generalized linear models (GLM; Nelder and Wedderburn, 1972). As GLMs typically have several parameters – one parameter per predicting covariate plus an intercept and potentially a dispersion parameter – direct specification of the parameters’ joint prior is often difficult. However, our prior predictive approach can handle this situation elegantly.

In case of a GLM, our elicitation method requires the selection of sets of covariate values for which the expert is comfortable to express probability judgements about plausible realizations of $\operatorname{\boldsymbol{Y}}$. More formally, for each set of covariates $\operatorname{\mathbf{x}}_{j}$ $=$ $[x_{j,1}\cdots x_{j,C}]$, $j=1,\ldots,J$, the expert provides probability judgements $\operatorname{\mathbf{p}}_{j}=[p_{j,1}\cdots p_{j,n_{j}}]$ with $\sum_{i_{j}=1}^{n_{j}}p_{j,i_{j}}=1$, where $n_{j}$ is the partition size for covariate set $j$ implying the partition $\mathbf{A}_{j}=\{A_{j,1},$ $\ldots$, $A_{j,n_{j}}\}$. Under the assumption of the judgement $\operatorname{\mathbf{p}}_{j}$ being pairwise conditionally independent, we can express the likelihood function of $\alpha$ and $\operatorname{\boldsymbol{\lambda}}$ as

$\displaystyle\mathcal{D}(\operatorname{\mathbf{p}}_{1},\ldots,\operatorname{\mathbf{p}}_{J}|$

$\displaystyle\alpha,\operatorname{\boldsymbol{\lambda}})=\dfrac{\Gamma(\alpha)^{J}}{\prod\limits_{j=1}^{J}\prod\limits_{i_{j}=1}^{n_{j}}\Gamma(\alpha\hskip 2.27626pt\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}},\operatorname{\mathbf{x}}_{j}})}\prod_{j=1}^{J}\prod_{i_{j}=1}^{n_{j}}p_{j,i_{j}}^{\alpha\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}},\operatorname{\mathbf{x}}_{j}}-1}$

(7)

where $\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}},\operatorname{\mathbf{x}}_{j}}$ is the prior predictive probability for the set $A_{j,i_{j}}$ related to covariate set $\operatorname{\mathbf{x}}_{j}$.

Importantly, there is no need for the partitions themselves or their size to be the same throughout the sets of covariate values: For each $j$, the expert can create any partition they are most comfortable with making judgements about. This feature provides much more freedom to the expert in expressing their knowledge of the data compared to alternative methods. For example, to obtain a prior distribution for logistic regression model, the method of Bedrick et al. (1997) requires the user to provide a fixed number of probabilities just enough to make the Jacobians appearing in their method invertible.

Using the machinery above requires obtaining the probability judgements $\operatorname{\mathbf{p}}$ from the user. The method itself is general, and can be used as part of any practical Bayesian modelling workflow when linked to any particular elicitation interface. We have implemented an extension of the SHELF interface (Oakley and O’Hagan, 2019) as a reference, by replacing the direct parameter elicitation components with variants that query the user for the prior predictive probabilities. This readily provides practical elicitation methods for the user to specify probabilities by utilizing probability quantiles or roulette chips. This means that probability ratios for events are provided and then individual probabilities are recovered under the natural constraint $\sum_{i=1}^{n}p_{i}=1$. Hence, the user can choose the way of providing information they feel most comfortable with. Besides graphical interfaces, the elicitation can be carried out by the modeller interviewing a domain expert. Experienced modellers may also choose to simply express some particular priors via providing $\operatorname{\mathbf{p}}$ while designing the model.

In this section we highlight the steps to obtain the prior predictive probability, by rewriting it as a expected value w.r.t. the prior distribution as follows. Given that the probabilistic models, $\pi_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\theta}})$ and the prior $\pi_{\operatorname{\boldsymbol{\theta}}}$ are positive functions, we can rearrange the order of the integration (See Folland, 2013, Fubini-Tonelli theorem). Hence we have

	$\displaystyle\mathbb{P}_{A|\operatorname{\boldsymbol{\lambda}}}$	$\displaystyle:=\displaystyle\int_{A}\pi_{\operatorname{\boldsymbol{Y}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\lambda}})\ \mathrm{d}\operatorname{\mathbf{y}}=\displaystyle\int_{A}\displaystyle\int_{\Theta}\pi_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\theta}})\pi(\operatorname{\boldsymbol{\theta}}|\operatorname{\boldsymbol{\lambda}})\ \mathrm{d}\operatorname{\boldsymbol{\theta}}\mathrm{d}\operatorname{\mathbf{y}}$
		$\displaystyle\stackrel{{\scriptstyle\textrm{Fub.}}}{{=}}\displaystyle\int_{\Theta}\displaystyle\int_{A}\pi_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\theta}})\pi(\operatorname{\boldsymbol{\theta}}|\operatorname{\boldsymbol{\lambda}})\ \mathrm{d}\operatorname{\mathbf{y}}\mathrm{d}\operatorname{\boldsymbol{\theta}}=\displaystyle\int_{\Theta}\mathbb{P}_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}({\operatorname{\boldsymbol{Y}}\in A|\operatorname{\boldsymbol{\theta}}})\pi(\operatorname{\boldsymbol{\theta}}|\operatorname{\boldsymbol{\lambda}})\ \mathrm{d}\operatorname{\boldsymbol{\theta}}$
		$\displaystyle=\operatorname{\mathbb{E}}_{\operatorname{\boldsymbol{\theta}}}\big{(}\mathbb{P}_{\operatorname{\boldsymbol{Y}}|\operatorname{\boldsymbol{\theta}}}({\operatorname{\boldsymbol{Y}}\in A|\operatorname{\boldsymbol{\theta}}})\big{)}.$		(14)

Here we show the approximate behaviour of the precision parameter $\alpha$ for the general case when covariates are present. The simplification to other cases in straightforward. Recall the likelihood function of $\operatorname{\boldsymbol{\lambda}}$ given expert data reads,

$\displaystyle\mathcal{D}$

$\displaystyle(\operatorname{\mathbf{p}}_{1},\ldots,\operatorname{\mathbf{p}}_{J}|\alpha,\operatorname{\boldsymbol{\lambda}})=\dfrac{\Gamma(\alpha)^{J}}{\prod\limits_{j=1}^{J}\prod\limits_{i_{j}=1}^{n_{j}}\Gamma(\alpha\hskip 2.27626pt\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}})}\prod_{j=1}^{J}\prod_{i_{j}=1}^{n_{j}}p_{j,i_{j}}^{\alpha\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}-1}.$

(15)

Consider the Stirling’s approximation²²2This is a precise approximation. to the $\Gamma(\cdot)$ function given by,

$\displaystyle\Gamma(x)\approx\sqrt{\dfrac{2\pi}{x}}\left(\dfrac{x}{e}\right)^{x}.$

(16)

Rewriting the likelihood function in terms of the above approximation and removing terms that does not depend on $\alpha$ with a simplified notation we get,

	$\displaystyle\mathcal{D}(\operatorname{\mathbf{p}}|\alpha,\operatorname{\boldsymbol{\lambda}})$	$\displaystyle\approx\dfrac{\left(\sqrt{\dfrac{2\pi}{\alpha}}\left(\dfrac{\alpha}{e}\right)^{\alpha}\right)^{J}}{\prod\limits_{j,i_{j}}\sqrt{\dfrac{2\pi}{\alpha\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}}}\left(\dfrac{\alpha\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}}{e}\right)^{\alpha\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}}}\exp\left(\sum_{i,j}\alpha(\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}-1)\log p_{j,i_{j}}\right)$
		$\displaystyle\approx\dfrac{\alpha^{{}^{\sum_{j}n_{j}/2-J/2}}\prod\limits_{i,j}\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}^{1/2}}{\exp\left(\alpha\sum\limits_{i,j}\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}\log\ \dfrac{\mathbb{P}_{A_{j,i_{j}}|\operatorname{\boldsymbol{\lambda}}}}{p_{i,i_{j}}}\right)}$		(17)

Take the logarithm of the above function and the derivative w.r.t. $\alpha$. Setting it to zero and solving for $\alpha$ we obtain,

$\displaystyle\hat{\alpha}\approx\dfrac{\sum_{j}n_{j}/2-J/2}{\sum\limits_{j}KL(\mathbb{P}_{{}_{j}}||\operatorname{\mathbf{p}}_{j})}$

(18)

where the notation $\mathbb{P}_{{}_{j}}=[\mathbb{P}_{A_{j,1}|\operatorname{\boldsymbol{\lambda}}}\cdots\mathbb{P}_{A_{j,n_{j}}|\operatorname{\boldsymbol{\lambda}}}]^{\top}$ and $KL(P||Q)$ denotes the Kullback-Leibler divergence in this order.

The Fisher information matrix for the unknown hyperparameters can be obtained in closed-form by the fact that, in the original parametrisation of the Dirichlet distribution, the Fisher information is already known. In the original parametrisation and in its basic form, the probability density function reads

$\displaystyle\mathcal{D}(\operatorname{\mathbf{p}}|\alpha,\mathbb{P})$

$\displaystyle=\dfrac{\Gamma(\alpha)}{\prod_{i=1}^{n}\Gamma(\alpha\hskip 2.27626pt\mathbb{P}_{i})}\prod_{i=1}^{n}p_{i}^{\alpha\mathbb{P}_{i}-1}$

(19)

where $\mathbb{P}=[\mathbb{P}_{1}\cdots\mathbb{P}_{n}]^{\top}$. Also knowing that the Dirichlet distribution belongs the exponential family, the Fisher information matrix reads,

$\displaystyle H_{\mathbb{P}}=\alpha^{2}(\operatorname{\mathrm{diag}}(\psi^{\prime}(\alpha\mathbb{P}))-\psi^{\prime}(\alpha)\mathds{1}\mathds{1}^{\top}),$

(20)

whose inverse is given in closed-form as

$\displaystyle H_{\mathbb{P}}^{-1}=\tfrac{1}{\alpha^{2}}\bigg{(}\operatorname{\mathrm{diag}}(\psi^{\prime}(\alpha\mathbb{P}))^{-1}+\frac{\operatorname{\mathrm{diag}}(\psi^{\prime}(\alpha\mathbb{P}))^{-1}\mathds{1}\mathds{1}^{\top}\operatorname{\mathrm{diag}}(\psi^{\prime}(\alpha\mathbb{P}))^{-1}}{(1/\psi^{\prime}(\alpha)-\mathds{1}^{\top}\operatorname{\mathrm{diag}}(\psi^{\prime}(\alpha\mathbb{P}))^{-1}\mathds{1})}\bigg{)}$

(21)

where $\mathds{1}$ is $n\times 1$ vector with each component equals to $1$.

In the main paper, the vector of parameters $\mathbb{P}$ of the Dirichlet distribution is written as a function of $\operatorname{\boldsymbol{\lambda}}$. Using the change of variables for a new parametrisation (see Calderhead, 2012; Girolami and Calderhead, 2011, page 64, Section 3.2.5, equation 3.27), the Fisher information matrix with respect to $\operatorname{\boldsymbol{\lambda}}$ can be obtained directly (by passing any need of recalculating integrals) as,

$\displaystyle H_{\operatorname{\boldsymbol{\lambda}}}=\big{[}\tfrac{d}{d\lambda_{1}}\mathbb{P}\cdots\tfrac{d}{d\lambda_{M}}\mathbb{P}\big{]}^{\top}H_{\mathbb{P}}\big{[}\tfrac{d}{d\lambda_{1}}\mathbb{P}\cdots\tfrac{d}{d\lambda_{M}}\mathbb{P}\big{]}$

(22)

where the vector $\tfrac{d}{d\lambda_{m}}\mathbb{P}=\big{[}\tfrac{d}{d\lambda_{m}}\mathbb{P}_{1}\cdots\tfrac{d}{d\lambda_{m}}\mathbb{P}_{n}\big{]}^{\top}$ (the Jacobian matrix). Note that $H_{\mathbb{P}}$ is invertible and positive-definide, so as $H_{\operatorname{\boldsymbol{\lambda}}}$. Hence $H_{\operatorname{\boldsymbol{\lambda}}}$ is also invertible and its cholesky decomposition is stable to compute.

For the case where $\mathbb{P}_{j}$ does not have closed-form expression we can estimate $\mathbb{P}_{j}$ and its derivatives w.r.t $\operatorname{\boldsymbol{\lambda}}$ using the reparametrisation gradients and automatic differentiation. The main idea is to find a pivotal function (see Casella and Berger, 2001, page 427, Section 9.2.2) and obtain Monte-Carlo estimates of $\mathbb{P}_{j}$ and gradients $d/d\lambda_{m}\mathbb{P}_{j}$ with low computational cost according to Figurnov et al. (2018) and Mohamed et al. (2019).

With a simplified notation, recall the prior distribution $\pi_{\operatorname{\boldsymbol{\theta}}}$ and that the prior predictive probability can be rewritten as a expected value

$\displaystyle\mathbb{P}_{A|\operatorname{\boldsymbol{\lambda}}}=\operatorname{\mathbb{E}}_{\operatorname{\boldsymbol{\theta}}}\big{(}\mathbb{P}({\operatorname{\boldsymbol{Y}}\in A|\operatorname{\boldsymbol{\theta}}})\big{)}$

(24)

which depends on $\operatorname{\boldsymbol{\lambda}}$. Here the expression $\mathbb{P}({\operatorname{\boldsymbol{Y}}\in A|\operatorname{\boldsymbol{\theta}}})$ depends only on $\operatorname{\boldsymbol{\theta}}$. Then, find a pivotal function $X=T(\operatorname{\boldsymbol{\theta}})$ such that the distribution of $X$ does not depend on $\operatorname{\boldsymbol{\lambda}}$. We then can rewrite the expectation,

$\displaystyle\mathbb{P}_{A|\operatorname{\boldsymbol{\lambda}}}=\operatorname{\mathbb{E}}_{X}\big{(}\mathbb{P}({\operatorname{\boldsymbol{Y}}\in A|T^{-1}_{X}(\operatorname{\boldsymbol{\lambda}}}))\big{)}$

(25)

The gradients can be computed interchanging the order of integration and derivation,

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}\lambda_{m}}\mathbb{P}_{A|\operatorname{\boldsymbol{\lambda}}}=\operatorname{\mathbb{E}}_{X}\left(\dfrac{\mathrm{d}}{\mathrm{d}\lambda_{m}}\mathbb{P}({\operatorname{\boldsymbol{Y}}\in A|T^{-1}_{X}(\operatorname{\boldsymbol{\lambda}}}))\right).$

(26)

Where $T^{-1}_{X}(\cdot)$ is a inverse function of $T$ and depends on $X$ and $\operatorname{\boldsymbol{\lambda}}$. The important notion here is that there is no need for resampling $X$ since the distribution $\pi_{X}(\cdot)$ is free of $\operatorname{\boldsymbol{\lambda}}$ by definition.

The probabilistic model for observed data (stature of male human being) is specified as follows,

	$\displaystyle Y_{t}|\operatorname{\boldsymbol{\theta}},b$	$\displaystyle\sim\mathcal{W}(h(t;\operatorname{\boldsymbol{\theta}}),b)$
	$\displaystyle b$	$\displaystyle\sim\mathcal{G}(a_{0},b_{0})$
	$\displaystyle\theta_{d}$	$\displaystyle\stackrel{{\scriptstyle i.i.d}}{{\sim}}\mathcal{LN}(a_{d},b_{d})$		(36)

where $Y_{t}$ is univariate $S=1$ and denotes the stature of the human being at time $t$. The parameters of the growth-model $h(t;\operatorname{\boldsymbol{\theta}})$ are denoted as $\operatorname{\boldsymbol{\theta}}=[\theta_{1}\ \theta_{2}\ \theta_{3}\ \theta_{4}\ \theta_{5}]=[h_{1},h_{t_{*}},t_{*},s_{0},s_{1}]^{\top}$, where $h_{1}$ is the average height of an adult human, $h_{t_{*}}$ is the average high for the event ”growth-spurt” (Preece and Baines, 1978), $t_{*}$ is when that event happens, $s_{0}$ and $s_{1}$ are constants from the model. The parameter $b$ controls the variance of the variable $Y_{t}$ around $h(t;\operatorname{\boldsymbol{\theta}})$. Large the values of $b$ less variance around the $h(t;\operatorname{\boldsymbol{\theta}})$ and vice-versa. $\mathcal{W}$, $\mathcal{G}$ and $\mathcal{LN}$ stands for respectively, Weibull, Gamma and log-Normal distributions.

We used the Weibull distribution in the mean-variance parametrisation which means that the probability distribution of $Y_{t}|\operatorname{\boldsymbol{\theta}},b$ is given by,

$\displaystyle\pi_{Y_{t}|\operatorname{\boldsymbol{\theta}},b}(y)=b\hskip 1.42271pt\tfrac{\Gamma(1+1/b)}{\exp\left(h(t;\operatorname{\boldsymbol{\theta}})\right)}\hskip 1.42271pt\left(y\hskip 1.42271pt\tfrac{\Gamma(1+1/b)}{\exp\left(h(t;\operatorname{\boldsymbol{\theta}})\right)}\right)^{b-1}\exp\left(-\left(y\hskip 1.42271pt\tfrac{\Gamma(1+1/b)}{\exp\left(h(t;\operatorname{\boldsymbol{\theta}})\right)}\right)^{b}\right)$

(37)

The other distribution used for the prior are used in their standard parametrisation scale-shape for Gamma and mean-variance for log-Normal distribution. The vector of hyperparameters is $\operatorname{\boldsymbol{\lambda}}=\{a_{m},b_{m},m=0,\ldots,5\}$. The human-growth model obtained by Preece and Baines (1978) and given in Section 2, Model 1 in their paper. In our notation this growth-model reads

$\displaystyle h(t;\operatorname{\boldsymbol{\theta}})=h_{1}-\dfrac{2(h_{1}-h_{t_{*}})}{\exp[s_{0}\ (t-t_{*})]+\exp[s_{1}\ (t-t_{*})]}.$

(38)