text:
# Combining cosmological datasets: hyperparameters and Bayesian evidence

## 1 Introduction

It is now common practice in cosmology to estimate the values of cosmological parameters by a joint analysis of a number of different datasets. The standard technique for performing such an analysis is to assume that the datasets are statistically independent and so take the joint likelihood function for the parameters to be given simply by the product of the individual likelihood functions for each separate dataset. The joint likelihood function can then be used in the standard way to determine the optimal values of the parameters and their associated errors.

As discussed by Lahav et al. (2000; hereinafter Paper I), however, there exists some freedom in the relative ‘weight’ that may be given to each dataset in the analysis (see also Godwin & Lynden-Bell 1987; Press 1996). The assignment of weights often occurs when two or more of the datasets are inconsistent, and is usually made in a somewhat ad-hoc manner. Typically some datasets are excluded from the analysis, and hence given a weight of zero, while the remainder are analysed jointly with equal weights. Despite its widespread use, this procedure has many unsatisfactory features, not least of which is the subjectivity associated with the choice of which datasets to include, and which to discard. As advocated in Paper I, a more objective procedure for assigning weights to the datasets is provided by the introduction of hyperparameters. This device allows the statistical properties of the data themselves to determine the weight given to each dataset in the analysis.

In Paper I, a method was presented for introducing hyperparameters into the analysis of datasets for which each likelihood function had a particular simple form. This technique was then applied to the estimation of the Hubble parameter $`h`$ from several sets of observations of the power spectrum of cosmic microwave background (CMB) anisotropies. It was shown that apparently discrepant datasets could be analysed jointly to provide a consistent estimate of $`h`$, together with an associated error. This approach has also recently been applied to the joint analysis of the baryon mass fraction in clusters and cepheid-calibrated distances (Erdogdu, Ettori & Lahav 2002).

In this paper, we extend the work of Paper I to accommodate more general situations. In section 2, we review the standard approach to Bayesian parameter estimation and discuss the role of the Bayesian evidence in model selection. In section 3, we present a general account of hyperparameters and their use in joint estimation of parameters. We also discuss how to use the Bayesian evidence to decide whether or not the data support the inclusion of hyperparameters in the first instance. In section 4, we consider the use of hyperparameters in the weighting of datasets and, in section 5, we discuss the common special case in which the likelihood function for each dataset is Gaussian. In section 6, we illustrate the useful properties of the hyperparameters approach by applying the technique to the toy problem of fitting a straight-line to a number of datasets. A cosmological application of the hyperparameters approach is discussed in section 7, in which we perform a joint estimate of the physical baryon density $`\mathrm{\Omega }_bh^2`$ and the scalar spectral index $`n`$ using the most recent sets of observations of the CMB power spectrum. Finally, our conclusions are presented in section 8.

## 2 Bayes’ theorem and evidence

Suppose the totality of our data is represented by the data vector D and we are interested in estimating the values of the parameters $`𝜽`$ in some underlying model of the data. The standard approach to this problem is to use Bayes’ theorem

$$\mathrm{Pr}(𝜽|\text{D})=\frac{\mathrm{Pr}(\text{D}|𝜽)\mathrm{Pr}(𝜽)}{\mathrm{Pr}(\text{D})},$$
(1)
which gives the posterior distribution $`\mathrm{Pr}(𝜽|\text{D})`$ in terms of the likelihood $`\mathrm{Pr}(\text{D}|𝜽)`$, the prior $`\mathrm{Pr}(𝜽)`$ and the evidence $`\mathrm{Pr}(\text{D})`$ (which is also often called the marginalised likelihood).

### 2.1 Parameter estimation

For the purpose of estimating parameters, one usually ignores the normalisation factor $`\mathrm{Pr}(\text{D})`$ in Bayes’ theorem, since it does not depend on the parameters $`𝜽`$. Thus, one instead works with the ‘unnormalised posterior’

$$\overline{\mathrm{Pr}}(𝜽|\text{D})\mathrm{Pr}(\text{D}|𝜽)\mathrm{Pr}(𝜽),$$
(2)
where we have written $`\overline{\mathrm{Pr}}`$ to denote the fact that the ‘probability distribution’ on the left-hand side is not normalised to unit volume. In fact, it is also common to omit normalising factors, that do not depend on the parameters $`𝜽`$, from the likelihood and the prior. As we shall see below, however, if one wishes to calculate the Bayesian evidence for a particular model, the likelihood and the prior must be properly normalised such that $`\mathrm{Pr}(𝜽)𝑑𝜽=1`$ and $`\mathrm{Pr}(\text{D}|𝜽)𝑑\text{D}=1`$. We will therefore assume here that the necessary normalising factors have been retained.

Strictly speaking, the entire (unnormalised) posterior is the Bayesian inference of the parameters values. Unfortunately, if the dimension $`M`$ of the parameter space is large, it is often numerically unfeasible to calculate $`\overline{\mathrm{Pr}}(𝜽|\text{D})`$ on some $`M`$-dimensional hypercube. Thus, particularly in large problems, it is usual to present one’s results in terms of the ‘best’ estimates $`\widehat{𝜽}`$, which maximise the (unnormalised) posterior, together with some associated errors. These errors are usually quoted in terms of the estimated covariance matrix

$$\text{C}=\left[\mathrm{ln}\overline{\mathrm{Pr}}(𝜽|\text{D})|_{𝜽=\widehat{𝜽}}\right]^1.$$
(3)
or as confidence limits on each parameter $`a_i`$ $`(i=1,\mathrm{},M)`$, obtained from the one-dimensional marginalised (unnormalised) posterior distributions

$$\overline{\mathrm{Pr}}(a_i|\text{D})=\overline{\mathrm{Pr}}(𝜽|\text{D})𝑑\stackrel{ˇ}{𝜽},$$
(4)
where $`d\stackrel{ˇ}{𝜽}=da_1\mathrm{}d\stackrel{ˇ}{a}_i\mathrm{}da_M`$ denotes that the integration is performed over all other parameters $`a_j`$ $`(ji)`$.

The estimates $`\widehat{𝜽}`$ are most often obtained by an iterative numerical minimisation algorithm. Indeed, standard numerical algorithms are generally able to locate a local (and sometimes global) maximum of this function even in a space of large dimensionality. Similarly, the covariance matrix of the errors can be found straightforwardly by first numerically evaluating the Hessian matrix $`\overline{\mathrm{Pr}}(𝜽|\text{D})`$ at the peak $`𝜽=\widehat{𝜽}`$, and then calculating (minus) its inverse.

### 2.2 Bayesian evidence and model selection

The standard technique outlined above produces inferences of the parameter values for a given model of the data, but it does not provide a mechanism for deciding which one of a set of alternative models is most suitable for describing the data. This problem may be addressed, however, using the Bayesian evidence $`\mathrm{Pr}(\text{D})`$. Very readable introductions to this topic are given by Bishop (1995) and Sivia (1996).

Although the evidence term is usually ignored in the process of parameter estimation, it is central to selecting between different models for the data. For illustration, let us suppose we have two alternative models (or hypotheses) for the data D; these hypotheses are traditionally denoted by $`H_0`$ and $`H_1`$. Let us assume further that the model $`H_0`$ is characterised by the parameter set $`𝜽`$, whereas $`H_1`$ is described by the set of parameters $`\mathit{\varphi }`$. For the model $`H_0`$, the probability density for an observed data vector D is given by

$$\mathrm{Pr}(\text{D}|H_0)=\mathrm{Pr}(\text{D}|𝜽)\mathrm{Pr}(𝜽)𝑑𝜽,$$
(5)
where, on the left-hand side, we have made explicit the conditioning on $`H_0`$. Similarly, for the model $`H_1`$,

$$\mathrm{Pr}(\text{D}|H_1)=\mathrm{Pr}(\text{D}|\mathit{\varphi })\mathrm{Pr}(\mathit{\varphi })𝑑\mathit{\varphi }.$$
(6)
In either case, we see that the evidence is given by the average of the likelihood function with respect to the prior. Thus, a model will have a larger evidence if more of its allowed parameter space is likely, given the data. Conversely, a model will have a small evidence if there exist large areas of the allowed parameter space with low values of the likelihood, even if the likelihood function is strongly peaked and the corresponding model predictions agree closely with the data. Hence the value of the evidence naturally incorporates the spirit of Ockham’s razor: a simpler theory, having a more compact parameter space, will generally have a larger evidence than a more complicated theory, unless the latter is significantly better at explaining the data.

The question of which is the models $`H_0`$ and $`H_1`$ is prefered is thus answered simply by comparing the relative values of the evidences $`\mathrm{Pr}(\text{D}|H_0)`$ and $`\mathrm{Pr}(\text{D}|H_1)`$. The hypothesis having the larger evidence is the one that should be accepted. Although not widely used in cosmology, the idea of model selection using evidence ratios has been considered previously in this context by Jaffe (1996) and, more recently, by John & Narlikar (2001).

### 2.3 The Gaussian approximation to the posterior

Unfortunately, the evaluation of an evidence integral, such as (5), is a challenging numerical task. From (2), we see that

$$\mathrm{Pr}(\text{D}|H_0)=\overline{\mathrm{Pr}}(𝜽|\text{D})𝑑𝜽,$$
(7)
and so the evidence may only be evaluated directly if one can calculate $`\overline{\mathrm{Pr}}(𝜽|\text{D})`$ over some hypercube in parameter space, which we noted earlier is often computational unfeasible. Nevertheless, if the data are conclusive, we would expect this (unnormalised) posterior to be sharply peaked about the position of the maximum $`𝜽=\widehat{𝜽}`$. In this case, we may approximate this function by performing a Taylor expansion about $`𝜽=\widehat{𝜽}`$. Working with the log-posterior, and keeping only terms up to second-order in $`𝜽`$, leads to the Gaussian approximation of the (unnormalised) posterior distribution, which reads

$$\overline{\mathrm{Pr}}(𝜽|\text{D})\overline{\mathrm{Pr}}(\widehat{𝜽}|\text{D})\mathrm{exp}\left[\frac{1}{2}(𝜽\widehat{𝜽})^\mathrm{T}\text{C}^1(𝜽\widehat{𝜽})\right],$$
(8)
where $`\text{C}^1`$ is the estimated inverse covariance matrix of the parameters and is given by

$$\text{C}^1=\mathrm{ln}\mathrm{Pr}(𝜽|\text{D})|_{𝜽=\widehat{𝜽}}=\mathrm{ln}\overline{\mathrm{Pr}}(𝜽|\text{D})|_{𝜽=\widehat{𝜽}}.$$
(9)
Substituting the form (8) into (7), we thus find that an approximation to the value of the evidence is given by

$$\mathrm{Pr}(\text{D}|H_0)(2\pi )^{M/2}|\text{C}|^{1/2}\mathrm{Pr}(\widehat{𝜽})\mathrm{Pr}(\text{D}|\widehat{𝜽})$$
(10)
where $`M`$ is the number of parameters of interest $`𝜽`$ and we have rewritten $`\overline{\mathrm{Pr}}(\widehat{𝜽}|\text{D})`$ using (2). Since the estimation of parameter values and their associated errors already requires one to calculate all the quantities on the right-hand side of (10), we see that this approximate evidence may be evaluated with no extra work. We note, however, that for (10) to hold, the prior and likelihood must be correctly normalised, such that $`\mathrm{Pr}(𝜽)𝑑𝜽=1`$ and $`\mathrm{Pr}(\text{D}|𝜽)𝑑\text{D}=1`$. We also note, in this case, that consideration of the ratio of evidences is a natural generalisation of the standard likelihood-ratio approach to model comparison.

### 2.4 Markov-Chain Monte-Carlo methods

We note, in passing, that the approach to parameter estimation and the approximate evaluation of evidences outlined above may soon become obsolete. With the advent of faster computers and efficient algorithms, it has recently become numerically feasible to sample directly from the posterior distribution using Monte-Carlo Markov-Chain (MCMC) techniques (see Knox, Christensen & Skordis 2001). This allows one trivially to obtain one-dimensional marginalised posteriors for each parameter of interest $`a_i`$, and also enables the direct numerical evaluation of evidence integrals. Clearly, the MCMC technique has enormous potential for the future estimation of cosmological parameters.

## 3 Hyperparameters

In this paper, we wish to construct a robust technique for performing a joint estimation of cosmological parameters from combined datasets. The basic idea behind the approach presented here is to introduce additional hyperparameters $`𝜶`$ into the Bayesian inference problem. In other words, we extend our parameter vector to include not only the parameters of interest $`𝜽`$, but also the hyperparameters $`𝜶`$. These hyperparameters are analogous to ‘nuisance’ parameters that often arise in the standard approach to parameter estimation. In our case, however, they are not present in our model a priori, but we have chosen to introduce them in order to allow extra freedom in the parameter estimation process.

### 3.1 Marginalisation over hyperparameters

As with any set of nuisance parameters, we must integrate out (or marginalise over) the hyperparameters $`𝜶`$ in order to recover the posterior distribution of our parameters of interest $`𝜽`$. Thus, we have

$`\mathrm{Pr}(𝜽|\text{D})`$ $`=`$ $`{\displaystyle \mathrm{Pr}(𝜽,𝜶|\text{D})𝑑𝜶}`$ (11)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{Pr}(\text{D})}}{\displaystyle \mathrm{Pr}(\text{D}|𝜽,𝜶)\mathrm{Pr}(𝜽,𝜶)𝑑𝜶}.`$
In this paper, we shall assume that the parameters of interest $`𝜽`$ and the hyperparameters $`𝜶`$ are independent, so that

$$\mathrm{Pr}(𝜽,𝜶)=\mathrm{Pr}(𝜽)\mathrm{Pr}(𝜶).$$
(12)
Substituting this form of the prior into (11), we immediately recover Bayes’ theorem (1), where the form of the likelihood function in the presence of hyperparameters is

$$\mathrm{Pr}(\text{D}|𝜽)=\mathrm{Pr}(\text{D}|𝜽,𝜶)\mathrm{Pr}(𝜶)𝑑𝜶.$$
(13)
Indeed, under the assumption (12), this expression for the likelihood embodies the complete hyperparameters technique.

Let us now turn our attention to the form of this likelihood function. Assuming that the totality of data D consists of $`N`$ independent datasets $`\text{D}_k`$ $`(k=1,\mathrm{},N)`$, we have

$$\mathrm{Pr}(\text{D}|𝜽,𝜶)=\underset{k=1}{\overset{N}{}}\mathrm{Pr}(\text{D}_k|𝜽,𝜶).$$
(14)
Moreover, in this paper, we will make the additional simplifying assumption that the $`k`$th term in the product on the right-hand side of (14) depends only on $`\alpha _k`$. In this case, it reduces to

$$\mathrm{Pr}(\text{D}|𝜽,𝜶)=\underset{k=1}{\overset{N}{}}\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k).$$
(15)
We shall also assume that the individual hyperparameters $`\alpha _i`$ $`(i=1,\mathrm{},N)`$ are themselves independent, so that $`\mathrm{Pr}(𝜶)=\mathrm{Pr}(\alpha _1)\mathrm{Pr}(\alpha _2)\mathrm{}\mathrm{Pr}(\alpha _N)`$. Thus the expression (13) for the likelihood becomes

$$\mathrm{Pr}(\text{D}|𝜽)=\underset{k=1}{\overset{N}{}}\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)\mathrm{Pr}(\alpha _k)𝑑\alpha _k,$$
(16)
which is simply the product of the individual likelihoods $`P(\text{D}_k|𝜽)`$ for each dataset, after marginalisation over the hyperparameter. This form can be substituted into Bayes’ theorem (1) to obtain the posterior $`\mathrm{Pr}(𝜽|\text{D})`$, which can then be used to obtain the best estimates $`\widehat{𝜽}`$ of the parameters and their associated errors.

Finally, it is worth noting that we may use (16) to write the standard approach to parameter estimation, in which no hyperparameters are introduced, as a special case of the hyperparameters technique. This is achieved by fixing the hyperparameters $`\alpha _k`$ to have particular values $`\alpha _k^0`$ $`(k=1,\mathrm{},N)`$, which is easily accommodated in the above formalism by assigning the priors

$$\mathrm{Pr}(\alpha _k)=\delta (\alpha _k\alpha _k^0),$$
(17)
where $`\delta (z)`$ is the Dirac delta function. Most often, the hyperparameters are introduced in such a way that $`\alpha _k^0=1`$ (for all $`k`$) corresponds to the standard approach, in which hyperparameters are absent. In other words, the standard likelihood function for the $`k`$th dataset may be denoted by $`\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)`$.

### 3.2 Hyperparameters and evidence

So far, we have not addressed the question of whether the data support the introduction of hyperparameters in the first instance. For example, if a collection of different datasets are all mutually-consistent, then one might not consider it wise to introduce the hyperparameters $`𝜶`$. Indeed, their introduction could lead to larger uncertainties on the estimated values of the parameters of interest $`𝜽`$, since one has to perform a marginalisation over $`𝜶`$. On the other hand, if several of the datasets are not in good agreement, the use of hyperparameters is necessary, in order to obtain statistically meaningful results.

In fact, Bayes’ theorem itself allows us to make an objective decision regarding whether the data warrant the introduction of hyperparameters. As discussed in section 2.2, this may be achieved using the Bayesian evidence $`\mathrm{Pr}(\text{D})`$. For illustration, let us consider two models for the data as follows: the model $`H_0`$ does not include hyperparameters, whereas the model $`H_1`$ assigns a free hyperparameter to each dataset. In either case, the probability of obtaining the observed data D is

$`\mathrm{Pr}(\text{D})`$ $`=`$ $`{\displaystyle \mathrm{Pr}(\text{D},𝜽,𝜶)𝑑𝜶𝑑𝜽}`$ (18)
$`=`$ $`{\displaystyle \mathrm{Pr}(\text{D}|𝜽,𝜶)\mathrm{Pr}(𝜽,𝜶)𝑑𝜶𝑑𝜽}`$
$`=`$ $`{\displaystyle \mathrm{Pr}(\text{D}|𝜽)\mathrm{Pr}(𝜽)𝑑𝜽},`$
where, in the last line, the likelihood $`\mathrm{Pr}(\text{D}|𝜽)`$ is given by (16). For $`H_0`$, the priors are simply $`\mathrm{Pr}(\alpha _k)=\delta (\alpha _k1)`$ and the integral simplifies accordingly, whereas for $`H_1`$ the priors will, in general, have some more complicated form.

Denoting the resulting evidence values for our two hypotheses by $`\mathrm{Pr}(\text{D}|H_0)`$ and $`\mathrm{Pr}(\text{D}|H_1)`$ respectively, the question of whether the data warrant the introduction of hyperparameters, is now answered simply by comparing the relative values of $`\mathrm{Pr}(\text{D}|H_0)`$ and $`\mathrm{Pr}(\text{D}|H_1)`$. The hypothesis having the larger evidence is the one that should be accepted. If $`\mathrm{Pr}(\text{D}|H_1)\mathrm{Pr}(\text{D}|H_0)`$, this indicates that the datasets are all mutually consistent and that the introduction of hyperparameters is not warranted. Conversely, the condition $`\mathrm{Pr}(\text{D}|H_1)\mathrm{Pr}(\text{D}|H_0)`$ is strong indication that the datasets are not mutually consistent, and that hyperparameters are necessary in order to obtain meaningful statistical results.

## 4 Weighting of datasets

In the above discussion, we have specialised to the case where a single hyperparameter is associated with each dataset. We have not yet, however, made explicit how this hyperparameter enters the form of the modified likelihood $`\mathrm{Pr}(\text{D}|𝜽,\alpha _k)`$ for each dataset. Neither have we fixed the form of the prior $`\mathrm{Pr}(\alpha _k)`$ on each hyperparameter. As mentioned in the Introduction, an obvious use of hyperparameters in cosmological parameter estimation is in the weighting of the different datasets being analysed. We now consider this application in more detail.

### 4.1 The prior on each hyperparameter

Let us first turn our attention to the prior $`\mathrm{Pr}(\alpha _k)`$ on each hyperparameter. As discussed in Paper I, since $`\alpha _k`$ is a scale parameter, we might adopt the ‘non-informative’ or Jeffrey’s prior, which is uniform in the logarithm of the parameter. Thus, we have

$$\overline{\mathrm{Pr}}(\alpha _k)=\frac{1}{\alpha _k},$$
(19)
(or, more generally, one might take $`\mathrm{Pr}(\alpha )=1/\alpha ^n`$, where $`n`$ is a non-negative integer). The Jeffrey’s prior does, however, cause some difficulties, since it is improper and so cannot be normalised (as indicated). This is not a problem for obtaining estimates of the parameter values $`\widehat{𝜽}`$ and their associated errors, since these quantities do not depend on the overall normalisation of the posterior. However, the use of an improper prior makes it impossible to calculate evidences, and so we are unable to assess the whether the data warrant the introduction of the hyperparameters.

We must ask ourselves, however, whether we are truly ignorant of the value of each hyperparameter $`\alpha _k`$. In fact, this is rather unlikely. Given the nature of the data-weighting problem, we might suppose that the expectation value of $`\alpha _k`$ is unity. In other words, in the first instance, we expect that the experimental data have been correctly analysed and that the results are free from any systematic bias or misquoted random errors. Thus, we need to assign a prior probability distribution $`\mathrm{Pr}(\alpha _k)`$ subject to the constraint that $`E[\alpha _k]=1`$. As shown by Jaynes (1957a,b), building on the work of Shannon (1948), the only consistent way of assigning the prior probability distribution $`\mathrm{Pr}(\alpha _k)`$ is by maximising the ‘entropy’ functional

$$S[\mathrm{Pr}(\alpha _k)]=_0^{\mathrm{}}\mathrm{Pr}(\alpha _k)\mathrm{ln}\mathrm{Pr}(\alpha _k)d\alpha _k,$$
(20)
subject to the normalisation constraint $`_0^{\mathrm{}}\mathrm{Pr}(\alpha _k)𝑑\alpha _k=1`$ and the constraint $`E[\alpha _k]=_0^{\mathrm{}}\alpha _k\mathrm{Pr}(\alpha _k)𝑑\alpha _k=1`$ on the expectation value. This is a straightforward problem in the calculus of variations, and has the simple solution

$$\mathrm{Pr}(\alpha _k)=\mathrm{exp}(\alpha _k),$$
(21)
which is, of course, properly normalised. In this paper, we shall use this exponential form of the prior.

We note in passing, however, that in some cases one may not only have $`E[\alpha _k]=1`$, but also some a priori expectation on the variance of $`\alpha _k`$, i.e. some limit on the range of weights that should be assigned to each dataset. In this case, we need to assign a prior probability distribution $`\mathrm{Pr}(\alpha _k)`$ by maximising (20) subject to the constraints $`_0^{\mathrm{}}\mathrm{Pr}(\alpha _k)𝑑\alpha _k=1`$, $`E[\alpha _k]=1`$ and $`V[\alpha _k]=\sigma ^2`$ (say). Unfortunately, this problem has no solution as stated. Nevertheless, a closed-form solution is easily obtained if one does not restrict $`\alpha _k`$ to be non-negative, but instead allows it to take any value between $`\mathrm{}`$ and $`\mathrm{}`$ (thereby modifying the limits on the integral in (20) and in the normalisation of $`\mathrm{Pr}(\alpha _k)`$). In this case, the calculus of variations problem is easily solved to obtain the Gaussian form

$$\mathrm{Pr}(\alpha _k)=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left[\frac{(\alpha _k1)^2}{2\sigma ^2}\right],$$
(22)
which is again properly normalised. Clearly, this form is not strictly applicable to weights, which are required to be non-negative. Nevertheless, provided $`\sigma 1`$, the above Gaussian form will be a good approximation to the true prior.

### 4.2 The likelihood function for each dataset

If the hyperparameters $`\alpha _k`$ are to act as weights on the different datasets, we specify the log-likelihood function for each dataset to have the form

$$\mathrm{ln}\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)=\alpha _k\mathrm{ln}\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)\mathrm{ln}Z_k(𝜽,\alpha _k),$$
(23)
where $`\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)`$ is the standard likelihood function (in the absence of hyperparameters) and $`Z_k`$ is a normalisation factor that ensures $`\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)𝑑\text{D}_k=1`$. The expression (23) for the log-likelihood clearly corresponds to the likelihood function itself having the form

$$\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)=\frac{\left[\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)\right]^{\alpha _k}}{Z_k(𝜽,\alpha _k)},$$
(24)
from which we see that the normalisation factor has the explicit form

$$Z_k(𝜽,\alpha _k)=\left[\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)\right]^{\alpha _k}𝑑\text{D}_k.$$
(25)
In particular, we note that $`Z_k(𝜽,1)=1`$. From (23), we see that the standard approach to joint parameter estimation corresponds to taking $`\alpha _k=1`$ for those datasets that are included in the analysis, and $`\alpha _k=0`$ for those that are discarded. Using the prior (21), and the form of the likelihood for each dataset given in (24), we find that the new likelihood function for the $`k`$ dataset is given by

$$\mathrm{Pr}(\text{D}_k|𝜽)=_0^{\mathrm{}}\frac{\left[\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)\right]^{\alpha _k}}{Z_k(𝜽,\alpha _k)}e^{\alpha _k}𝑑\alpha _k.$$
(26)
### 4.3 The full likelihood and posterior distribution

Since the datasets $`\text{D}_k`$ $`(k=1,\mathrm{},N)`$ are assumed to be independent, the full likelihood function is given by

$$\mathrm{Pr}(\text{D}|𝜽)=\underset{k=1}{\overset{N}{}}\mathrm{Pr}(\text{D}_k|𝜽),$$
(27)
where $`\mathrm{Pr}(\text{D}_k|𝜽)`$ is given by (26). The posterior distribution $`\mathrm{Pr}(𝜽|\text{D})`$ of the parameters of interest is then given by Bayes’ theorem (1).

For the purposes of illustration, let us suppose finally that the prior $`\mathrm{Pr}(𝜽)`$ on the parameters of interest is uniform. In order that we may subsequently calculate evidences, we must ensure that this prior is properly normalised. Suppose the prior is zero outside some (large) region $``$ is the $`M`$-dimensional parameter space. Thus, we may define the prior as

$$\mathrm{Pr}(𝜽)=\{\begin{array}{cc}1/V()& \text{if }𝜽\text{ lies in }\hfill \\ 0& \text{otherwise},\hfill \end{array}$$
(28)
where $`V()=_{}\mathrm{Pr}(𝜽)𝑑𝜽`$ is the ‘volume’ of the region $``$.

The resulting posterior distribution $`\mathrm{Pr}(𝜽|\text{D})`$ can then be calculated over some $`M`$-dimensional hypercube in parameter space, or used in the standard way to obtain the estimates $`\widehat{𝜽}`$ and their associated errors. We note that, since the hyperparameters $`\alpha _k`$ have been marginalised over in (27), they do not have specific values that can be quoted at the end of the analysis. Nevertheless, it can be useful to know which values of $`\alpha _k`$ were most favoured, and hence which datasets were given a large or small weight in the analysis. Thus, at each point in the space of parameters $`𝜽`$, we define the ‘effective’ weight $`\alpha _k^{\mathrm{eff}}(𝜽)`$ to be that which maximises the corresponding individual likelihood $`\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)`$ at that point. Of course, the most relevant set of such quantities are those evaluated at the point $`𝜽=\widehat{𝜽}`$.

## 5 Gaussian likelihood functions

Our analysis has thus far been presented in a general form, which may be applied to a wide range of problems in which one wishes to weight different datasets. It is quite common, however, for the likelihood function to take the form of a multivariate Gaussian in the data. In this case, for the $`k`$th dataset, we have

$$\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k=1)=\frac{1}{(2\pi )^{n_k/2}|\text{V}_k|^{1/2}}\mathrm{exp}\left(\frac{1}{2}\chi _k^2\right),$$
(29)
where

$$\chi _k^2=(\text{D}_k𝝁_k)^\mathrm{T}\text{V}_k^1(\text{D}_k𝝁_k).$$
(30)
In these expressions, $`n_k`$ is the number of items in the $`k`$th dataset, $`𝝁_k`$ is the expectation value of the quantities $`\text{D}_k`$ and $`\text{V}_k`$ is their covariance matrix. In general, both $`𝝁_k`$ and $`\text{V}_k`$ may depend on the parameters of interest $`𝜽`$. (note that, in general, the likelihood function is not a multivariate Gaussian in the parameter space $`𝜽`$, unless $`\text{V}_k`$ does not depend on the parameters $`𝜽`$ and $`𝝁_k`$ depends only linearly on them).

We note that, even when the form of the likelihood is not Gaussian, one can often perform a coordinate transformation to a new set of data variables $`\text{z}_k`$, for which the likelihood is (approximately) Gaussian (see, for example, Bond, Jaffe & Knox 2000). Moreover, as shown by Bridle et al. (2001) in the context of CMB observations, if the data $`\text{D}_k`$ originally obey a Gaussian likelihood function of the form (29), then one can perform analytic marginalisations over calibration and ‘beam’ uncertainties that yield a new likelihood function which is also of Gaussian form, but with a modified ‘covariance’ matrix $`\text{V}_k^{}`$.

In the standard approach to estimating the parameters $`𝜽`$, one assumes a model $`H_0`$, in which no hyperparameters are introduced. Thus the full likelihood function is given simply by

$$\mathrm{Pr}(\text{D}|𝜽,H_0)=\underset{k=1}{\overset{N}{}}\frac{1}{(2\pi )^{n_k/2}|\text{V}_k|^{1/2}}\mathrm{exp}\left(\frac{1}{2}\chi _k^2\right),$$
(31)
where, on the left-hand side, we have made explicit the conditioning on $`H_0`$. One then substitutes this expression into (2) to obtain the corresponding (unnormalised) posterior $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$.

Alternatively, we may adopt the model $`H_1`$, in which a free hyperparameter $`\alpha _k`$ is assigned as a weight to each dataset. Defining the modified likelihood function for each dataset $`\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)`$ by (24), and using (25) to evaluate the normalistion factor $`Z_k(𝜽,\alpha _k)`$, one finds

$$\mathrm{Pr}(\text{D}_k|𝜽,\alpha _k)=\frac{1}{(2\pi )^{n_k/2}|\text{V}_k|^{1/2}}\alpha _k^{n_k/2}\mathrm{exp}\left(\frac{1}{2}\alpha _k\chi _k^2\right),$$
(32)
This result may then be substituted into the expression (26) to obtain the new likelihood function for the $`k`$th dataset, which is given by

$$\mathrm{Pr}(\text{D}_k|𝜽)=\frac{1}{(2\pi )^{n_k/2}|\text{V}_k|^{1/2}}_0^{\mathrm{}}\alpha _k^{n_k/2}e^{\alpha _k({\scriptscriptstyle \frac{1}{2}}\chi _k^2+1)}𝑑\alpha _k.$$
(33)
The integral over $`\alpha _k`$ can be performed easily, using the definition of the Gamma function $`\mathrm{\Gamma }(n)=_0^{\mathrm{}}x^{n1}e^x𝑑x`$, and one finds

$$\mathrm{Pr}(\text{D}_k|𝜽)=\frac{2\mathrm{\Gamma }\left(\frac{n_k}{2}+1\right)}{\pi ^{n_k/2}|\text{V}_k|^{1/2}}(\chi _k^2+2)^{\left({\scriptscriptstyle \frac{n_k}{2}}+1\right)},$$
(34)
which is properly normalised. Thus, assuming the datasets to be independent, the full likelihood (27) is given by

$$\mathrm{Pr}(\text{D}|𝜽,H_1)=\underset{k=1}{\overset{N}{}}\frac{2\mathrm{\Gamma }\left(\frac{n_k}{2}+1\right)}{\pi ^{n_k/2}|\text{V}_k|^{1/2}}(\chi _k^2+2)^{\left({\scriptscriptstyle \frac{n_k}{2}}+1\right)},$$
(35)
where we have made explicit the conditioning on $`H_1`$. As above, this form of the likelihood can then be substituted into (2) to obtain the corresponding (unnormalised) posterior $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$.

### 5.1 Evaluation of the posterior and evidences

Ideally one would wish to calculate the full (unnormalised) posteriors $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ on some hypercube in parameter space. In this way, the location of the (global) maximum is obtained immediately, and the presence of multiple peaks in the posterior(s) is readily observed. Moreover, marginalised distributions may be trivially calculated, in order to place confidence limits on individual parameter values. Assuming the uniform prior (28), we see from (31) and (35) that the evaluation of the posterior for the models $`H_0`$ and $`H_1`$ respectively requires similar functions to be evaluated. In other words, if one has an algorithm for calculating $`|\text{V}_k|`$ and $`\chi _k^2`$, as required for the evaluation of the standard posterior (model $`H_0`$), then one can immediately evaluate the hyperparameters posterior (model $`H_1`$).

An additional advantageous feature of calculating the full (unnormalised) posteriors $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ on a hypercube in parameter space is that it allows the immediate evaluation of the evidence in each case, which are given by

$`\mathrm{Pr}(\text{D}|H_0)`$ $`=`$ $`{\displaystyle \overline{\mathrm{Pr}}(𝜽|\text{D},H_0)𝑑𝜽},`$ (36)
$`\mathrm{Pr}(\text{D}|H_1)`$ $`=`$ $`{\displaystyle \overline{\mathrm{Pr}}(𝜽|\text{D},H_1)𝑑𝜽}.`$ (37)
The ratio of these quantities can then be used to decide whether or not the inclusion of hyperparameters is warranted by the data.

### 5.2 Estimation of parameter values

Unfortunately, evaluation of the full posterior distribution(s) on a hypercube in parameter space is only feasible when the number of parameters $`M`$ is small (although, as mentioned in section 2.3, this restriction may be overcome using MCMC sampling techniques). For larger problems, the estimates $`\widehat{𝜽}`$ of the parameters of interest are usually obtained by maximising the posterior. Assuming the uniform prior (28), for the standard model $`H_0`$ (without hyperparameters) this corresponds to minimising the function

$$2\mathrm{ln}\mathrm{Pr}(𝜽|\text{D},H_0)=\underset{k=1}{\overset{N}{}}\left(\mathrm{ln}|\text{V}_k|+\chi _k^2\right)+c,$$
(38)
where $`c`$ is a constant. This is, of course, a very familiar result. For the model $`H_1`$, however, which contains hyperparameters, one must instead minimise

$$2\mathrm{ln}\mathrm{Pr}(𝜽|\text{D},H_1)=\underset{k=1}{\overset{N}{}}\left[\mathrm{ln}|\text{V}_k|+(n_k+2)\mathrm{ln}(\chi _k^2+2)\right]+c,$$
(39)
where $`c`$ is a (different) constant.

From (38), for the model $`H_0`$, the parameters estimates $`\widehat{𝜽}_0`$ may be obtained by solving

$$\underset{k=1}{\overset{N}{}}(\mathrm{ln}|\text{V}_k|+\chi _k^2)=0,$$
(40)
whereas for the model $`H_1`$, the estimates $`\widehat{𝜽}_1`$ are found by solving

$$\underset{k=1}{\overset{N}{}}[\mathrm{ln}|\text{V}_k|+(n_k+2)\mathrm{ln}(\chi _k^2+2)]=0.$$
(41)
In the last case, however, we may write

$$\mathrm{ln}(\chi _k^2+2)=(\chi _k^2)/(\chi _k^2+2).$$
(42)
Thus, if one is able to evaluate the derivatives $`\mathrm{ln}|\text{V}_k|`$ and $`\chi _k^2`$ for $`k=1,\mathrm{},N`$, to obtain the standard parameter estimates $`\widehat{𝜽}_0`$, one may easily obtain the estimates $`\widehat{𝜽}_1`$.

We note, in passing, that in the special case where the covariance matrices $`\text{V}_k`$ do not depend on the parameters $`𝜽`$, the expression (39) is very similar to that presented in Paper I. In that paper, an improper Jeffrey’s prior was assumed on the hyperparameters, and it was found that the parameter estimates $`\widehat{𝜽}_1`$ were given by minimising the quantity

$$\underset{k=1}{\overset{N}{}}n_k\mathrm{ln}(\chi _k^2).$$
(43)
Unfortunately, the corresponding posterior in this case cannot be normalised, as a result of using an improper prior. Nevertheless, in the limit where the number of data items $`n_k`$ in each dataset is large, we would expect $`\chi _k^2`$ also to be large. Thus, in this limit, minimisation of the expressions (39) and (43) respectively (in the case where the $`\text{V}_k`$ are independent of the parameters $`𝜽`$) should give almost identical parameter estimates $`\widehat{𝜽}_1`$.

As mentioned in section 4.3, since we have marginalised over the hyperparameters $`\alpha _k`$ $`(k=1,\mathrm{},N)`$ in (39), we cannot quote a value for them at the end of the analysis, but we can determine an ‘effective’ value $`\alpha _k^{\mathrm{eff}}(𝜽)`$ for each hyperparameter at any point in parameter space. In each case, this is given by the value of $`\alpha _k`$ that maximises the modified likelihood $`\mathrm{Pr}(\text{D}|𝜽,\alpha _k)`$ at that point. Differentiating the expression (32) with respect to $`\alpha _k`$ and setting the result equal to zero, we find

$$\alpha _k^{\mathrm{eff}}(𝜽)=\frac{n_k}{\chi _k^2(𝜽)}.$$
(44)
Clearly, this expression is most meaningful when evaluated at the point $`𝜽=\widehat{𝜽}`$.

Finally, it is worth noting the case in which minimisation of the functions (38) and (39) gives the same parameters estimates $`\widehat{𝜽}`$. We see that the expressions (40) and (41) are identical when $`\chi _k^2(\widehat{𝜽})=n_k`$ for all $`k`$. From (44), we see that this corresponds to the case in which $`\alpha _k^{\mathrm{eff}}(\widehat{𝜽})=1`$ for all $`k`$.

### 5.3 Error estimates and approximate evidences

The estimated covariance matrix of the parameter errors is given by (3). From (38) and (39), for the models $`H_0`$ and $`H_1`$ respectively, we have

$`\text{C}_0^1`$ $`=`$ $`\frac{1}{2}{\displaystyle \underset{k=1}{\overset{N}{}}}[\mathrm{ln}|\text{V}_k|+\chi _k^2]_{𝜽=\widehat{𝜽}_0},`$ (45)
$`\text{C}_1^1`$ $`=`$ $`\frac{1}{2}{\displaystyle \underset{k=1}{\overset{N}{}}}[\mathrm{ln}|\text{V}_k|+(n_k+2)\mathrm{ln}(\chi _k^2+2)]_{𝜽=\widehat{𝜽}_1}.`$ (46)
Since we may write

$$\mathrm{ln}(\chi _k^2+2)=\frac{(\chi _k^2+2)\chi _k^2(\chi _k^2)(\chi _k^2)^\mathrm{T}}{(\chi _k^2+2)^2},$$
(47)
we see that, if one can calculate the functions $`\mathrm{ln}|\text{V}_k|`$ and $`\chi _k^2`$ (together with $`\chi _k^2`$, which is required for obtaining the parameter estimates) to give the standard (inverse) covariance matrix $`\text{C}_0^1`$, one may easily calculate $`\text{C}_1^1`$.

Once the (inverse) covariance matrices of the errors have been calculated, the evaluation of (approximate) evidences is straightforward. Using the approximate expression (10) for the evidence, and assuming the uniform prior (28), we see that the ratio of evidences for the two models $`H_0`$ and $`H_1`$ is given by

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}\frac{|\text{C}_1|^{1/2}\mathrm{Pr}(\text{D}|\widehat{𝜽}_1,H_1)}{|\text{C}_0|^{1/2}\mathrm{Pr}(\text{D}|\widehat{𝜽}_0,H_0)}.$$
(48)
If this ratio is much less than unity, one concludes that the datasets are mutually consistent and do not support the introduction of hyperparameters; one should then quote the results of the standard analysis, given by the estimates $`\widehat{𝜽}_0`$ and the covariance matrix $`\text{C}_0`$. If the ratio is much greater than unity, however, then this indicates that inconsistencies do exist between (some of) the datasets and that the inclusion of hyperparameters is warranted; it is then more appropriate to quote the estimates $`\widehat{𝜽}_1`$ and the covariance matrix $`\text{C}_1`$.

## 6 A toy problem

To illustrate the potential usefulness of the hyperparameters approach to weighting datasets, in this section we apply the method to the classic parameter estimation problem of fitting a straight line through a set of data points. In particular, we will show how the use of hyperparameters as weights can overcome the common problems of inaccurately quoted error-bars and the presence of systematic errors in the measurements (see also Bridle 2000).

Let us suppose that the true underlying model for some process is the straight line

$$y(x)=mx+c,$$
(49)
where the slope $`m=1`$ and the intercept $`c=1`$. We assume that two independent sets of measurements $`\text{D}_1`$ and $`\text{D}_2`$ are made, each of which contains five data points, so $`n_1=5`$ and $`n_2=5`$. In simulating each dataset, the $`x`$-values are drawn at random from a uniform distribution between zero and unity, and we assume that these values are known precisely. The corresponding $`y`$-values are independently distributed about their true model values according to a Gaussian distribution of known variance $`\sigma _k`$ for each data set. Thus, the likelihood function for the $`k`$th dataset has the Gaussian form given by (29) and (30), where $`(𝝁_k)_j=m(\text{x}_k)_j+c`$ and $`\text{V}_k=\text{diag}(\sigma _k^2,\mathrm{},\sigma _k^2)`$. Thus, in this simple example, the covariance matrix does not depend on the parameters $`m`$ and $`c`$.

We now consider three different scenarios, which illustrate the general useful properties the hyperparameters technique for weighting datasets. We note that, for this simple toy problem, it is clear from a plot of the data points if inconsistencies exist between two datasets, and we will see that the hyperparameters approach verifies our initial suspicions where appropriate. Indeed, this toy model was chosen so that inconsistencies revealed by the hyperparameters method are transparent from a plot of the data. It should be remembered, however, that for realistic cosmological datasets and theoretical models, the situation is usually much more complicated. In this case, one is often unable visually to discern any inconsistencies in the data, but these may nevertheless be deduced by using the hyperparameters technique. This is illustrated in section 7.

### 6.1 Accurate error-bars and no systematic error

In our first case, both data sets $`\text{D}_1`$ and $`\text{D}_2`$ are drawn from the correct model $`m=1`$ and $`c=1`$. The noise rms for each dataset is $`\sigma _1=\sigma _2=0.1`$, and these (correct) values are used in the parameter estimation process. The two simulated datasets are shown in Fig. 1 (left panel), in which the dataset $`\text{D}_1`$ is denoted by the filled circles, whereas $`\text{D}_2`$ corresponds to the open squares. The true underlying model is plotted as a solid line.

Assuming an appropriately normalised uniform prior on the parameters in the range $`m:02`$ and $`c:02`$, and using the expressions (31) and (35) for the likelihood in each case, we may calculate the (unnormalised) posteriors $`\overline{\mathrm{Pr}}(m,c|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(m,c|\text{D},H_1)`$ corresponding to the standard and hyperparameters approaches respectively (where D denotes the combination of the datasets $`\text{D}_1`$ and $`\text{D}_2`$). These posterior distributions are plotted in Fig. 1 (middle and right panels), in which the contours contain 68, 95 and 99 per cent respectively of the total probability; the true parameter values are indicated by a solid circle.

In both cases $`H_0`$ and $`H_1`$, the posterior distributions appear very similar and contain the true parameter values at about the 68 per cent confidence level. However, since the prior and the likelihood function in each case are properly normalised, the greyscale units on the plots in Fig. 1 are not arbitrary. Indeed, the evidence in each case can be calculated directly by integrating numerically under these distributions, and we find

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}=0.54.$$

This indicates that standard approach is very slightly preferred or, equivalently, that the introduction of hyperparameters is marginally disfavoured by the data. An evidence ratio of order unity does not, however, provide a strong indication that one case is preferred over the other. We may conclude that, in this case, our ability to infer the parameter values $`m`$ and $`c`$ from the data is essentially unaffected by the introduction of hyperparameters.

It is interesting to compare the exact evidence ratio above with that obtained using the approximate expression (48), in which each posterior distribution is approximated by a Gaussian centred at its peak. Since the function (49) is linear in the parameters $`m`$ and $`c`$, and the covariance matrices $`\text{V}_k`$ do not depend on these parameter, the posterior $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ is, in fact, truly Gaussian. We found, in this case, that the approximate expression (10) for the evidence did indeed agree with that obtained by direct numerical integration. In the hyperparameters case, however, the posterior $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ is not Gaussian. Nevertheless, we found that the approximate expression (10) underestimated the true evidence value by only 2 per cent, and the approximate evidence ratio was found to be 0.53.

Finally, we note the effective weight assigned to each dataset at the peak $`\widehat{𝜽}_1`$ of the hyperparameters posterior $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$. These values are given by (44) evaluated at this point, and we find $`\alpha _1^{\mathrm{eff}}(\widehat{𝜽}_1)=1.79`$ and $`\alpha _2^{\mathrm{eff}}(\widehat{𝜽}_1)=2.64`$.

### 6.2 Inaccurate error-bars and no systematic error

In this case, the datasets $`\text{D}_1`$ and $`\text{D}_2`$ are identical to those used in the previous subsection, but we assume that the quoted errors-bars on the dataset $`\text{D}_1`$ have been underestimated by a factor of 5, while the error-bars on the dataset $`\text{D}_2`$ are quoted correctly. Thus, in the parameter estimation procedure, we assume the (incorrect) values $`\sigma _1=0.02`$ and $`\sigma _2=0.1`$. The data points and their reported error-bars are shown in Fig. 2 (left panel).

The resulting (unnormalised) posteriors $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ are shown in the middle and right-hand panels of the figure. In this case, the two posteriors are very different. In the standard approach $`H_0`$, the posterior distribution is tightly constrained about its maximum as a result of underestimating the errors on dataset $`\text{D}_1`$. Indeed, this posterior is virtually indistinguishable from that calculated from the dataset $`\text{D}_1`$ alone. The true parameter values are excluded at a confidence level that far exceeds the 99 per cent limit. In the hyperparameters case, however, the posterior is much broader and resembles the corresponding distribution in Fig. 1. In particular, the true parameter values are comfortably contained within the 95 per cent confidence limit.

Integrating under the posterior distributions directly, the exact evidence ratio is found to be

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}=2.1\times 10^4,$$

which clearly implies that the data favour the introduction of hyperparameter weights. Using the expression (48), the approximate evidence ratio is $`1.6\times 10^4`$, which shows that the Gaussian approximation to the hyperparameters posterior is once again reasonably accurate.

At the peak $`\widehat{𝜽}_1`$ of the hyperparameters posterior, we find the effective weights assigned to the two datasets $`\text{D}_1`$ and $`\text{D}_2`$ to be $`\alpha _1^{\mathrm{eff}}(\widehat{𝜽}_1)=0.12`$ and $`\alpha _2^{\mathrm{eff}}(\widehat{𝜽}_1)=1.28`$. Thus, the first dataset (with error-bars underestimated by a factor of 5) has been assigned an appropriate smaller statistical weight.

### 6.3 Accurate error-bars and a systematic error

We now consider the introduction of a systematic error into the dataset $`\text{D}_1`$. This is simulated by drawing this dataset from an incorrect straight-line model, for which the parameter values $`m`$ and $`c`$ differ from unity. Dataset $`\text{D}_2`$, however, is still drawn from the correct straight-line model with $`m=1`$ and $`c=1`$. We shall, in fact, consider two separate complementary cases. In the first case, we introduce a systematic error in the direction of the natural degeneracy line in the $`(m,c)`$-plane, whereas, in the second case, the introduced systematic error is orthogonal to the natural degeneracy line. In both cases, we assume that the error-bars on each dataset are quoted accurately as $`\sigma _1=0.1`$ and $`\sigma _2=0.1`$ respectively.

#### 6.3.1 Case 1

In our first case, dataset $`\text{D}_1`$ is drawn from a straight line with slope $`m=0`$ and intercept $`c=1.5`$. The datasets $`\text{D}_1`$ and $`\text{D}_2`$ are shown in Fig. 3 (left panel), together with the underlying straight-line models from which each is drawn. The resulting posterior distributions $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ are shown in the middle and right-hand panels of the figure. Once again, the two posteriors are very different. In the standard approach, $`H_0`$, the posterior distribution peaks between the two sets of true values. In spite of the fact that the two sets of parameter values define a direction along the natural degeneracy line in the $`(m,c)`$-plane, neither is contained within the 99 per cent confidence contour, and so both models are excluded at a high significance level. In the hyperparameters case, however, the posterior is much further extended along the natural degeneracy line. In particular, we note that the distribution is bimodal, with each peak lying close to one of the true sets of parameter values. Thus, the hyperparameters indicates the presence of two underlying models for the data, which signals an inconsistency between the two datasets. This could be interpreted as one (or both) of the datasets containing a systematic error.

The exact evidence ratio in this case is found to be

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}=11.6,$$

which gives a reasonably robust indication that the data favour the introduction of hyperparameters. Using the expression (48), the approximate evidence ratio is given by 4.8. The reason for the large inaccuracy in this case is that the Gaussian approximation to the bimodal hyperparameters posterior is clearly rather poor. In fact, as might be expected in this case, the Gaussian approximation underestimates the true value of the evidence by about a factor of two.

Although the hyperparameters posterior is bimodal, the global maximum $`\widehat{𝜽}_1`$ occurs at the peak close to the parameter values from which the dataset $`\text{D}_1`$ was drawn. At this peak, we find the effective weights assigned to the two datasets $`\text{D}_1`$ and $`\text{D}_2`$ are $`\alpha _1^{\mathrm{eff}}(\widehat{𝜽}_1)=3.11`$ and $`\alpha _2^{\mathrm{eff}}(\widehat{𝜽}_1)=0.18`$, which indicates (correctly) that the first dataset has been assigned a larger statistical weight at this point in parameter space. However, at the subsidiary peak $`\widehat{𝜽}_1^{}`$ located near the parameter values from which the dataset $`\text{D}_2`$ was drawn, we find $`\alpha _1^{\mathrm{eff}}(\widehat{𝜽}_1^{})=0.11`$ and $`\alpha _2^{\mathrm{eff}}(\widehat{𝜽}_1)=6.41`$, and so the roles of the datasets have been reversed.

#### 6.3.2 Case 2

In our second case, dataset $`\text{D}_1`$ is drawn from a straight line with slope $`m=0.7`$ and intercept $`c=0.7`$. The datasets $`\text{D}_1`$ and $`\text{D}_2`$ are shown in Fig. 4 (left panel), together with the underlying straight-line models from which each is drawn. The the middle and right-hand panels of the figure show the resulting posterior distributions $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$. As in Case 1, the standard approach produces a posterior distribution that peaks between the two true sets of parameters values, excluding both at a high significance level. We note, in this case, that the two sets of true parameter values define a direction orthogonal to the natural degeneracy line in the $`(m,c)`$-plane. In the hyperparameters case, however, the posterior is again bimodal, with each peak lying close to one of the true sets of parameter values. Thus, once again the hyperparameters approach signals the presence of two underlying models for the data and hence an inconsistency between the datasets.

The exact evidence ratio is found to be

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}=5.9\times 10^3,$$

which strongly implies that the data favour the introduction of hyperparameters. In this case, the approximate evidence ratio (48) is given by $`2.4\times 10^3`$. As in Case 1, the Gaussian approximation underestimates the true value of the evidence for the hyperparameters posterior by about a factor of 2, as a result of it being bimodal.

Once again the effective weights clearly show which dataset is dominating at each peak of the hyperparameters posterior. At the global maximum $`\widehat{𝜽}_1`$, we find $`\alpha _1^{\mathrm{eff}}(\widehat{𝜽}_1)=0.066`$ and $`\alpha _2^{\mathrm{eff}}(\widehat{𝜽}_1)=7.96`$, whereas at the subsidiary maximum $`\widehat{𝜽}_1^{}`$ we obtain $`\alpha _1^{\mathrm{eff}}(\widehat{𝜽}_1^{})=3.31`$ and $`\alpha _2^{\mathrm{eff}}(\widehat{𝜽}_1)=0.079`$.

## 7 A cosmological illustration

In this section, we illustrate the hyperparameters technique for weighting datasets by applying it to recent measurements of the CMB power spectrum provided by the BOOMERANG (Netterfield et al. 2001), MAXIMA (Lee et al. 2001) and DASI (Halverson et al. 2001) experiments. We also include the 8 data points for large-scale normalisation of the spectrum provided by the COBE satellite. The flat band-powers reported by each of these experiments are plotted in Fig 5.

The solid line in the plot corresponds to the predicted CMB power spectrum for a spatially-flat inflationary CDM model with no tensor contribution and $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_bh^2=0.02`$, $`h=0.7`$, $`n=1`$, $`\tau =0`$ and $`Q=18`$ $`\mu `$K.

From Fig. 5, we see immediately that there is good agreement between the datasets, all of which are broadly consistent with the plotted theoretical CMB power spectrum. One would therefore expect the inclusion of and marginalisation over hyperparameter weights not to be warranted by the data. It should also be remembered, however, that the points and error-bars plotted in the figure do not include the calibration and beam uncertainties associated with each experiment. Indeed, since these uncertainties can be as large as 20 per cent in some cases, one would naively expect the envelope of the experimental data to be much wider, resulting in poorer agreement with the theoretical spectrum. As we will see below, this naive observation is supported by the proper calculation of Bayesian evidences.

### 7.1 No calibration or beam uncertainty

We first analyse the data without including calibration and beam errors, i.e. as plotted in Fig. 5. For each dataset $`\text{D}_k`$ $`(k=1,2,3,4)`$, we assume the uncertainties on the flat band-powers are described by a multivariate Gaussian of the form given in (29). For the COBE and DASI datasets, we use the publicly-available window functions and full covariance matrix $`\text{V}_k`$ for each dataset. In the absence of the equivalent information for the BOOMERANG and MAXIMA experiments, we assume a top-hat window function for the spectrum in each bin and neglect correlations between bins. This approach is a good approximation of the correct one (see de Bernardis et al. 2001) and does not affect our conclusions.

Denoting the totality of the resulting data by the vector D, the next step in the analysis is to calculate the unnormalised posterior distributions $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$, corresponding to the standard and hyperparameters approaches respectively, for some set of cosmological parameters $`𝜽`$. For the purposes of illustration, we assume a spatially-flat Universe ($`\mathrm{\Omega }_k=1`$), with no tensor contribution to the CMB spectrum, and take the parameter vector $`𝜽`$ to consist of 5 parameters, namely the normalisation $`Q`$, the Hubble parameter $`h`$, the matter density $`\mathrm{\Omega }_\mathrm{m}`$, the physical baryon density $`\mathrm{\Omega }_\mathrm{b}h^2`$ and the scalar spectral index $`n`$. Since the parameter space is only 5-dimensional, it is feasible to evaluate the unnormalised posteriors over a hypercube. Thus, we assume (suitably normalised) uniform priors on the parameters, with the ranges $`12<Q<24`$ $`\mu `$K, $`0.5<h<0.9`$, $`0.1<\mathrm{\Omega }_\mathrm{m}<0.5`$, $`0.01<\mathrm{\Omega }_\mathrm{b}h^2<0.04`$ and $`0.8<n<1.2`$. The corresponding likelihood functions for the hypotheses $`H_0`$ and $`H_1`$ have the forms (31) and (35) respectively. The posterior distribution is calculated at 20 points in the directions $`Q`$, $`h`$ and $`\mathrm{\Omega }_m`$ in parameter space, and at 10 points in the directions $`n`$ and $`\mathrm{\Omega }_bh^2`$. The theoretical power spectrum corresponding to each model was calculated using CAMB (Lewis, Challinor & Lasenby 2000).

Since $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ are each calculated over a hypercube in the 5-dimensional parameter space, one can calculated the evidence integrals (36) and (37) directly, and one finds the exact evidence ratio to be

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}=0.05.$$

This indicates that the data do not support the inclusion of and subsequent marginalisation over a free hyperparameter weight $`\alpha _k`$ $`(k=1,2,3,4)`$ for each dataset. In other words, as expected, the datasets are statistically consistent both with one another and with the range of theoretical models with which they have been compared.

For illustration, let us consider the constraints placed by the current CMB data on the parameters $`n`$ and $`\mathrm{\Omega }_bh^2`$ with and without the inclusion of hyperparameter weights. The limits placed on these parameters by the standard and hyperparameters approaches respectively are easily obtained by marginalising the corresponding posteriors $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ over the remaining parameters $`Q`$, $`h`$ and $`\mathrm{\Omega }_m`$. The resulting distributions are shown in Fig. 6, together with corresponding the 68 and 95 per cent confidence contours. Although the above evidence ratio shows that the data do not favour marginalisation over hyperparameter weights, we see that their inclusion does not significantly affect the constraints placed on the parameters. This illustrates that, for statistically consistent datasets, the standard and hyperparameters techniques may be use interchangeable without degrading the constraints on cosmological parameters. In particular, we note that the current CMB data are consistent with the scale invariant spectrum $`n=1`$ and the nucleosynthesis constraints on $`\mathrm{\Omega }_bh^2`$. Indeed, by performing an additional marginalisation over $`n`$ or $`\mathrm{\Omega }_bh^2`$ we obtain (to two significant figures) the one-dimensional 68-per cent confidence intervals $`0.019<\mathrm{\Omega }_bh^2<0.023`$ and $`0.92<n<1.00`$ for both the case $`H_0`$ and $`H_1`$.

### 7.2 Marginalisation over calibration and beam uncertainties

So far our analysis has ignored the calibration and beam uncertainties associated with each dataset. Any meaningful analysis of these data must, however, include these effects of these ‘nuisance’ parameters. As with hyperparameter weights, the correct approach is to assign some prior to these parameters and then marginalise over them. For calibration and beam uncertainties, one does have a priori knowledge of both the expectation value of each parameter and the range of values it might take (i.e. its variance). Thus, as shown in section 4.1, it is appropriate to adopt independent Gaussian priors on the calibration and beam uncertainties. As discussed in detail by Bridle et al. (2001), for Gaussian likelihood functions of the form (29), one can again perform this marginalisation analytically. Moreover, the resulting likelihood function after marginalisation is also of the form (29), but with a modified covariance matrix $`\text{V}_k^{}`$.

Using these modified datasets, one can then calculate the posteriors $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_0)`$ and $`\overline{\mathrm{Pr}}(𝜽|\text{D},H_1)`$ over the same hypercube in the 5-dimensional parameter space as used above, and evaluate the evidence integrals (36) and (37) directly. In this case, the evidence ratio is found to be

$$\frac{\mathrm{Pr}(\text{D}|H_1)}{\mathrm{Pr}(\text{D}|H_0)}=0.08,$$

which illustrates that the inclusion of and marginalisation over hyperparameter weights is once more unwarranted. The corresponding marginalised posteriors in the $`(n,\mathrm{\Omega }_bh^2)`$-plane, with and without marginalisation over are again almost identical and closely resemble those plotted in Fig. 6. One obtains the one-dimensional 68-per cent confidence intervals $`0.018<\mathrm{\Omega }_bh^2<0.024`$ and $`0.90<n<0.99`$ for the case $`H_0`$ (without hyperparameter weights) and $`0.018<\mathrm{\Omega }_bh^2<0.024`$ and $`0.90<n<1.00`$ for the case $`H_1`$ (with hyperparameter weights). Once again, we see that the inclusion of and marginalisation over hyperparameter weights has little effect on the parameter constraints imposed by the data.

### 7.3 Evidence for calibration and beam uncertainties

Although the evidence ratios given above show that the current CMB power spectrum data clearly do not require the use of hyperparameter weights in their analysis, it is also interesting to compare the relative evidence values for the cases with and without marginalisation over calibration and beam uncertainties (for both $`H_0`$ and $`H_1`$). The evidence values corresponding to each of the 4 cases we have considered are given in Table 1, where they have been rescaled so that the evidence is unity for the case where no hyperparameter weights are introduced and no marginalisation over calibration and beam uncertainties is performed. We see that the largest evidence value is obtained for the case in which no marginalisation is performed over hyperparameter weights or calibration and beam uncertainties. While one might expect from the data plot in Fig. 5 that the introduction of and marginalisation over hyperparameter weights might not be favoured, it is surprising that the marginalisation over the calibration and beam uncertainties also reduces the evidence for the data. This means that, when compared with the range theoretical models discussed above, the probabilty of obtaining the observed data is highest if one assumes no uncertainty exists in the calibration and beam properties adopted by each experiment, rather than marginalising over the reported uncertainties in these parameters. Thus we see that the evaluation of evidences has provided quantitative support for our earlier naive observation that the CMB datasets plotted in Fig. 5 agree remarkably well given the 5–20 per cent calibration uncertainties alone (which are not plotted). We note here simply that this eventuality is unlikely to occur by chance.

## 8 Conclusions

We have presented a general account of the use of hyperparameters in the analysis of cosmological datasets. In particular, we have concentrated on applying the hyperparameters technique to the problem of weighting different datasets in a joint analysis aimed at estimating some set of cosmological parameters. The basic approach is to assign a free hyperparameter weight to each dataset and then perform a marginalisation over these hyperparameters. This method allows the statistical properties of the data themselves to determine the effective weight given to each dataset and is in sharp contrast to the more common practice of excluding certain datasets from the analysis, which hence are assigned a weight of zero, and analysing the remaining datasets with equal weights.

Assuming the expected value of the hyperparameter weight on each dataset to be unity, we find that the prior on each hyperparameter should be of exponential form, rather than the more commonly used improper Jeffrey’s prior. In either case, the marginalisation over the hyperparameters may be performed analytically for Gaussian likelihood functions. Since the exponential prior is properly normalised, in this case one may also calculate the Bayesian evidence for the data given the hyperparameters model. This evidence value can then be compared with the corresponding evidence for the data in the absence of hyperparameters, in order to determine whether the data warrant the inclusion of hyperparameters in the first instance. In each case, the evidence may be calculated either by direct integration, for low-dimensionality parameter spaces, or approximated straightforwardly by assuming each posterior to be Gaussian near its peak.

The hyperparameter approach to weighting datasets is illustrated by applying it to the classic toy model of fitting a straight to a number of datasets. We find the hyperparameters technique correctly infers the existence of systematic errors and/or misquoted random errors in the datasets. In such cases, the evidence ratio for the hyperparameters and standard approaches clearly indicates that the data warrant the introduction of the hyperparameters. Nevertheless, in the case where no systematic errors exists, and the random errors are accurately quoted, the evidence ratio correctly indicates that hyperparameters should not be included in the analysis.

Finally, the hyperparameters technique is applied to the latest measurements of the cosmic microwave background (CMB) power spectrum by the BOOMERANG, MAXIMA and DASI experiments, together with the large-scale normalisation of the power spectrum provided by the COBE satellite. The evidence ratio between the hyperparameters and standard approach shows that the current CMB datasets do not warrant the inclusion of hyperparameters weights and a subsequent marginalisation over them. In other words, the datasets are statistically consistent both with one another and with the range of theoretical models with which they were compared. Nevertheless, the inclusion of hyperparameters in shown not to lessen the constraints imposed by the data on cosmological parameters.

The analysis is repeated for the case in which one first marginalises over the beam and calibration uncertainties associated with each dataset, and again the evidence ratio shows that the inclusion of hyperparameter weights is unwarranted. Of more interest is the comparison of evidences with and without marginalisation over calibration and beam uncertainties. This shows that the evidence for the data is largest if one assumes no uncertainty in the calibration and beam properties, but instead adopt the values corresponding to those implied by the published CMB power spectrum points for each dataset. This is unlikely to occur by chance.

### ACKNOWLEDGMENTS

MPH thanks Steve Gull and Anthony Lasenby for some illuminating discussions concerning Bayes’ theorem and probability theory. The authors also thank Antony Lewis for providing the $`C_{\mathrm{}}`$-grid and Anze Slozar for the use of his $`C_{\mathrm{}}`$-grid interpolation routines. SLB acknowledges Selwyn College, Cambridge, for support in the form of the Trevelyan Research Fellowship.


format: text

text:
# Remnants of the Sagittarius Dwarf Spheroidal Galaxy around the young globular cluster Palomar 12

## 1 Introduction

The picture of building the galactic halo from merging “fragments” resembling dwarf galaxies, which Searle & Zinn (1978, hereafter SZ) proposed on the basis of the properties of the Milky Way globular clusters, is often regarded as the local manifestation of the hierarchical galaxy formation theory whereby dwarf galaxies were the first to form and their subsequent mergers created large galaxies (e.g. Moore et al. 1998; Navarro, Frenk & White 1995). The Milky Way continues to be one of the best places to test this picture, and much of the recent work in this area has focused on the identification of stellar streams left behind in the halo after the accretion of low-mass satellite galaxies (e.g. Yanny et al. 2000; Ibata et al. 2001, Dohm-Palmer et al. 2001; Vivas et al. 2001). The prototypical case of accretion is the present-day tidal destruction of the Sgr dwarf spheroidal (dSph) galaxy, which has deposited at least 4 and probably 5 (Dinescu et al. 2000) globular clusters in the halo and produced a long stellar stream that wraps around the sky. One of the Sgr clusters, the very luminous cluster M54, and the most luminous globular cluster in the Galaxy, $`\omega `$Cen, are suspected to be the nucleus of Sgr (Sarajedini & Layden 1995) and a now extinct dwarf galaxy, respectively (Lee et al. 1999). van den Bergh (2000) has argued that the “young” globulars located in the outer halo might be the nuclei of extinct dSph galaxies and has suggested that searches be made for the vestiges of their parent systems. We report here our survey of a large field surrounding Palomar 12 (Pal 12)

Pal 12 is a remote globular cluster located at a distance of 19 kpc from the Sun. It is among the youngest and most metal rich of the halo globular clusters, which has fueled speculation that has been accreted from another galaxy (e.g., Zinn 1993). On the basis of its age and radial velocity, Lin & Richer (1992) argued that Pal 12 may have been captured from the Magellanic Clouds. The measurement of its proper motion and the determination of its orbit around the Milky Way by Dinescu et al. (2000) have shown that Pal 12 probably did not originate in the Clouds but in the Sgr dSph galaxy (as suggested first by Irwin 1999). In this paper, we report the detection of a very low density stellar system in the same direction and, to within the errors, the same distance as Pal 12.

## 2 OBSERVATIONS AND DATA REDUCTION

Pal 12 was observed in $`B`$ and $`R`$ Johnson–Cousins filters with the Wide Field Camera (WFC) at the prime focus of the 2.5 m Isaac Newton Telescope (INT) at the Roque de los Muchachos Observatory in June 2001. The WFC holds four 4096 $`\times `$ 2048 pixel EEV CCDs with pixel size $`0\stackrel{}{\mathrm{.}}33`$, which provides a total field of about $`35\mathrm{}\times 35\mathrm{}`$. The field was centered at 10$`\mathrm{}`$ South of the center of Pal 12 with the purpose of including the extra-tidal regions around the cluster (Field 1). In addition, a field situated 2$`\mathrm{°}`$ North of the cluster’s center was also taken with the purpose of gauging the field-star contamination (Field 2). However, this field is not far enough away to avoid the large stellar system that we have discovered in the environs of Pal 12. Another control field at Northern galactic latitute, which was observed during the same run, has been used instead.

Bias and flatfield corrections were made with IRAF. DAOPHOT and ALLSTAR (Stetson 1994) were then used to obtain the instrumental photometry of the stars. For the final photometric list, we selected stars with $`\sigma <0.20`$, $`1<SHARP<1`$ and $`0<CHI<2`$ as provide by ALLSTAR. These criteria reject extended objects, so the background contamination is expected to be only stellar-shaped objects. The atmospheric extinction and the transformations to the standard Johnson–Cousins photometric system were obtained from 52 measurements of 20 standard stars from the Landolt (1992) list. The photometric conditions during the observing run were stable and produced very small zero point errors for the photometric transformation: $`\pm `$ 0.020 mag in $`B`$ and $`\pm `$ 0.016 in $`R`$. The offsets between the different WFC chips were obtained from observations of several standard fields in each chip, and then removed with a precision of better than 0.01 mag. Taking all uncertainties into account, we estimate the total zero-point errors of the photometry are $`0.025`$ in both filters.

## 3 THE LOW SURFACE BRIGHTNESS SYSTEM

Figure 1a shows the color-magnitude diagram (CMD) for the WFC field centered on Pal 12 (Field 1). In addition to the main-sequence (MS) of the cluster, a MS-like feature is observed at bluer color, overlapping the MS of the cluster at $`BR0.8`$ and $`V>21.0`$. This feature is also clearly seen in the CMD of the field situated $`2\mathrm{°}`$ North of the center of Pal 12 (Field 2).

This unexpected feature of the Pal 12 diagram is more evident in Figure 1b, which shows the CMD of the extra-tidal field of Pal 12. This includes the part of Field 1 beyond the tidal radius of the cluster (i.e. for $`r>17\mathrm{}`$) and one-half of the area of Field 2. These regions have almost indistinguishable CMDs, with the only difference that the density of stars in the MS feature is larger in the cluster’s field (Field 1). The CMD of the control field situated at (l,b)=$`(28.7,42.2)`$ in the North is shown in Figure 1c. The total area of this field is the same as that plotted in Fig. 1b. The absence of the MS feature in the control field is the most striking difference between these Northern and Southern hemisphere regions. The control field has slightly smaller l and $`b`$ than the fields in the direction to Pal 12. If the halo is symmetric about the galactic plane and has smooth density contours, either spherical or flattened towards the plane, the control field and not the fields in the direction to Pal 12 should have a higher density of halo stars. This is clearly not the case (compare the regions bordered by $`22>^{\prime \prime }V^{\prime \prime }>18`$ & $`0.6<BR<2.0`$ in fig. 1b and 1c).

To check the significance of this detection, we constructed luminosity functions of the stars in the “V” range 18.0-22.4, where the photometric errors are $``$ 0.05 and $``$ 0.03 for the field near Pal 12 and the control field, respectively. The color range was restricted to $`0.6`$B-R$`1.1`$, which should encompass the MS and the subgiant branch of an old stellar system. Its large width ensures that this comparison is insensitive to reddening differences, which are expected to be very small ($`\mathrm{\Delta }`$E(B-V)$`0`$, Schlegel et al. 1998, SFD). Figure 2 shows that for “V”$`<20.6`$ there is little difference between the field near Pal 12 and the control field. Hovever, for $`20.8`$“V”$`22.4`$, there is an excess of stars in the field near Pal 12 over the control field, and for every 0.2 mag. bin in this interval, this excess corresponds to several standard deviations.

Figure 2 also compares the luminosity functions with that of a model stellar population, which was constructed from the isochrone calculations of Girardi et al. (2000) for a metal abundance of Z=0.001 (\[Fe/H\] =-1.3) and an age of 12.6 Gyr. We imposed the same color boundaries as above and transformed to apparent “V” magnitude by adopting the same distance modulus and interstellar extinction as Pal 12. To model the field stars, we fit a quadratic equation to the luminosity function of the stars in the control field. The solid line in Figure 2 is the sum of this function plus the one for the 12.6 Gyr stellar population. It is normalized so that it matches the observed number of stars in the field near Pal 12 for “V”$`<19.0`$. The number of stars in the 12.6 Gyr system was varied until a rough match was obtained between the model function and the observed points for “V”$`>21.5`$. The similarity between this model luminosity function and the observed one over the whole range of “V” suggests that the excess of stars over the control field is caused by a real stellar system at approximately the same distance as Pal 12. Unfortunately, the number of stars is too small to tightly constrain the luminosity function, and firm conclusions about the age or composition of the system cannot be drawn from this match.

To estimate the surface brightness of this system, we used the control field for an estimate of the field star contamination and produced a decontaminated CMD following the procedure of Gallart, Aparicio & Vílchez (1996). Over the range $`20.5`$“V”$`22.9`$ and within the same color range as above, there is an excess of $`293\pm 15`$ MS stars over the background. The comparison of this number with that found in a decontaminated CMD of the main body of the Sgr dSph ( Martínez-Delgado et al. 2002) yields an estimate of $`\mathrm{\Sigma }`$= 30.5$`\pm `$ 0.2 mag arcsec-2. While undoubtedly this system has a very low surface brightness (LSB), we caution that this value is a rough estimate and also note that there is probably a gradient in $`\mathrm{\Sigma }`$ across the field. This value is the same found by Mateo et al. (1998) in their outermost detection of the Sgr stream, 34$`\mathrm{°}`$ SE from the center of Sgr.

A comparison between Figures 1a and 1b shows that the extra-tidal stellar population is bluer than that of Pal 12 and displays a significant width in color, indicating the presence of a range of age and/or metallicity and/or some depth along the line of sight, as expected of a dwarf galaxy. The simplest explanation is that this MS feature is not due to Pal 12 but to a separate stellar population that is part of a dSph galaxy or a tidal stream from one. The presence of this feature throughout the area covered by Field 2 puts a lower limit of 0.9 kpc on the extension of the system, assuming that it is at the same distance as Pal 12 (19.0 kpc; see below). This size is compatible with the width of a tidal stream crossing the field (e.g. Sgr tidal stream; Newberg et al. 2002) or the projected tidal radius of a typical dSph at the cluster’s distance.111 A typical dSph with a tidal radius of $``$ 1 kpc would have angular radius of 3 degrees at this distance..

While the comparison in Figure 2 between the model and observed luminosity functions illustrates that the LSB system could be at the same distance as Pal 12, it is important to examine this question in more detail. In Figure 1b, the MS appears to terminate at “V”=$`20.48\pm 0.10`$. The MS turnoff (i.e., the bluest point of the MS) is probably fainter than this (e.g., the turnoff of the 12.6 Gyr population in Fig. 2 occurs at 20.65, $`M_V=4.14`$ Girardi et al. 2000), but this point cannot identified with certainty given the small number of stars. It seems unlikely that it could be fainter than “V”=20.9. Adopting this value, a somewhat brighter turnoff luminosity ($`M_V=3.9`$), and correcting for the interstellar extinction (E(B-V)=$`0.037\pm 0.002`$, SFD), we obtain 23.8 kpc for an upper limit on the distance, if the LSB system is composed of a very old stellar population. For a lower limit under the same assumption, we adopt 20.48 as the turnoff and a fainter turnoff luminosity ($`M_V=4.2`$), which yields a distance of 17.1 kpc. It is likely that the LSB population is a mixture of ages and compositions, which is much harder to model. Under the assumption that its stellar population resembles the main body of the Sgr dSph galaxy, we can use the decontaminated CMDs of the LSB system and Sgr (see above) to help set the distance modulus of the LSB system. These CMDs suggest that the MS feature seen in Figure 1b is only the densest part of the MS and that its “termination” at “V”=20.48 corresponds to point in Sgr where $`M_V3.7`$ (Layden & Sarajedini 2000). This yields a distance of 21.6 kpc for the LSB system, with an uncertainty of $`2`$ kpc. The distance to Pal 12 is $`19.0\pm 0.9`$ kpc (Rosenberg et al. 1998, but using E(B-V) from SFD), which places it near the middle of the range of estimates for the distance to the LSB system. If they are unrelated, this is a very remarkable coincidence.

Dinescu et al. (2000) have shown on the basis of the orbit of Pal 12 that the Sgr dSph galaxy is likely to be its parent galaxy. Recent surveys of Sgr (Mateo, Olszewski & Morrison 1998; Ibata et al. 2001; Majewski et al. 1999; Yanny et al. 2000; Martínez-Delgado et al. 2001; Martínez-Delgado, Gómez-Flechoso & Aparicio 2002) have shown that this galaxy forms a giant stream that wraps completely around the Milky Way in an almost polar orbit, in good agreement with the predictions of theoretical models (Gómez-Flechoso, Fux & Martinet, 1999; Johnston, Sigurdsson & Hernquist 1999 ; Helmi & White 2001). The position of Pal 12 on the sky ($`40\mathrm{°}`$ from the Sgr’s main body) is very close to that predicted for the Sgr Southern stream (see Martínez-Delgado et al. 2001). The predicted distance for the Sgr’s tidal stream at the cluster’s position is $`19`$ kpc (Martínez-Delgado et al. 2002), in good agreement with the distance obtained above for th LSB system.

The CMD of the LSB system also resembles those of previous detections of the Sgr Southern stream. Majewski et al. (1999) reported a MS turnoff at V=21 in a field situated at (l,b)= (11$`\mathrm{°}`$, -40$`\mathrm{°}`$), which is $`20\mathrm{°}`$ from Pal 12. In the Sloan Digitized Sky Survey (SDSS), Newberg et al. (2002) found a similar structure (named S167-54-21.5) in the CMD of a long, narrow region centered on the celestial equatorial and $`15\mathrm{°}<\alpha <0\mathrm{°}`$, that could be also part of the Sgr tidal stream.

As we noted above, the MS of the LSB system is bluer than the MS of Pal 12 and extends to approximately the same bright magnitude. These characteristics suggest that it is more metal poor than Pal 12 ($`[Fe/H]`$=-1.0; Brown, Wallerstein & Zucker 1997) and may be as old or older than it, since lower metallicity isochrones have brighter turnoffs for a given age. The MS feature has intermediate color between the $`[Fe/H]`$=-1.7 and -0.7 isochrones of Bertelli et al. (1994), which yields a mean value of $`1.2`$. This is similar to the metallicity of the metal-poor component of the main body of Sgr ($`[Fe/H]`$= -1.3; Layden & Sarajedini 2000) and the giants in the Northern stream of Sgr (Dohm-Palmer et al. 2001).

## 4  DISCUSSION

Our detection of a LSB stellar system around Pal 12 strengthens considerably the contention of Dinescu et al. (2000) that Pal 12 originated in the Sgr dSph galaxy. The disruption of Sgr has therefore released at least 5 globular clusters to the Galactic halo (Ibata, Gilmore & Irwin 1994; Irwin 1999), and additional ones may remain undetected among the globular clusters in the outer halo. Two of the Sgr clusters (M54 and Ter 8) are similar in age to the oldest globular clusters in the Milky Way. A third (Arp 2) is slightly younger than these other two (Buonanno et al. 1998; Layden & Sarajedini 2000). The much younger Sgr cluster (Ter 7) and Pal 12 (Buonanno et al. 1998; Rosenberg et al. 1998) are several Gyrs younger than the other three Sgr clusters, and they are the youngest and most metal-rich globular clusters in the galactic halo. Sgr resembles the Fornax dSph galaxy in having 5 globular clusters, although their histories of cluster formation are significantly different. Only one of the Fornax clusters is clearly younger than the other four, and this cluster (cluster 4, Buonanno et al. 1999) is both older and less metal rich than the youngest Sgr clusters (Ter 7 and Pal 12).

Spectroscopic observations by Brown et al. (1997) have shown that Pal 12 is almost unique among halo globular clusters by having $`[\alpha /Fe]`$=0.0 instead of $`0.3`$. This relatively low $`\alpha `$ abundance is a sign of metal enrichment by Type Ia supernovae (Brown et al. 1997). The measurement of the abundance ratios in several red giants in Sgr by Smecker-Hane & McWilliam (1999) have shown that two metal-poor stars ($`[Fe/H]<1`$) are $`\alpha `$ enhanced while 9 others with $`[Fe/H]>1`$ have $`[\alpha /Fe]0.0`$. The chemical composition of Pal 12 ($`[Fe/H]=1`$, $`[\alpha /Fe]`$=0.0, Brown et al. 1997) is also consistent with membership in Sgr.

The three “young” Sgr clusters (Arp 2, Ter 7, Pal 12) constitute one quarter of the “young” globular clusters in the Galactic halo (van den Bergh 2000). Many of the other young halo clusters (perhaps all), and at least some of the old clusters may have originated in dwarf galaxies that later merged with the Milky Way. The stellar streams from Sgr show that this process not only added globular clusters to the halo, but also contributed stars of different ages and metallicities, as predicted by the SZ scenario.

We are grateful to M. A. Gómez-Flechoso, for making her Sgr model available to us and for several fruitful suggestions. We also thanks to R. García-López, Steve Majewski and M. Shetrone for many useful comments. This work is based on observations made with the 2.5 m INT operated by the Isaac Newton Group in the Spanish Observatorio del Roque de Los Muchachos of the Instituto de Astrofísica de Canarias. RZ was supported by grant AST 00-98428 from the National Science Foundation.

Figure 1. panel a) Color-magnitude diagram for the WFC field centered on Pal 12 (Field 1). The narrow main-sequence (MS) of the cluster is clearly delineated. Notice that it is surrounded by a sparser and wider MS-like feature at bluer color; panel b) displays the CMD for an extra-tidal field of $`35\mathrm{}\times 35\mathrm{}`$, half of which is from Field 1 but outside of the cluster’s tidal radius ( $`r>17\mathrm{}`$) and the other half is from Field 2, situated $`2\mathrm{°}`$ North of Pal 12. The MS feature is most clearly observed at $`BR0.8`$ and $`21<V<23`$; panel c) shows the CMD of our control field, which is the same size as the field shown in 1b.

Fig. 2. Comparison of the luminosity functions of the extra-tidal field of Pal 12 and the control field. The solid line is the model LF for a stellar population with age 12.6 Gyr and Z=0.001, adopting the distance and reddening of Pal 12. The error bars are $`\pm \sqrt{N}`$, where $`N`$ is the number of stars in each bin.