An Objective Prior from a Scoring Rule

Abstract

In this paper, we introduce a novel objective prior distribution levering on the connections between information, divergence and scoring rules. In particular, we do so from the starting point of convex functions representing information in density functions. This provides a natural route to proper local scoring rules using Bregman divergence. Specifically, we determine the prior which solves setting the score function to be a constant. Although in itself this provides motivation for an objective prior, the prior also minimizes a corresponding information criterion.

1. Introduction

A major drawback of objective priors, such as Jeffreys prior [1] and the reference prior [2], is that in many cases, they are improper. The case for objective priors has been made, see for example [3], yet while for a parameter that is defined over a bounded interval, such as $(0,1)$, it is generally possible to derive objective prior distributions that are proper, this is not the case for parameters on $(0,\infty )$ or $(-\infty ,\infty )$. The literature provides many examples where improper prior distributions cannot be suitably employed; such as Bayes factors, mixture models and hierarchical models, to name but a few. Methods have been proposed to overcome these obstacles, for example, Intrinsic Bayes Factors [4] and Fractional Bayes factors [5] or reparametrizing mixture models [6]. However, these types of results are generally valid for a limited number of specific conditions. Additionally, improper prior distributions are not too suitable to be employed where large numbers of parameters are involved as it would be difficult to establish properness of the full posterior distribution.

The idea of this paper is to present a novel objective prior distribution for continuous parameter spaces by considering the connection between information, divergence and scoring rules. In particular, the proposed prior can be defined over $(0,\infty )$ and $(-\infty ,\infty )$, the latter by extending the former, and it has the appealing property of being proper.

Recently, Ref. [7], introduced a new class of objective prior which solved a differential equation of the form $S(q,q^{\prime },q^{″})=0$, where S is a score function and the solution q acts as the prior distribution, and where $q^{\prime }$ and $q^{″}$ are, respectively, the first and second derivative of q. The solution is also shown to minimize an information criterion.

There are two well-known relations that connect information, proper local scoring rules and divergences. The most famous of which links Shannon information, Kullback–Leibler divergence and the log–score, given by

$$\int plog\frac{p}{q}=\int plogp+\int p(-logq),$$

(1)

where p and q are two densities, and integrals will be generally defined with respect to the Lebesgue measure. The term on the left-hand-side of (1) is the Kullback–Leibler divergence [8] between p and q, the first term on the right-hand-side is the Shannon information associated with density p, and the second term is the expectation of the log-score function.

Another way to connect information, divergence and proper local scoring rules, involves Fisher divergence, Fisher information, and the Hyvärinen score function ([9]):

$$\int p{\left(\frac{p^{\prime }}{p}-\frac{q^{\prime }}{q}\right)}^2=\int \frac{{\left(p^{\prime }\right)}^2}{p}+\int p\left(2\frac{q^{″}}{q}-{\left(\frac{q^{\prime }}{q}\right)}^2\right),$$

where the final term has been obtained using an integration by parts, requiring $p{q}^{\prime }/q$ to vanish at the boundary values. In general, these relationships can be expressed as

$$D(p,q)=I(p)+\int pS(q),$$

(2)

where D denotes the divergence, I the measure of information and S the score, and clearly $I\left(p\right)=-\int pS\left(p\right)$.

Recently, in [10], a new class of score function was introduced, where the starting point is the property of the score function, which is that

$$p=arg\underset{q}{min}\int pS\left(q\right),$$

for all densities p. In other words, a score is said to be proper if the above is minimized by the choice of $q=p$. Let us consider the well-known log-score, $S\left(q\right)=-logq\left(x\right)$. Then, we have that it satisfies the above property, since for any density p it is that $\int plog(p/q)\ge 0$, with equality only when $q\equiv p$. As such, we have that the log-score is a proper score. Furthermore, a score is said to be local if it only depends on q through the density value $q\left(x\right)$. See [10,11]. It must be noted that the log-score is the only proper score to be local.

If we consider the Hyvärinen score function [9], which is given by

$$S(q,q^{\prime },q^{″})=2\frac{q^{″}}{q}-{\left(\frac{q^{\prime }}{q}\right)}^2,$$

which we note depends on $(q,q^{\prime },q^{″})$, i.e., q and the first two derivatives, as such it is not local in the above sense. However, the locality condition can be weakened [10] by allowing the score to depend on a finite number m of derivatives. Therefore, the Hyvärinen score will be an order–2 proper local scoring rule.

More generally, if a proper score depends on m derivatives, then it will be defined an order–m local scoring rule. The theory in support of this is based on the fact that the minimizer of $\int pS\left(q\right)$ is p, and this can be investigated using variational analysis. The relevant Euler–Lagrange equation of order two being

$$S+q\frac{\partial S}{\partial q}-\frac{d}{dx}q\frac{\partial S}{\partial q^{\prime }}+\frac{d^2}{d{x}^2}q\frac{\partial S}{\partial q^{″}}=0.$$

(3)

The corresponding general case of (3) is given as Equation (18) in [10]. Throughout this paper we will focus on the case $m=2$, since this is where we draw our prior from. The Appendix A provides the expression for a general m.

In [10], the solution to Equation (3), is proposed using properties of differential operators and 1–homogeneous functions. Recall that a 1–homogeneous function f is such that $f(x,\lambda q,\lambda q^{\prime })=\lambda f(x,q,q^{\prime })$ for any $\lambda >0$. In particular, the Hyvärinen score arises with $f(x,q,q^{\prime })={\left(q^{\prime }\right)}^2/q$ and

$$S\left(q\right)=-\frac{\partial f}{\partial q}+\frac{d}{dx}\frac{\partial f}{\partial q^{\prime }}.$$

(4)

Furthermore, Refs. [10,11] characterize all local and proper scoring rules of order $m=2$. With this respect, as an additional interesting result, in the Appendix A we present the characterization using measures of information and the Bregman divergence [12]. The benefits of the proposed approach are that complicated mathematical analysis is avoided and the derivation of the local rule is made explicit.

Following [10,11] and the novel derivation of their results using Bregman divergence, which is the focus of the Appendix A, information, divergence and scores can be obtained as follows: For some convex function $\alpha :\mathbb{R}\to \mathbb{R}$,

• Divergence: Given the result (A10), we obtain

  $$D(p,q)=\int p\alpha (p^{\prime }/p)-\int p\frac{\partial \varphi }{\partial q}-\int p^{\prime }\frac{\partial \varphi }{\partial q^{\prime }},$$

  where

  $$\frac{\partial \varphi }{\partial q}=\alpha (q^{\prime }/q)-(q^{\prime }/q){\alpha }^{\prime }(q^{\prime }/q)\mathrm{and}\frac{\partial \varphi }{\partial q^{\prime }}={\alpha }^{\prime }(q^{\prime }/q).$$
  Using integration by parts on the right most integral, and assuming that $[p·\partial \varphi /\partial q^{\prime }]$ vanishes at the extremes of the integral,

  $$D(p,q)=\int p\alpha (p^{\prime }/p)+\int p\left\{\frac{d}{dx}{\alpha }^{\prime }(q^{\prime }/q)-\alpha (q^{\prime }/q)+(q^{\prime }/q){\alpha }^{\prime }(q^{\prime }/q)\right\}.$$
• Information: This follows from the divergence, and from (A4), and is given by

  $$I\left(p\right)=\int p\alpha (p^{\prime }/p).$$
• Score: Again, from the form of the divergence and (A8), this is given by

  $$S(q,q^{\prime },q^{″})=\frac{d}{dx}{\alpha }^{\prime }(q^{\prime }/q)-\alpha (q^{\prime }/q)+(q^{\prime }/q){\alpha }^{\prime }(q^{\prime }/q).$$

  (5)

The score $S(q,q^{\prime },q^{″})$ in (5) generalizes the Hyvärinen score, which arises when $\alpha \left(u\right)=u^2$.

The paper is organized as follows. Section 2 introduces the proposed objective prior. Section 3 includes a thorough simulation study, and an application to mixture models that involves both simulated and real data. In Section 4 we have discussed another critical scenario where improper priors may resolve in improper posteriors, i.e., assigning an objective prior to the variance parameter in a hierarchical model. The supporting theory is presented in the Appendix A. In Appendix A.1 we use Bregman divergence to obtain general forms for score functions and associated divergences and following on from this in Appendix A.2 we detail how we use Bregman divergences to obtain a divergence between probability density functions using their first derivatives, and show how to obtain score functions from these divergences. Appendix A.3 provides the general case using m derivatives. Finally, in Appendix A.4 we make the connection with our derivations of scores and that of [10].

2. New Objective Prior

Ref. [7] proposed constructing objective prior distributions on parameter spaces by solving equations of the kind $S\left(q\right)=0$. Specifically, they used a weighted mixture of the log-score and the Hyvärinen score functions. Please note that the sole use of the log-score function would result in the uniform prior, which is not appropriate in many cases and may yield improper posterior distributions. On the other hand, a weighted combination of the two score functions yields a differential equation given by

$$-wlogq\left(x\right)+\frac{q^{″}\left(x\right)}{q\left(x\right)}-\frac{1}{2}{\left(\frac{q^{\prime }\left(x\right)}{q\left(x\right)}\right)}^2=0,$$

(6)

where q denotes the prior density and $w>0$ the weight balancing the two score functions.

Solutions to the differential Equation (6) can be found for different spaces, and constraints on the shape of q can be considered; so, to have a prior density with desirable behavior, such as monotone, convex, log–concave and more.

We have already seen that the Hyvärinen score arises with $\alpha \left(u\right)=u^2$; see (5). An important property an objective prior distribution may be required to have is a heavy tail and it is this type of prior we are seeking to obtain in a formal way through an appropriate choice of $\alpha$. We will consider such on $(0,\infty )$. Mirroring the Hyvärinen score, we adopt $\alpha \left(u\right)=u^{-2}$ with $u=-q^{\prime }/q$, and q a decreasing density on $(0,\infty )$. In this case, Equation (5) becomes $6u^{\prime }/u^4-3/u^2,$ which, by setting to 0, becomes $u^{\prime }=\frac{1}{2}{u}^2.$ The solution is easily seen to be $u\left(x\right)=-2/(a+x),$ for some constant a. In this case, the prior on the parameter space $(0,\infty )$ is

$$q\left(x\right)=\frac{a}{{(a+x)}^2}.$$

(7)

To obtain a density function on $(0,\infty )$; i.e., q is non-negative, we choose $a>0$, as we are permitted to do through the constant of integration from the differential equation.

Interestingly, the prior in (7), is a Lomax distribution [13] with scale parameter a and shape parameter equal to 1. Recalling that the Lomax distribution can be directly connected to the Pareto Type I and Pareto Type II distributions, its heavy-tailed nature is immediately obvious.

Figure 1 shows the prior with $a=1$.

Making the connection more directly with the theory set out in the paper with $\alpha \left(u\right)=u^{-2}$, we have $\varphi (u,v)=u^3/v^2$ which is easy to show satisfies $\varphi =u\partial \varphi /\partial u+v\partial \varphi /\partial v$. Then using (A8) we obtain

$$S(q,q^{\prime },q^{″})=3{\left(\frac{q}{q^{\prime }}\right)}^2\left\{2\frac{q{q}^{″}}{{\left(q^{\prime }\right)}^2}-3\right\}.$$

(8)

Setting this to zero; i.e., $2q{q}^{″}=3{\left(q^{\prime }\right)}^2$, this can be solved and the solution is precisely of the form $a/{(a+x)}^2$. We now write this all out in a theorem.

Theorem 1.

Let $\varphi (p,p^{\prime })=p^3/{\left(p^{\prime }\right)}^2$ be the convex function appearing in (A3); i.e., $\varphi (u,v)=u\alpha (v/u)$ with $\alpha \left(\xi \right)={\xi }^{-2}$. Then ϕ is convex for either $\xi <0$ or $\xi >0$. The Euler equation associated with this ϕ; i.e., $d/dx\partial \varphi /\partial p^{\prime }=\partial \varphi /\partial p$, yields

$$6p^3{p}^{″}/{\left(p^{\prime }\right)}^4-6p^2/{\left(p^{\prime }\right)}^2=3p^2/{\left(p^{\prime }\right)}^2,$$

the solution to which can be written as $S(p,p^{\prime },p^{″})=0$ where S is the corresponding score function (8).

To obtain the corresponding prior on $(-\infty ,\infty )$ through symmetry about 0, we obtain

$$q\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}a{(x-a)}^{-2},\hfill & \mathrm{for}x<0\hfill \\ \frac{1}{2}a{(x+a)}^{-2},\hfill & \mathrm{for}x\ge 0.\hfill \end{array}\right.$$

(9)

Here we motivate the natural objective choice for the constant a as 1. We note that a is a scale parameter in the above and as such the most fundamental transformation to be considered here is $\varphi =1/x$. This, for example, would take variance to precision (and vice versa). For the prior in (7) to be invariant, i.e., $p_a\left(\varphi \right)=a/{(a+\varphi )}^2$, we need to have $a=1$, since

$$\begin{array}{c}\hfill p_a\left(\varphi \right)=\frac{a}{{(a+1/\varphi )}^2}\left|\frac{1}{{\varphi }^2}\right|=\frac{a}{{(a\varphi +1)}^2},\end{array}$$

which yields $p_a\left(x\right)=a/{(a+x)}^2$ iff $a=1$. All the illustrations that follow have been made taking this choice for a. Although other transformations are clearly available, the dominance of the reciprocal transform is sufficient motivation to select the $a=1$.

First Examples

The first simulation study was to make inference on a scale parameter; specifically, the standard deviation of a normal density with mean $\mu =0$ and standard deviation $\sigma \in (0,\infty )$. We compare prior (7) with Jeffreys prior, i.e., $\pi \left(\sigma \right)\propto 1/\sigma$. We took 250 samples of size $n=30$ and $n=100$, obtained the posterior distributions using standard MCMC methods (Metropolis-Hastings with normal random walk proposal, 6000 iterations, with a burn–in of 1000 and a thinning of 10) and computed the following two indexes. The root Mean Squared Error (MSE) divided by the true parameter value from the sample mean;

$$MSE=\frac{\sqrt{\mathbb{E}{(\widehat{\sigma }-\sigma )}^2}}{\sigma },$$

where $\widehat{\sigma }$ is the posterior mean, and the coverage of the 95% posterior credible interval for $\sigma$.

Table 1. Posterior MSE for $\sigma$ to compare Jeffreys prior against the prior based on scores for data generated by a normal density with mean $\mu =0$ and unknown variance. This has been obtained on 250 samples of size $n=30$ and $n=100$ and with standard deviation $\sigma =\{0.25,0.50,1,2,5,10,20\}$.

	$\mathit{n}=30$		$\mathit{n}=100$
$\mathit{\sigma }$	Jeffreys Prior	Score Prior	Jeffreys PRIOR	Score Prior
0.25	0.1368	0.1371	0.0723	0.0720
0.5	0.1367	0.1333	0.0721	0.0721
1	0.1365	0.1329	0.0721	0.0716
2	0.1349	0.1326	0.0719	0.0716
5	0.1347	0.1333	0.0722	0.0720
10	0.1349	0.1331	0.0723	0.0716
20	0.1353	0.1337	0.0722	0.0718

shows the results for the MSE for $\sigma =\{0.25,0.50,1,2,5,10,20\}$, where we see little difference between the performance of the two priors. As one would expect, the MSE is lower for the largest sample size. However, the important point is that the score prior is proper, an important property.

The coverage of the 95% posterior credible interval is shown in

Table 2. Posterior coverage of the 95% credible interval for $\sigma$ to compare Jeffreys prior against the prior based on scores for data generated by a normal density with mean $\mu =0$ and unknown variance. This has been obtained on 250 samples of size $n=30$ and $n=100$ and with standard deviation $\sigma =\{0.25,0.50,1,2,5,10,20\}$.

	$\mathit{n}=30$		$\mathit{n}=100$
$\mathit{\sigma }$	Jeffreys Prior	Score Prior	Jeffreys PRIOR	Score Prior
0.25	0.96	0.95	0.91	0.91
0.5	0.95	0.94	0.92	0.91
1	0.95	0.94	0.91	0.91
2	0.95	0.95	0.93	0.91
5	0.95	0.94	0.93	0.91
10	0.94	0.94	0.90	0.92
20	0.95	0.95	0.90	0.93

, where we can also see very similar behavior between the two priors.

To illustrate the frequentist properties of the prior in (9), we have compared it to a flat prior, $\pi \left(\mu \right)\propto 1$, in making inference for a location parameter of a log–normal density with unknown $\mu$ and known scale parameter $\sigma =1$. Similar to the previous case, we have drawn 250 samples of size $n=30$ and $n=100$ and computed the MSE and coverage of the 95% posterior credible interval. The values of $\mu$ considered were from the set $\{0,1,5,50,100\}$.

Table 3. Posterior MSE for $\mu$ to compare Jeffreys prior against the prior based on scores for data generated from a log-normal density with unknown location parameter $\mu$ and known scale parameter $\sigma =1$. This has been obtained on 250 samples of size $n=30$ and $n=100$ and with standard deviation $\mu =\{0,1,1,5,50,100\}$.

	$\mathit{n}=30$		$\mathit{n}=100$
$\mathit{\sigma }$	Jeffreys Prior	Score Prior	Jeffreys PRIOR	Score Prior
0	0.0290	0.0208	0.0114	0.0007
1	0.0325	0.0481	0.0086	0.0091
5	0.0325	0.0324	0.0086	0.0086
10	0.0325	0.0323	0.0086	0.0086
50	0.0325	0.0323	0.0086	0.0087
100	0.0325	0.0324	0.0086	0.0087

shows the MSE for the two priors, where we see that apart for a small difference for $\mu =0$, the two priors appear to perform in a very similar fashion.

The coverage of the 95% posterior credible interval for $\mu$ is shown in

Table 4. Posterior coverage of the 95% credible interval for $\mu$ to compare Jeffreys prior against the prior based on scores for data generated from a log-normal density with unknown location parameter $\mu$ and known scale parameter $\sigma =1$. This has been obtained on 250 samples of size $n=30$ and $n=100$ and with standard deviation $\mu =\{0,1,1,5,50,100\}$.

	$\mathit{n}=30$		$\mathit{n}=100$
$\mathit{\sigma }$	Jeffreys Prior	Score Prior	Jeffreys PRIOR	Score Prior
0	0.96	0.98	0.92	1.00
1	0.94	0.92	0.97	0.98
5	0.94	0.94	0.97	0.98
10	0.94	0.94	0.97	0.98
50	0.94	0.94	0.97	0.98
100	0.94	0.94	0.97	0.98

, where we can see a very similar behavior for the two priors, with an exception for $\mu =0$, although the two coverage levels are perfectly acceptable.

The general conclusion for the two experiments above, is that the prior obtained via $\alpha \left(u\right)=1/u^2$ exhibits tails which are sufficiently heavy to generate optimal frequentist performance even for large parameter values. These properties are comparable to those obtained by Jeffreys prior, which is well- known for being the objective prior yielding good frequentist properties of the posterior. The advantage with our prior is that it is always proper.

3. Mixture Models

A major area of challenge for objective priors is finite mixture models, where observations are assumed to be generated by the following model;

$$f\left(y\right)=\sum _{l=1}^k{\omega }_l{f}_l(y\mid {\theta }_l),\sum _{l=1}^k{\omega }_l=1,$$

(10)

for densities $\left(f_l(·|{\theta }_l)\right)$, where ${\theta }_l$ is vector of parameters characterizing the densities. Given the ill-defined nature of the model in (10), Ref. [6], the use of improper priors for the parameters is not acceptable. In particular, if we consider densities $f_l$ to be location-scale distributions, Ref. [6] shows that under certain circumstances, Jeffreys priors cannot be used, due to their improperness. For example, if all parameters are unknown (i.e., weights, location and scale parameters), then Jeffreys prior yields improper posteriors. Even in more restrictive circumstances the use of improper priors if troublesome; if we consider only the location parameters unknown, then Jeffreys prior yields improper posteriors for $k>2$. The above represents severe limitations in Bayesian analysis. Therefore, the objective priors proposed in this paper represent a clear solution to the above problem, avoiding reparameterization or addition of extra parameters, as proposed, for example in [6], the latter resulting in an increased uncertainty.

3.1. Single Sample

In this first simulation study, we illustrate the performance of the proposed prior on the following three-component mixture model, from which we have drawn a sample of size $n=200$,

$$0.25\mathrm{N}(0,{1.2}^2)+0.65\mathrm{N}(-10,1)+0.10\mathrm{N}(7,{0.8}^2).$$

(11)

In terms of prior distributions, we have assumed a symmetric Dirichlet prior with concentration parameters equal to one, and for the means and standard deviations of the Gaussian components the proposed prior on $(-\infty ,\infty )$ and on $(0,\infty )$, respectively. We have also assumed independent information, so the priors for the parameters of the component have the following form,

$$\pi (\mathit{\mu },\mathit{\sigma })=\prod _{l=1}^3{\pi }_l\left({\mu }_l\right)\times \prod _{l=1}^3\pi \left({\sigma }_l\right),$$

where $\mathit{\mu }=({\mu }_1,{\mu }_2,{\mu }_3)$ and $\mathit{\sigma }=({\sigma }_1,{\sigma }_2,{\sigma }_3).$ The histogram of the sample, together with the true density, is shown in Figure 2.

The analysis uses a Metropolis–within–Gibbs algorithm with normal random walk proposal, with a total of 60,000 iterations, a burn–in of 10,000 and a thinning of 100. In

Table 5. Posterior means and 95% credible intervals (in brackets) for a sample of size $n=200$ from model (11).

Component	$\mathit{\omega }$	$\mathit{\mu }$	$\mathit{\sigma }$
1	0.26	−10.2	1.1
	(0.21, 0.33)	(−10.5, −9.9)	(0.9, 1.4)
2	0.66	0.0	1.3
	(0.58,0.72)	(−0.1, 0.2)	(1.1, 1.5)
3	0.08	6.7	0.9
	(0.04, 0.12)	(6.2, 7.2)	(0.6, 1.4)

we have reported the posterior means and the 95% posterior credible intervals for the parameters of the mixture model, where we note that the true values are within the posterior credible intervals.

3.2. Repeated Sampling

To have a more thorough understanding of the proposed prior implementation, we have performed some experiment on repeated sampling, taking into consideration different scenarios, which include different sample sizes and model structure. We have limited the analysis to mixture of normal densities, but it is obvious that due to the properness of the prior, its implementation can be extended to any family of densities for the mixture components. We computed the posterior distribution for $M=20$ replications of sample of size $n=(50,100,200)$ of mixture models with number of components $k=(3,4,5)$. For these illustrations, we reported the results for the means and the standard deviations of each component, as they are estimated using the proposed prior. The models used are as follows:

$$\begin{array}{ccc}\hfill \frac{1}{3}\mathrm{N}(-10,1)& +& \frac{1}{3}\mathrm{N}(0,{0.8}^2)+\frac{1}{3}\mathrm{N}(7,{1.2}^2)\hfill \\ \hfill \frac{1}{4}\mathrm{N}(-10,1)& +& \frac{1}{4}\mathrm{N}(-3,0.9)+\frac{1}{4}\mathrm{N}(0,0.8)+\frac{1}{4}\mathrm{N}(7,1.2)\hfill \\ \hfill \frac{1}{5}\mathrm{N}(-10,1)& +& \frac{1}{5}\mathrm{N}(-3,0.9)+\frac{1}{5}\mathrm{N}(0,0.8)+\frac{1}{5}\mathrm{N}(3,1)+\frac{1}{5}\mathrm{N}(7,1.2).\hfill \end{array}$$

Please note that we have not chosen variable weights as these are not associated with a proposed prior.

Figure 3 shows the boxplots of the posterior means for the $\mathit{\mu }$ of the mixture components, while Figure 4 shows the boxplots of the posterior means for the standard deviations $\mathit{\sigma }$. As one would expect, the larger the sample size the less variability in the repeated estimates, for the same number of components. Keeping the sample size fixed, we notice an increase in the variability of the estimates as the number of components grows, which is also an expected result.

3.3. Real Data Analysis

In this section, we analyze the well-known galaxy data set, which contains the velocity of 82 galaxies in the Corona Borealis region. To support a particular theory about the formation of galaxies, the analysis aims to estimate the number of stellar populations. This is a benchmark data set, well investigated in the literature, for example in [14,15,16], among others. We consider the galaxies velocities as random variables distributed according to a mixture of k normal densities. The estimation of the number of components has proved to be delicate, with estimates ranging from 3 to 7, depending on factors, such as the priors for the parameters and the Bayesian method used. The histogram of the data set with a superimposed density is presented in Figure 5.

To select the number of components in the mixture of normal densities, we have fitted models with $k=(2,3,4,5,6,7,8)$ components, computing the Deviance Information Criterion (DIC) under each model. The results are reported in

Table 6. Deviance Information Criterion for mixture with number of components $k=(2,3,4,5,6,7,8)$.

k	2	3	4	5	6	7	8
DIC	476.72	425.20	371.48	413.68	446.45	446.54	458.81

. We notice that according to the computed index, we identify as best model the mixture with $k=4$ components, which is in line with [16], and slightly more conservative than [15,16], where the number of components with non-zero weight is 5.

Table 7. Posterior means of the parameters for the $k=5$ components. The results are reported by descending order of the weights.

Component	$\mathit{\omega }$	$\mathit{\mu }$	$\mathit{\sigma }$
1	0.32	19.16	0.71
2	0.31	23.59	3.46
3	0.29	22.01	4.59
4	0.10	9.31	0.58

shows the posterior means for the parameters of the 4 components estimated.

4. Variance Parameters in Hierarchical Models

In this section, we discuss a well-known implementation of a hierarchical model that is proposed, for example, in [17]. The basic two-level hierarchical model is as follows:

$$\begin{array}{ccc}\hfill y_{ij}& \sim & N(\mu +{\alpha }_j,{\sigma }_y^2),i=1,\dots ,n,j=1,\dots ,J\hfill \\ \hfill {\alpha }_j& \sim & N(0,{\sigma }_{\alpha }^2),j=1,\dots ,J.\hfill \end{array}$$

This model has three parameters, namely $\mu$, ${\sigma }_y$ and ${\sigma }_{\alpha }$. However, out interest for this paper is in ${\sigma }_{\alpha }$ only, noting that “regular” objective priors can be used on the remaining parameters, such as $\pi (\mu ,{\sigma }_y)\propto 1$, for example. Although being improper, this prior yields a proper posterior on the parameters.

The actual concern is on the variance (scale) parameter ${\sigma }_{\alpha }$, as if we were to put an improper prior on it, then the corresponding posterior, most likely, would be improper as well. To compare the proposed prior, we assign an inverse-gamma prior on the variance with parameters set so to define a very sparse probability distribution. This is recommended, for example, in [18], where the prior is $\pi \left({\sigma }_{\alpha }^2\right)\sim IG(\epsilon ,\epsilon )$, with $\epsilon >0$ sufficiently small. We do not discuss in detail the appropriateness of the above choice, or other alternatives; the reader can refer to [17], for example, for a through discussion.

The method used to obtain the posterior samples is a Metropolis-within-Gibbs in both cases, with 40,000 iterations, a burn-in of 20,000 and a thinning of 10.

The data consists of $J=8$ educational testing experiments, where the parameters ${\alpha }_1,\dots ,{\alpha }_8$ represent the relative effects of Scholastic Aptitude Test coaching programs in different schools. In this example, the parameter ${\sigma }_{\alpha }$ represents the between-schools variability (standard deviation) of the effects.

Table 8. Observed effects $\left(y_j\right)$ of special preparation on SAT scores on eight randomized experiments, where ${\sigma }_j$ are the standard errors of effect estimate.

School	${\mathit{y}}_{\mathit{j}}$	${\mathit{\sigma }}_{\mathit{j}}$
A	28	15
B	8	10
C	−3	16
D	7	11
E	−1	9
F	1	11
G	18	10
H	12	18

shows the data.

We have compared the marginal posteriors $\pi \left({\sigma }_{\alpha }^2\right|\mathbf{y})$ obtained using the inverse-gamma prior with $\epsilon =1$ and the proposed prior in (7) with $a=1$. The histograms of the marginal posteriors are in Figure 6, where we note similar results. The statistics of the posteriors are reported in

Table 9. Posterior statistics for the marginal distribution of ${\sigma }_{\alpha }^2$ for the schools problem.

Prior	Mean	95% C.I.
Inverse-Gamma	3.8	(0.3, 23.9)
Score-based	2.3	(0.2, 57.8)

, where we note a less-informative distribution when the proposed prior is employed. This is expected, as the inverse-gamma distribution is considered a relatively informative one [17].

5. Discussion

In this paper, we have derived a class of objective prior distributions that have the appealing properties of being proper and heavy-tailed. These have been obtained by exploiting a straightforward approach to the construction of score functions (here proposed). In detail, using convex function $\alpha (·)$ we can find the score function with first two derivatives using (5). The Hyvärinen score arises with $\alpha \left(u\right)=u^2$; whereas we have used $\alpha \left(u\right)=u^{-2}$ and used it to construct objective prior distributions using methodology introduced in [7].

The class of prior is heavy-tailed, behaving as $1/x^2$ for large $\left|x\right|$; this result is immediately obvious as the prior on $(0,\infty )$ is a Lomax distribution with shape 1. In this respect, it behaves similar to standard objective priors but comes without the problems of being improper. The benefits of using a proper prior is that the posterior is automatically proper and so does not need to be checked.

We have showed that when compared to Jeffreys prior on simulated data, the frequentist performance of the prior distribution derived from score functions are nearly equivalent. In addition, we have showed that on both simulated and real data, the proposed prior is suitable to be used in a key scenario where improper priors (e.g., Jeffreys and reference) are not suitable (or are yet to be found). We have also illustrated the prior on a common problem for hierarchical models, i.e., assigning an objective prior for the variance parameter.

As a final point, we briefly discuss the case where a prior is needed on a multidimensional parameter space. Therefore, say we have a model with k parameters, i.e., $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_k)$, where ${\theta }_j\in {\Theta }_j$, for $j=1,\dots ,k$. We also assume that the uni-dimensional space for each parameter is either $(0,\infty )$ or $(-\infty ,\infty )$. Assuming k relatively large, besides some specific statistical models such as regression models or graphical models, a common practice to assign objective priors on $\mathit{\theta }$ is as follows:

$$\pi \left(\mathit{\theta }\right)={\pi }_1\left({\theta }_1\right)\times \dots \times {\pi }_k\left({\theta }_k\right).$$

In other words, parameters are assumed to be independent a priori, so the join prior distribution is represented by the product of the marginal priors on each parameter. We can then set ${\pi }_j\left({\theta }_j\right)$ to be either (7) or (9), for $j=1,\dots ,k$, depending on ${\Theta }_j$.