text:
# Implications for unified schemes from the quasar fraction and emission-line luminosities in radio-selected samples

## 1 Introduction

It has been recognised for some time that there is a strong positive correlation between the extended-radio luminosities and narrow-emission-line luminosities of 3C radio sources (Baum & Heckman, 1989; Rawlings et al., 1989; McCarthy, 1993), and by combining the 3C sample with the 7C Redshift Survey (7CRS), Willott et al. (1999) found that these correlations were not primarily due to redshift. This suggests that the sources of narrow-line emission and radio jets are linked, with explanations for this link ranging from the effects of environment (Dunlop & Peacock, 1993), to jet-disk symbiosis (Rawlings & Saunders, 1991; Falcke & Biermann, 1995).

Unified schemes for radio galaxies and radio quasars propose that they are the same objects viewed at different angles between their radio axis and their line of sight. An obscuring torus is invoked to hide the nucleus at large angles to the jet axis i.e. in radio galaxies. Emission from narrow-line regions is believed to be broadly independent of the jet-axis orientation, as it is emitted from a region larger than the extent of the torus. The opening angle of the torus $`\mathrm{\Theta }_{\mathrm{trans}}`$ marks the transition from the object being viewed as a quasar to a radio galaxy. We therefore expect for simple, arguably naive, unified schemes, where $`\mathrm{\Theta }_{\mathrm{trans}}`$ is constant, that the distribution of emission-line luminosities should be similar for radio quasars and radio galaxies. Jackson & Browne (1990) found that \[OIII\] luminosities of radio quasars are $`510`$ times more luminous than radio galaxies. In similar comparisons of radio quasars and radio galaxies, Hes, Barthel & Fosbury (1993) found no difference using \[OII\], and Jackson & Rawlings (1997) found no difference in \[OIII\] at high redshift. In summary, some studies comparing the narrow emission line strengths of radio quasars and radio galaxies are seemingly in agreement with the predictions of the simplest unified schemes, whilst others are not.

These seemingly contradictory, results can perhaps be understood in the context of a ‘receding torus’ model for AGN. This model (Lawrence, 1991; Hill, Goodrich & DePoy, 1996) proposes that the inner radius $`r`$ of the obscuring torus is determined by the radius at which dust sublimates, scaling as $`rL_{\mathrm{phot}}^{0.5}`$. Assuming that the half-height $`h`$ of the torus is independent of $`L_{\mathrm{phot}}`$, then the half-opening-angle of the torus, $`\theta =\mathrm{tan}^1r/h`$, will be larger in the higher-$`L_{\mathrm{phot}}`$ objects. This means that more luminous objects, with higher $`L_{\mathrm{phot}}`$, are more likely to be viewed as quasars so that quasar fraction rises systematically with $`L_{\mathrm{phot}}`$. It also means that orientation-independent quantities which scale with $`L_{\mathrm{phot}}`$ will be higher on average for quasars than for radio galaxies. Lawrence (1991) found that the fraction of narrow-line objects in the low-frequency-selected 3CR sample decreased with increasing radio luminosity, and that narrow-line objects have weaker \[OIII\] at a fixed radio luminosity. He argued that this was inconsistent with the simplest unified schemes but that it could be explained by a cone-angle dependence on luminosity, e.g a receding-torus model. Simpson (1998) argued that the \[OIII\] emission line is a much better indicator of $`L_{\mathrm{phot}}`$ than \[OII\], explaining why differences between quasars and radio galaxies are more obvious in \[OIII\] than \[OII\]. Simpson (2003) revisited the arguments of Simpson (1998) and corrected a small error (compare Fig. 2 of Simpson, 2003 with Fig. 4 of Simpson, 1998). He concluded that, allowing for a spread in $`h`$, the receding torus model predicts that radio quasars should be a factor of a few brighter in $`L_{\mathrm{phot}}`$ than radio galaxies in samples exhibiting a wide range of quasar fractions.

Willott et al. (2002a) found that radio quasars have higher submillimetre luminosities by a factor of $`4`$ than radio galaxies of the same radio luminosity and redshift, a factor which cannot be reduced below $`2`$ by accounting for possible synchrotron contamination. This result supports the idea that the simplest unified schemes, where $`\mathrm{\Theta }_{\mathrm{trans}}`$ is constant, are not an adequate description of the relationship between radio galaxies and quasars: submillimetre emission comes from cool dust grains in optically thin regions, and therefore radiates isotropically; moreover, such emission could not be obscured by a torus even if it were emitted from regions close to the nucleus. Submillimetre luminosity is therefore an orientation-independent quantity which might scale closely with $`L_{\mathrm{phot}}`$ either because the cool dust is heated directly by the quasar or because it is heated by a starburst whose luminosity scales with $`L_{\mathrm{phot}}`$. We conclude that the submillimetre study of Willott et al. (2002a) is in quantitative agreement with the receding torus model of Simpson (2003).

Other arguments for a receding-torus-like model include the evidence that the fraction of lightly reddened $`z1`$ 3C quasars decreases with increasing radio luminosity in agreement with the higher fraction of lines of sight expected to graze the torus at lower $`L_{\mathrm{phot}}`$ (at fixed $`h`$) in the receding torus model (Hill et al., 1996; Simpson, Rawlings & Lacy, 1999). Also the near-infrared spectral indices of quasars from the 3CR sample are correlated with luminosity, whereas the optical spectral indices are uncorrelated with the quasar luminosity or orientation, so that the strength of the ‘big red bump’, relative to the ionizing continuum, appears to be less in the more luminous objects (Simpson & Rawlings, 2000; Simpson, 2003).

However, there is significant scatter in the relationship between radio and emission-line luminosities and a receding torus is not the only way of explaining the differences in the emission-line properties of radio quasars and radio galaxies or the luminosity-dependence of the quasar fraction (the fraction of objects that show broad emission lines). Two-population radio luminosity functions (RLFs) have been used to provide a best-fit to the 3CRR, 6CE and 7CRS radio source redshift surveys and radio source counts (Willott et al., 2001). It is possible that a two-population model with a simple unified scheme in one population, combined with the effects of scatter, could mimic the effects of a receding torus in producing both emission-line differences between radio galaxies and radio quasars and the gradual rise in quasar fraction with emission-line luminosity. Willott et al. (2000) found a drop in the quasar fraction of AGN at low luminosity, postulating that their results are consistent with either the emergence of a second population of low-luminosity radio sources, which lack a well-fed nucleus, or a gradual decrease in $`\mathrm{\Theta }_{\mathrm{trans}}`$ with decreasing radio luminosity.

Traditionally, radio galaxies have been divided into two populations, FRI and FRII, (Fanaroff & Riley, 1974), based on radio structure but correlated with radio luminosity such that above $`\mathrm{log}_{10}(L_{151}/\mathrm{W}\mathrm{Hz}^1\mathrm{sr}^1)25`$, objects are FRII, whereas below that critical $`L_{151}`$, they are FRI. Chiaberge, Capetti & Celotti (2002) found that the FRII population is inhomogeneous and that not all of them can be unified with quasars. They find a low-radiative-efficiency accretion, weak or absent broad-line emission and a lack of a significant nuclear absorbing structure for weak-jet, low-ionisation narrow-line galaxies. For broad-line objects and obscured high-ionization narrow-line galaxies they see or infer intense ionizing emission, powerful jets and a torus-like absorber. In low-radio-frequency-selected samples, the fraction of objects with observed broad lines changes rapidly from $`0.4`$ above $`\mathrm{log}_{10}(L_{151}/\mathrm{W}\mathrm{Hz}^1\mathrm{sr}^1)26.5`$ to $`0.1`$ below this. The less luminous population from Willott et al. (2001) is composed of FRIs and FRIIs with weak/absent emission lines and their more luminous population of strong-emission-line FRII radio galaxies and quasars. A two-population model like this is motivated by the fact that the presence or absence of a quasar nucleus (as shown by emission-line strength) seems likely to be connected to the properties of a central engine of the radio source. Kaiser & Alexander (1997), and others before them, have argued that radio structure is influenced by the environment. It seems likely that a high jet power is necessary for a highly collimated non-dissipative jet, thus all high luminosity sources are FRII, but as lower jet powers are reached, a jet is more likely to disrupt resulting in a FRI source. The exact environmental density and homogeneity would determine the radio luminosity threshold of FRI/FRII within the lower luminosity population. Ledlow & Owen (1996) found that the FRI/FRII division is proportional to the square of the optical luminosity of the host galaxy, which is plausibly related to such environmental effects, although other causes have been proposed. This simple picture is shown to be incomplete by the discovery by Blundell & Rawlings (2001) of an optically powerful quasar with FRI radio structure. It also seems that there is evidence that some FRI sources in 3C have hidden quasar nuclei (Cao & Rawlings, 2003), and one clear example of a 3CR FRI with broad lines (3C386, Simpson et al., 1996) is known. In this paper we follow Hine & Longair (1979), Laing et al. (1994) and Jackson & Wall (1999) in assuming that the FRI/FRII division is strongly influenced by environmental effects and is much less fundamental than a division based on accretion properties of the central engines.

It is well known (e.g. Dunlop & Peacock 1990) that the evolution in comoving space density of radio galaxies depends on both $`L_{151}`$ and $`z`$. In two-population models this can be expressed as differential density evolution of these populations (Willott et al., 2001). Snellen & Best (2001) have criticized this approach on the basis of the artefacts it produces when the populations join.

The primary motivation of this work was to investigate in a quantitative manner the differences between the emission-line properties of radio-loud quasars and radio galaxies alongside the luminosity dependence of the quasar fraction. In this paper we will quantify the three-dimensional distribution in \[OIII\] emission-line luminosity $`L_{[\mathrm{OIII}]}`$, $`151`$ MHz luminosity $`L_{151}`$ and redshift $`z`$ using a new method applied to complete samples. A principal components analysis (PCA) is used to find the axis which causes the most differentiation between the objects, and it is possible that it can identify parameters which are more physically meaningful than the attributes of each data point (radio and emission-line luminosity) used in the PCA. We carry out a PCA and define a generalized luminosity function (GLF) in this new parameter space. The GLF is then used to make proper comparisons between receding-torus and two-population models, allowing for cosmic evolution of the populations.

The paper is organized as follows. In Sec. 2 the data are described. In Sec. 3, the principal components analysis is described. In Sec. 4, simple one-population GLFs are described to enable a comparison between unified schemes with and without a receding torus. The effects of using two-population GLFs are presented in Sec. 5, and in Sec. 6 we compare our results with previous studies and discuss possible meanings for the new parameters found by the PCA. We assume throughout that $`\mathrm{H}_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1,\mathrm{\Omega }_\mathrm{M}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$.

## 2 Data

### 2.1 Classification of radio sources

Most theories postulate the existence of a central accreting supermassive black hole in all radio sources, so any division between radio quasars and radio galaxies is bound to be problematic. Operationally, however, the division seems reasonably clean. Previous prescriptions for making this division (e.g. Jackson & Rawlings, 1997; Willott et al., 1998) have defined radio quasars as objects whose integrated optical luminosities are dominated by point sources (rather than the host galaxy). For most other radio sources, defined as radio galaxies, there is no direct detection of a compact nuclear source or broad emission lines, although of course both could be present but obscured. Following Hine & Longair (1979) and Laing et al. (1994), most authors tend to sub-divide radio galaxies into those with high-excitation narrow-line spectra (HEGs) and those with low-excitation narrow-line spectra (LEGs). We follow Jackson & Rawlings (1997) in defining LEGs as objects with (rest-frame) \[OIII\] equivalent widths of $`<10\mathrm{\AA }`$, \[OII\]/\[OIII\] ratios $`>1`$, or both.

Amongst HEGs, there are a number of cases where classification as a radio quasar is arguably more natural than classification as a radio galaxy. This number has tended to increase as observational methods have improved, and become more varied and sophisticated. Such objects, often called ‘weak quasars’ can be sub-divided into a number of distinct categories: (i) objects with unobscured, broad-line optical nuclei which are insufficiently luminous to outshine the host galaxy (3C382, Cohen et al., 1999; 3C390.3, Pérez et al., 1988; 5C7.17, 5C7.118, Willott et al., 1998); (ii) objects with lightly veiled (dust obscuration $`A_\mathrm{V}15`$) broad-line nuclei seen via broad wings on the H$`\alpha `$ line (3C33.1, 3C67, 3C268.3, Laing et al., 1994; 3C22, 3C41, Simpson et al., 1999; 3C109, Cohen et al., 1999; 3C303, Eracleous & Halpern, 1994) or via broad Paschen lines (3C184.1, 3C219, 3C223, Hill et al., 1996), via broad wings on rest-frame UV lines (3C325, Grimes et al. in prep), or inferred from nuclear point sources seen in the thermal IR (3C65, Simpson et al., 1999; 3C79, 3C234, Simpson, Ward & Wall, 2000) or, at high spatial resolution, with the HST (3C455, Lehnert et al., 1999); and (iii) objects with broad lines apparent only in optically-polarized light, presumably scattered from a nucleus with $`A_\mathrm{V}10`$ (3C265, Dey & Spinrad, 1996; Simpson et al., 1999; 3C226, di Serego Alighieri, Cimatti & Fosbury, 1994). The thermal-IR observations of Simpson et al. (1999) show, nevertheless, that material of high ($`A_\mathrm{V}10`$) visual extinction is needed to hide any bright nuclei in most HEGs. It is highly plausible (c.f Cyg A, Ogle et al., 1997) that some of these highly obscured nuclei will eventually show scattered (i.e. polarized) broad lines in deep spectropolarimetry.

For the LEGs there seems little doubt that any quasar nuclei are of comparatively low bolometric luminosity, but whether or not they are naked (as argued by Chiaberge, Capetti & Celotti, 1999) or buried (as predicted by the receding torus model) is still open to debate. There is one clear example of a broad-line, optically-compact nucleus in a LEG (3C386, Simpson et al., 1996; Chiaberge, Capetti & Celotti, 1999).

With this bewildering variety of ways in which the terms radio quasars and radio galaxies can become confused, it is clear that any division into two classes needs to be carefully explained. Fig. 1 is a cartoon of a classification scheme based on both the intrinsic (unobscured) bolometric luminosity of nucleus and the visual extinction to the nucleus. It is worth emphasising that how one goes about meaningfully separating radio quasars and radio galaxies depends on the nature of the question being asked. Here, as discussed in detail in Sec. 1, we are principally interested in comparing a two-population-like model in which LEGs have naked, but very low luminosity, nuclei with a receding-torus-like model in which $`L_{\mathrm{Bol}}`$ controls $`A_V`$ so that LEGs necessarily have high $`A_V`$. The natural division for this problem, illustrated in Fig. 1, separates objects with any direct evidence for a quasar nucleus, veiled or not and weak or not, from those in which there is no such evidence. In most figures we will plot weak quasars (WQs), which are significantly ($`A_V\stackrel{>}{}1`$) veiled and/or low luminosity, with a separate symbol. Whether or not one defines objects with broad lines only in scattered light is a matter of debate because broad lines could, in theory, be scattered from highly obscured nuclei. Considering such objects, we somewhat arbitrarily classify 3C226 and 3C265 as weak quasars – in the latter case because a thermal-IR nucleus with $`A_\mathrm{V}10`$ is seen (Simpson, Rawlings & Lacy, 1999) – and 3C321 as a radio galaxy because although scattered broad H$`\alpha `$ is seen (Cohen et al., 1999), transmitted broad Paschen $`\alpha `$ is not (Hill et al., 1996).

### 2.2 Complete samples

A complete sample consists of every radio source in a certain area of sky brighter than a specified flux-density limit at the specified selection frequency. Ideally, all the sources are identified optically with a radio galaxy or quasar, the redshifts and emission-line luminosities are determined for each source and the sample sky area and flux-density limits are known. Radio ($`151`$ MHz) and emission line (\[OIII\]) data were used from the 3CRR sample of Laing, Riley & Longair (1983), the 6CE sample ( Rawlings, Eales & Lacy 2001, a revision of the 6C sample of Eales 1985) and the 7CRS sample (Willott et al., 2002b).

The 3CRR sample is selected with $`S_{178}10.9`$ Jy ($`S_{151}12.4`$ Jy assuming a spectral index of $`0.8`$), and with a sky area of $`4.23`$ steradians. Three sources are excluded, (3C231 because the radio emission is due to a starburst and not an AGN, and the flat-spectrum quasars, 3C345 and 3C454.3 as they are only in the sample because of Doppler boosting), leaving $`170`$ sources (Willott et al., 2001). Of these, $`39`$ are classified as quasars, $`20`$ as weak quasars, and $`111`$ are radio galaxies, (see Table 1 for a summary of all samples). The radio galaxies consist of $`88`$ FRIIs and $`23`$ FRIs, and unless otherwise specified, the term radio galaxies will henceforth refer to both FRIIs and FRIs.

The flux-density limits of the 6CE sample are $`2.0S_{151}3.93`$ Jy, and the sky area is $`0.103`$ sr, and $`58`$ sources are used (as one object 6C1036+3616 is occluded by a bright star and is excluded without bias), of which $`49`$ are radio galaxies and $`9`$ are radio quasars.

The 7C-I and 7C-II subsamples of the 7CRS are used, consisting of $`37`$ and $`39`$ sources, with $`0.0061`$ sr and $`0.0069`$ sr and $`S_{151}0.51`$ Jy and $`S_{151}0.49`$ Jy respectively. Two sources are excluded: $`3\mathrm{C}200`$ because it is already included in the 3C sample and $`5\mathrm{C}7.230`$ as it is a flat-spectrum object which is only in the sample because of Doppler boosting. Of these objects, $`49`$ are radio galaxies, $`23`$ radio quasars and $`2`$ weak quasars (5C7.17 and 5C7.118).

Even coverage of the radio luminosity-redshift plane would provide the best constraint with which to derive luminosity functions. However from Fig. 2 we can see that the complete samples described above have a very uneven distribution on the $`L_{151}z`$ plane, because the samples are flux-density limited and quite small. However, the 3CRR, 6CE and 7CRS samples do have very high spectroscopic completeness and there are reliable redshift estimates for the small fraction of objects without spectroscopic redshifts (Willott, Rawlings & Blundell, 2001a).

### 2.3 Emission-line data

The emission-line data for the 6CE sample are given in Rawlings, Eales & Lacy (2001) and for the 7CRS sample in Willott et al. (2002b). For the 3CRR survey, the emission-line data are taken primarily from Jackson & Rawlings (1997), Hirst, Jackson & Rawlings (2003), R. Laing and J. Wall (priv. comm.) following Willott et al. (1999).

Since the sources in the samples described above have redshifts in the range $`0<z<4`$, it is not possible to measure the flux of the same emission line in each object with optical spectroscopy and the availability of near-infrared spectroscopy is limited. The \[OIII\] $`\lambda 5007`$ line is chosen as the second luminosity in the model because it is an excellent indicator of the strength of the underlying non-stellar continuum (Simpson, 1998), and it is probably produced at radii beyond those obscured by the torus. Note that Willott et al. (1999) used \[OII\] because it was the line with the most measurements.

In cases where no \[OIII\] line flux was available, other narrow-line fluxes were used to estimate it. Most commonly, the other line flux used was \[OII\]. There are $`57`$ objects in the 3CRR sample with both \[OIII\] and \[OII\] measurements. This enabled a best-fitting relation to be determined of the form $`\mathrm{log}_{10}(L_{[\mathrm{OII}]})=a+b\mathrm{log}_{10}(L_{[\mathrm{OIII}]})`$, where $`a=4.74\pm 1.44,b=0.86\pm 0.04`$ (see Fig. 3). The ‘best-fit’ line was calculated using an algorithm that minimises the sum of the squares of the perpendicular distances from the data points to a line with an adjustable slope and intercept. One reason why this relation is not a proportionality may be the systematic changes of \[OII\]/\[OIII\] ratio with narrow-line luminosity expected if the effective ionization parameter changes as a function of $`L_{\mathrm{phot}}`$ (Saunders et al., 1989): HEGs, which dominate at high $`L_{\mathrm{phot}}`$, have a systematically lower \[OII\]/\[OIII\] ratio than LEGs, which dominate at low $`L_{\mathrm{phot}}`$. The most common other lines used were Lyman-$`\alpha `$, Mg II, \[NeIV\], \[NeV\], CIII\], H$`\alpha `$ and CII\] which were used to estimate \[OIII\] using the average line ratios quoted by McCarthy (1993).

Fourteen radio galaxies in the 3CRR sample and three radio galaxies in the 6CE had no emission-line data and these objects are identified in Fig. 4. These objects are the FRI radio galaxies 3C83.1B, 3C288, 3C296, 3C310, 3C314.1, 3C315, NGC6109, NGC6251, NGC7385 and 3C465, the FRII LEGs 3C427.1, 6C1143+3703 and 6C1159+3651, and the FRII HEGs 3C68.2, 4C74.16, 3C292 and 6C1129+3710. Emission-line data for these objects was obtained by generating a random value with a Gaussian distribution function with the mean given by Eqn. 4 and with the observed scatter.

Three 3CRR quasars (3C9, 3C432 and 3C454), six 6CE quasars (6C0824+3535, 6C0913+3907, 6C1148+3842, 6C1213+3504, 6C1220+3723 and 6C1255+3700) and twelve 7CRS quasars (5C6.5, 5C6.33, 5C6.34, 5C6.95, 5C6.160, 5C6.237, 5C6.279, 5C6.287, 7C0808+2854, 5C7.70, 5C7.87 and 7C0825+2930) had $`[\mathrm{OIII}]`$ emission-line strengths estimated by assuming a rest-frame equivalent width of $`30\mathrm{\AA }`$ (Hirst et al., 2003), and the optical frequency power-law index $`\alpha _{\mathrm{opt}}=0.44`$ for $`\lambda <5000\mathrm{\AA }`$ and $`\alpha _{\mathrm{opt}}=2.45`$ for $`\lambda >5000\mathrm{\AA }`$ (Vanden Berk et al, 2001). The distribution of quasars with estimates of $`L_{[\mathrm{OIII}]}`$ from equivalent widths is found to be compatible with the distribution of quasars with real measurements in each of the three samples and overall.

The line strengths for 6C0955+3844 and 6C1045+3513 were taken from Willott et al. (2002a). A new spectrum, redshift, classification and $`L_{[\mathrm{OII}]}`$ has been found for 3C325 . An upper limit for the $`[\mathrm{OIII}]`$ emission-line strength for 3C386 was taken from Simpson et al. (1996). All objects that had only an upper limit were treated as if the upper limit was the $`[\mathrm{OIII}]`$ luminosity, but this is not expected to have a significant effect on any result.

A complete updated list of measured 3CRR line luminosities and their references are available at http://www-astro.physics.ox.ac.uk/$``$cjw/3crr/3crr.html. A table of the 3CRR, 6CE and 7CRS emission-line strengths, from either measurements or estimates, is available at http://www-astro.physics.ox.ac.uk/$``$sr/grimes.html.

### 2.4 Radio source counts

The number of known radio sources at $`151`$ MHz is very much larger than the sources in the complete samples with redshifts described above and is a valuable extra constraint. The source counts have been determined from the 6C (Hales, Baldwin & Warner, 1988) and 7C (McGilchrist et al., 1990) surveys, which have much larger sky areas than the 6CE and 7CRS redshift surveys. The 7C source counts go as faint as $`0.1`$ Jy, and so provide a low-flux-density constraint on the GLF, even in the absence of redshift information.

## 3 Principal Components Analysis

PCA is a statistical technique that has been extensively used in astrophysics for the analysis of spectral data (e.g. Folkes et al. 1999; Madgwick et al. 2002) and for reducing the dimensionality of problems. It finds how a set of properties of a sample of objects are inter-related, by identifying sets of parameters that always correlate and are most discriminatory between each source in the sample. The strength of each new parameter is measured by the amount it contributes to the variance of the sample. Sub-samples of data are identified from any clustering in the space defined by the new axes of the PCA.

There have been many treatments of PCA in the literature (e.g. Francis & Wills 1999; Murtagh & Heck 1987). PCA searches for the best-fitting set of axes to replace the initial set of axes corresponding to the attributes of the data, using the squared deviation of the points from the axes as the goodness-of-fit criterion and by enforcing orthogonality. This process gives rise to an eigenvalue equation. The matrix formed from all the eigenvectors rotates the original basis to a new set of orthogonal axes.

Practically, it is necessary to standardize the variables under analysis, so that the components of the normalized attributes $`X=\{x_{ij}\}`$ are given by

$$x_{ij}=\frac{r_{ij}\overline{r}_j}{s_j\sqrt{n}},$$
(1)
where the index $`i`$ runs from $`1`$ to $`n=302`$, the total number of objects in the 3CRR, 6CE and 7CRS samples, and $`j=1,2`$, so that $`r_1=\mathrm{log}_{10}L_{151},r_2=\mathrm{log}_{10}L_{[\mathrm{OIII}]}`$, for each attribute $`\overline{r}`$ is the mean and $`s`$ the standard deviation.

### 3.1 Results of the PCA

The PCA was performed and the average value of $`\mathrm{log}_{10}L_{151}`$ was found to be $`26.73`$ with a spread of $`1.05`$ and the average value of $`\mathrm{log}_{10}L_{[\mathrm{OIII}]}`$ was $`35.46`$ with a spread of $`1.10`$. This gave rise to the eigenvalues and eigenvectors given in Table 2. As can be seen, the first eigenvector contributes $`93`$ per cent of the total scatter. Also, the principal eigenvector $`\alpha `$ shows that the dominant relation between $`\mathrm{log}_{10}L_{151}`$ and $`\mathrm{log}_{10}L_{[\mathrm{OIII}]}`$ is that they are positively correlated, (with an equal contribution from their normalized components because there are only two components), and the secondary relation is an anti-correlation, much smaller in magnitude than the main correlation.

We have now defined two new axes on which to define a GLF. The projections of the attributes on these new axes are given by $`\alpha `$ and $`\beta `$, where

$$\alpha _i=\frac{1}{\sqrt{2}}x_{i1}+\frac{1}{\sqrt{2}}x_{i2},$$
(2)
and similarly

$$\beta _i=\frac{1}{\sqrt{2}}x_{i1}\frac{1}{\sqrt{2}}x_{i2}.$$
(3)
From Fig. 4 we can see that the first transformation of the data (top to the centre plot) involves a normalization (Eqn. 1) so that the axes are expanded. The second transformation of the data (centre to bottom plot) is a rotation (Eqns. 2 and 3) of the data. There seems to be a cutoff in $`\alpha `$-space at $`0.10`$, so that below this value there is only one (weak) quasar. Also, the radio quasars are biased to lower values of $`\beta `$, because they have higher values of $`L_{[\mathrm{OIII}]}`$ for a given value of $`L_{151}`$.

Notice that the objects from the 3CRR, 6CE and 7CRS samples can all be described by the same principal components. Since it can be seen that for any flux-density-limited sample, $`L_{151}`$ and redshift are correlated (see Fig. 2), it might be expected that different samples may have different principal components. Tests have shown that these components do not vary significantly for various combinations of samples. The most obvious conclusion is that the correlation between emission-line and radio luminosities is much more important than any $`L_{151}z`$ correlations (Willott et al., 1999). It can be shown from the normalized components that the relationship between $`L_{[\mathrm{OIII}]}`$ and $`L_{151}`$ can be expressed as

$$\left(\frac{L_{[\mathrm{OIII}]}}{\mathrm{W}}\right)=3.4\times 10^7\times \left(\frac{L_{151}}{\mathrm{W}\mathrm{Hz}^1\mathrm{sr}^1}\right)^{1.045}$$
(4)
For comparison, Serjeant et al. (1998) found that there was a radio-optical correlation in steep-spectrum quasars of the form $`L_{\mathrm{phot}}L_{408}^{0.6\pm 0.1}`$. Willott et al. (1998) found a less steep slope of $`L_{\mathrm{phot}}L_{151}^{0.4\pm 0.1}`$ for $`24`$ radio loud quasars from 7CRS, and attributed the discrepancy at least in part to the fact that Serjeant et al. (1998) did not impose a limiting absolute magnitude in their definition of a quasar. Willott (2001) found a slope of $`1.00\pm 0.04`$ in the $`\mathrm{log}_{10}L_{[\mathrm{OII}]}\mathrm{log}_{10}L_{151}`$ relation. All these previous studies adopted $`\mathrm{\Omega }_\mathrm{M}=1.0,\mathrm{\Omega }_\mathrm{\Lambda }=0.0`$. The differences between these slopes can probably largely be attributed to the subtle differences between quasars and radio galaxies.

Values of $`\alpha `$ and $`\beta `$ for all objects in the samples are available from the table at http://www-astro.physics.ox.ac.uk/$``$sr/grimes.html.

## 4 Generalized luminosity function

In this section, we will define generalized luminosity functions (GLFs) assuming that all of the radio sources are drawn from the same population, for comparison with two-population GLFs, which will be presented in Sec. 5. A GLF $`\rho (\alpha ,\beta ,z)`$ is modelled as

$$\rho (\alpha ,\beta ,z)=\rho _010^{p(\alpha \alpha _{\mathrm{break}})}g(\beta )f(z),$$
(5)
where

$$f(z)=\mathrm{exp}\left(\frac{(zz_0)^2}{2z_1^2}\right),$$
(6)
and

$$g(\beta )=\mathrm{exp}\left(\frac{(\beta \beta _o)^2}{2\sigma _\beta ^2}\right),$$
(7)
and $`\alpha `$ is the first principal component, $`\beta `$ is the second principal component, $`\sigma _\beta `$ is the scatter, $`\rho _0`$ is the normalization and $`z_1`$ and $`z_0`$ define the Gaussian evolution in redshift $`f(z)`$ which fits the RLF in Willott et al. (2001). The GLF is defined as the number of sources per unit $`\alpha `$, per unit $`\beta `$ per unit volume. Gaussian scatter in a direction perpendicular to the best-fitting slope of $`\mathrm{log}_{10}L_{[\mathrm{OII}]}\mathrm{log}_{10}L_{151}`$ is seen by Willott (2001), so Gaussian scatter in $`\beta `$ is assumed. The function is a broken power-law in $`\alpha `$, so that $`p=p_1`$ when $`\alpha >\alpha _{\mathrm{break}}`$ and $`p=0`$ otherwise. The $`\beta `$-offset parameter $`\beta _o`$ accounts for the fact that a GLF $`\rho (\alpha ,\beta ,z)`$ derived from a symmetric distribution of objects in the $`(\alpha ,\beta )`$-plane will produce a distribution that has an offset in the $`\beta `$-direction for a flux-density limited sample. This arises because, for a given redshift, we select only the highest values of $`L_{151}`$. From the definition of $`\alpha `$, this means that we see only objects with low values of $`L_{[\mathrm{OIII}]}`$ for a given value of $`\alpha `$, and thus we see objects with less negative values of $`\beta `$ than would be seen in a sample that was not flux-density-limited. An analogous problem is that the distribution of sources with respect to spectral index in a complete radio-selected sample varies with the frequency at which the sample is selected, (Kellermann, 1964; Williams & Bridle, 1967).

We consider both a simple unified scheme with constant $`\mathrm{\Theta }_{\mathrm{trans}}`$ and a receding-torus scheme with

$$\mathrm{tan}\mathrm{\Theta }_{\mathrm{trans}}=(L_{[\mathrm{OIII}]}/L_0)^{\frac{1}{2}}\mathrm{tan}\mathrm{\Theta }_0,$$
(8)
where $`L_0`$ is taken as the median value of $`L_{[\mathrm{OIII}]}`$.

A factor of $`\mathrm{sin}\theta `$ accounts for the probability of detecting an object inclined at an angle $`\theta `$ to the line of sight. The detected object is a quasar if $`\theta <\mathrm{\Theta }_{\mathrm{trans}}`$ and a radio galaxy if $`\theta \mathrm{\Theta }_{\mathrm{trans}}`$. Therefore the GLF for radio galaxies is given by

$$\rho _{\mathrm{RG}}(\alpha ,\beta ,z)=\rho (\alpha ,\beta ,z)\times \mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}},$$
(9)
and the GLF for radio quasars is given by

$$\rho _\mathrm{Q}(\alpha ,\beta ,z)=\rho (\alpha ,\beta ,z)\times (1\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}}).$$
(10)
We will denote the GLF with the receding torus as 1R and with the simple unified scheme as 1S.

A maximum-likelihood analysis was performed to optimize the free parameters of the GLF. For both GLFs, there are eight free parameters: $`\sigma _\beta ,\rho _0,p_1,\alpha _{\mathrm{break}},\beta _o,z_0,z_1`$ and $`\mathrm{\Theta }_0`$. In order to maximise the likelihood we must minimize $`S`$, where $`S=2\mathrm{ln}(\text{likelihood})`$. Following Marshall et al. (1983), we find an expression for $`S`$,

$`S`$ $`=`$ $`2{\displaystyle \underset{i=1}{\overset{\mathrm{NQ}}{}}}\mathrm{ln}[\rho (\alpha _i,\beta _i,z_i)(1\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}})]`$
$``$ $`2{\displaystyle \underset{i=1}{\overset{\mathrm{NRG}}{}}}\mathrm{ln}[\rho (\alpha _i,\beta _i,z_i)\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}}]`$
$`+`$ $`2{\displaystyle \rho _\mathrm{Q}\mathrm{\Omega }\frac{\mathrm{d}V}{\mathrm{d}z}\left|J\right|\mathrm{dzd}\mathrm{log}_{10}\mathrm{L}_{151}\mathrm{d}\mathrm{log}_{10}\mathrm{L}_{[\mathrm{OIII}]}}`$
$`+`$ $`2{\displaystyle \rho _{\mathrm{RG}}\mathrm{\Omega }\frac{\mathrm{d}V}{\mathrm{d}z}\left|J\right|\mathrm{dzd}\mathrm{log}_{10}\mathrm{L}_{151}\mathrm{d}\mathrm{log}_{10}\mathrm{L}_{[\mathrm{OIII}]}},`$
where NQ is the number of radio quasars including weak quasars, NRG is the number of radio galaxies, $`\mathrm{\Omega }=\mathrm{\Omega }(\mathrm{log}_{10}L_{151},z)`$ is the sky area available from the samples for this value of $`z`$ and $`L_{151}`$, $`J`$ is the Jacobian matrix for the change of variables from principal components to $`\mathrm{log}_{10}L_{[\mathrm{OIII}]}`$ and $`\mathrm{log}_{10}L_{151}`$, and $`\mathrm{d}z\times \mathrm{d}V/\mathrm{d}z`$ is the differential co-moving volume element. In essence, there are two terms for quasars, and two for radio galaxies: one is the sum over all NQ quasars or NRG radio galaxies in the samples, and the last two terms are the integrals over the model being tested and should give $`2`$ NQ and $`2`$ NRG respectively for good fits.

The errors associated with the parameters were found by calculating the components of the Hessian matrix ($`S`$) at the location of the minimum, inverting this matrix to obtain the covariance matrix ($`[\sigma ^2]_{ij}=2[(S)^1]_{ij}`$), and taking the $`1\sigma `$ errors as given by the square root of the diagonal elements of this matrix (e.g. Sivia 1996).

To find relative probabilities for the GLFs with or without a receding-torus, the procedure set out by Sivia (1996) was used. The ratio of the posterior probabilities of GLF 1R with respect to GLF 1S, $`P_{\mathrm{RS}}`$, is given approximately by

$`P_{\mathrm{RS}}`$ $`=`$ $`{\displaystyle \frac{P(1\mathrm{R}|\mathrm{data})}{P(1\mathrm{S}|\mathrm{data})}}`$ (12)
$`=`$ $`{\displaystyle \frac{e^{\frac{S_{1\mathrm{R}|\mathrm{min}}}{2}}\sqrt{det(S_{1\mathrm{S}}|_{S_{\mathrm{min}}})}}{e^{\frac{S_{1\mathrm{S}|\mathrm{min}}}{2}}\sqrt{det(S_{1\mathrm{R}}|_{S_{\mathrm{min}}})}}}\times F,`$
where $`det(S_{1\mathrm{R}}|_{S_{\mathrm{min}}})`$ is the determinant of the Hessian matrix for GLF 1R, evaluated at $`S=S_{\mathrm{min}}`$. All of the free parameters are common to both GLFs, except for $`\mathrm{tan}\mathrm{\Theta }_{\mathrm{trans}}`$ which is coded in a different way in each GLF (see Equation 8) but the prior ranges on this parameter are equal, giving $`F=1`$.

### 4.1 Results

The maximum-likelihood routine was run for both the GLF with a receding torus 1R and the GLF with the simple constant-$`\mathrm{\Theta }_{\mathrm{trans}}`$ unified scheme 1S. A downhill simplex (Press et al., 1992) method was used to locate the best-fitting parameters, and these are presented in Table 3. GLF 1R was found to be more likely than GLF 1S by a factor of $`10^7`$. The difference in likelihood between 1R and 1S arises mainly from the differences in their unified scheme parameters, as most of the parameters describing the shape of the GLF are very similar. This naively implies that the simple non-luminosity-dependent unified scheme is strongly ruled out, as there exists an alternative physical model which is far more probable given the data. This has put on a quantitative basis the obvious (e.g. Lawrence 1991) statement that a constant-quasar-fraction model is inconsistent with the data. In Fig. 4, it is clear that the quasar fraction changes dramatically with $`\alpha `$.

Willott et al. (2001) derived three models for the RLF at $`151`$ MHz from the low-frequency-selected 3CRR, 6CE and 7CRS samples. The GLF in $`\alpha `$ and $`\beta `$ space, $`\rho (\alpha ,\beta ,z)`$ can be easily transformed to give a RLF $`\rho (\mathrm{log}_{10}L_{151},z)`$ for comparison

$`\rho (\mathrm{log}_{10}L_{151},z)={\displaystyle \rho (\alpha ,\beta ,z)\left|J\right|\frac{\mathrm{dV}}{\mathrm{dz}}\mathrm{d}\mathrm{log}_{10}\mathrm{L}_{[\mathrm{OIII}]}}.`$ (13)
Fig. 5 shows the Model B RLF from Willott et al. (2001), which has been adapted to a $`\mathrm{\Omega }_M=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,H_0=70`$ km s<sup>-1</sup> Mpc<sup>-1</sup> cosmology, and the RLFs derived from the GLFs 1R and 1S, which have almost exactly the same shape. Note that agreement with Willott et al. (2001) is only expected to be good above points which represent the flux-density limit of 7CRS. This is because only $`L_{151}z`$ plane data has been used to constrain the GLF, but the Willott et al. (2001) RLF is also constrained by source counts.

Simulations of the $`\mathrm{log}_{10}L_{[\mathrm{OIII}]}\mathrm{log}_{10}L_{151}`$ plane from GLFs 1R and 1S are shown in Fig. 6, and show that these simple models can reproduce the data reasonably well at high values of $`L_{151}`$ and $`L_{[\mathrm{OIII}]}`$. One of the key features of the data, namely a lack of quasars at low values of $`L_{[\mathrm{OIII}]}`$, is seen in the simulations from GLF 1R, whereas the simulations from 1S clearly do not cause a difference in emission-line luminosities between radio quasars and radio galaxies. It is therefore easy to see why there is such a huge difference in likelihoods.

We can compare the receding-torus parameter with Willott et al. (2000) and Simpson (1998). Fig. 7 shows the quasar fraction as a function of $`L_{[\mathrm{OIII}]}`$ for GLFs 1R and 1S, as well as two-population GLFs which will be introduced in Sec. 5. It is clear that a constant-transition-angle unified scheme cannot reproduce the quasar fraction data from the 3CRR, 6CE and 7CRS samples across the range of $`L_{[\mathrm{OIII}]}`$, while the one-population GLF with a receding torus provides a reasonably good fit to the data. Note that the quasar fraction curve derived in Willott et al. (2000) does not include FRI objects. The same conclusions can be drawn from Fig. 8, where the quasar fraction is plotted as a function of $`\alpha `$. The quasar fraction from GLF 1R also agrees reasonably well with the receding-torus type model parameters derived from 3CRR and 7CRS $`L_{[\mathrm{OII}]}`$ data from Willott et al. (2000).

The scatter in the $`\beta `$-direction, $`\sigma _\beta `$ can be related to the scatter in the $`L_{151}L_{[\mathrm{OIII}]}`$ relation by the factor $`s_2\sqrt{2n}`$, giving a value of $`0.58`$, which is similar to the scatter of $`0.6`$ in the radio-optical correlation found by Serjeant et al. (1998), and the scatter in the $`L_{[\mathrm{OII}]}L_{151}`$ correlation of $`0.54`$ found by Willott (2001).

## 5 A two-population model

The results of Sec. 4.1 imply that some mechanism is needed to cause the change in quasar fraction with $`\alpha `$, recalling that $`\alpha `$ is a measure of the strongly correlated radio and $`[\mathrm{OIII}]`$ luminosities. An alternative mechanism to the receding-torus scheme could be a two-population model composed of a high-$`\alpha `$ population of radio galaxies and radio quasars related by a unified scheme, in addition to a low-$`\alpha `$ population mostly composed of radio galaxies, but with a unified scheme that allows a small fraction of quasars, to account for objects such as 3C386. This would certainly cause an increase in quasar fraction with $`\alpha `$.

A two-population GLF was created by adding together a low-$`\alpha `$ population GLF ($`\rho _1`$) to a high-$`\alpha `$ population GLF ($`\rho _2`$). This method is similar to that of Willott et al. (2001) for the construction of the $`151`$ MHz RLF, and utilises the same evolutionary forms of both populations. Thus,

$$\rho _1=\rho _{10}10^{p_1(\alpha \alpha _{\mathrm{cut},1})}\mathrm{exp}\left[10^{q(\alpha \alpha _{\mathrm{cut},1})}\right]g(\beta )f_1(z),$$
(14)
where $`g(\beta )`$ is defined in Eqn. (7) and

| $`f_1(z)`$ | $`=`$ | $`(1+z)^k`$ | for | $`z<z_1`$ |
| --- | --- | --- | --- | --- |
|  | $`=`$ | $`(1+z_1)^k`$ | for | $`zz_1.`$ |

Also

$$\rho _2=\rho _{20}10^{p_2(\alpha \alpha _{\mathrm{cut},2})}\mathrm{exp}\left[10^{q(\alpha \alpha _{\mathrm{cut},2})}\right]g(\beta )f_2(z),$$
(15)
where

| $`f_2(z)`$ | $`=`$ | $`\mathrm{exp}\left[\frac{1}{2}(\frac{zz_{2\mathrm{a}}}{z_{2\mathrm{b}}})^2\right]`$ | for | $`z<z_{2\mathrm{a}}`$ |
| --- | --- | --- | --- | --- |
|  | $`=`$ | $`1`$ | for | $`zz_{2\mathrm{a}}`$. |

There is an exponential fall-off in the density of population $`1`$ for $`\alpha >\alpha _{cut}`$ and in the density of population $`2`$ objects when $`\alpha <\alpha _{cut}`$. A factor $`q=s_1\sqrt{n}/\sqrt{2}`$ is the factor by which the range in $`\alpha `$ is smaller than the range in $`\mathrm{log}_{10}L_{151}`$, and acts to make the exponential decline happen as in Willott et al. (2001). Fig. 9 shows the two-population GLF.

The first GLF we consider is a two-population scheme where the high-$`\alpha `$ population $`2`$ consists of radio quasars and radio galaxies unified by a receding-torus scheme, and the low-$`\alpha `$ population $`1`$ consists of radio galaxies and quasars with a simple unified scheme. Henceforth we will refer to this GLF as 2RS. We also consider the case where the both the low- and high-$`\alpha `$ populations are unified by a simple unified scheme and we will refer to this as the 2SS GLF. In addition, we consider a GLF for a two-population scheme which has quasars and radio galaxies unified by a receding torus model in both populations, which will be referred to as Model 2RR. The radio galaxy and quasar populations densities are given by

$`\rho _{\mathrm{RG}}(\alpha ,\beta ,z)`$ $`=`$ $`\rho _1\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}_1}+\rho _2\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}_2}.`$ (16)
$`\rho _\mathrm{Q}(\alpha ,\beta ,z)`$ $`=`$ $`\rho _1(1\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}_1})+\rho _2(1\mathrm{cos}\mathrm{\Theta }_{\mathrm{trans}_2}).`$
These GLFs are constrained by the 3CRR, 6CE and 7CRS datasets as before (Sec. 4.1) and, because of the increased complexity of the two-population GLFs, we need the additional constraint of the 6C and 7C source counts (described in Sec. 2.4). In addition, the Willott et al. (2001) RLF had the constraint of the local RLF. Here we use the more recently determined local ($`z0.25`$) RLF for AGN from the 2dF Galaxy Redshift Survey (Sadler et al., 2002). We simply use their normalization of the local RLF at $`\mathrm{log}_{10}(L_{151}/\mathrm{W}\mathrm{Hz}^1\mathrm{sr}^1)=24.0`$ (converted from $`1.4`$ GHz to $`151`$ MHz using $`\alpha _{\mathrm{rad}}=0.8`$) as a prior. Without the use of this we would have very little constraint on the faint end of the GLF at low values of redshift. The slope $`p_1`$ of the low-$`\alpha `$ population can similarly be fixed at a value $`6.4`$, derived using the results of Sadler et al. (2002).

The function to be minimised is now defined as

$$S_{\mathrm{total}}=A\mathrm{ln}(\chi _\mathrm{A}^2)+B\mathrm{ln}(\chi _\mathrm{B}^2)+C\mathrm{ln}(\chi _\mathrm{C}^2),$$
(17)
where the parameters $`A,B`$ and $`C`$ are ‘Lagrange multipliers’ or Bayesian ‘hyper-parameters’ (Lahav et al., 2000), which effectively control the relative weight applied to the data sets. The function to be minimised can be shown to be

$$S_{\mathrm{total}}=N_\mathrm{A}\mathrm{ln}(S)+N_\mathrm{B}\mathrm{ln}(\chi _{\mathrm{sc}}^2)+N_\mathrm{C}\mathrm{ln}(\chi _{\mathrm{LRLF}}^2),$$
(18)
where $`S`$ is the likelihood function for two-population models, analogous to $`S`$ defined in Eqn. (4), $`\chi _{\mathrm{sc}}^2`$ is the $`\chi ^2`$ value for the radio source counts and $`\chi _{\mathrm{LRLF}}^2`$ is the $`\chi ^2`$ value for the local RLF. $`N_\mathrm{A},N_\mathrm{B}`$ and $`N_\mathrm{C}`$ are the number of data points associated with the 3CRR, 6CE and 7CRS $`L_{[\mathrm{OIII}]}`$ and $`L_{151}`$ data ($`n=302`$ objects), the binned source count data ($`30`$ 6C and $`16`$ 7C bins) and the local RLF constraint ($`1`$) respectively. The $`\chi ^2`$ is given by

$$\chi ^2=\underset{i=1}{\overset{i=\mathrm{nbins}}{}}\left(\frac{f_{\mathrm{data},i}f_{\mathrm{model},i}}{\sigma _{\mathrm{data},i}}\right)^2,$$
(19)
and for the source counts

$$f_{\mathrm{data},i}=\left(\frac{\mathrm{d}N}{\mathrm{d}N_0}\right)_i,$$
(20)
where $`\mathrm{d}N_0=2400(S_{\mathrm{min}}^{1.5}S_{\mathrm{max}}^{1.5})`$ is the Euclidean source counts.

### 5.1 Results from the two-population GLFs

The maximum-likelihood routine, where $`S`$ is defined in a similar way to that in Eqn. 4, was performed on the two-population models, to find best-fitting values of the GLF parameters, and the results are given in Table 5. GLF 2RR, with a receding-torus model in both high- and low-$`\alpha `$ populations was the best-fitting model but, in terms of probability, it is only preferred over 2RS by a factor of $`2`$ (and over 2SS by a factor $`1.3`$). The scatter about the $`\mathrm{log}_{10}L_{[\mathrm{OIII}]}\mathrm{log}_{10}L_{151}`$ relation was found to be $`0.62`$.

The form of the best-fitting GLF can be seen in Fig. 9, which shows $`\rho (\alpha ,z)`$ (the GLF integrated over $`\beta `$) against $`\alpha `$. The two-population GLF is similar to the one-population GLF at high $`\alpha `$, but they differ at low values of $`\alpha `$, where the two-population model is additionally constrained by radio source counts. The RLF derived from the two-population GLF now agrees quite well with the Model B RLF from Willott et al. (2001) over the full range of $`L_{151}`$ and $`z`$ (Fig. 10). Snellen & Best (2001) note that the two-population RLFs of Willott et al. (2001) have prominent humps due to the different evolutionary forms used for the high- and low-luminosity populations, which are identical in form to those used for the high-and low-$`\alpha `$ populations in this study. They argue that FRI and FRII radio galaxies should not be treated as intrinsically different classes of object, and that evolution should simply be a function of radio power, so that the higher power FRII objects should undergo stronger evolution than FRIs. Since FRI objects are all drawn from the low-$`\alpha `$ population of the GLFs defined in this study, their evolution is in general weaker than FRIIs, which may be drawn from either population. The RLFs derived from the GLFs are smoother than the RLFs found by Willott et al. (2001) since at each value of $`L_{151}`$ objects are drawn from a band of objects in $`\alpha `$. It seems that the RLFs derived from GLFs naturally produce smoother results because the effects of scatter, if only partially, are taken into account. The 2RS, 2SS and 2RR two-population models all give very similar GLFs, so that the small likelihood differences probably arise from the differences in their unified schemes.

The best-fitting unified scheme parameters are compared with the quasar fraction as a function of $`L_{[\mathrm{OIII}]}`$ in Fig. 7 and as a function of $`\alpha `$ in Fig. 8. The effect of combining a simple unified scheme with a two-population model can be seen, as the two-population GLF 2SS fits the data much better than the one-population GLF 1S. The error bars on the quasar-fraction data are sufficiently large that, despite the very different shape of the quasar fraction curves, the non-luminosity-dependent schemes used in 2SS do not give a significantly worse fit than GLFs employing a receding-torus model in one or both populations. These GLFs, 2RS and 2RR, fit the quasar fraction data very well, differing only at low values of $`\alpha `$, where 2RR predicts that the quasar fraction goes to zero at very low $`\alpha `$. In comparison, 2RS and 2SS seem to over-predict slightly the number of objects at low $`\alpha `$, in order to be able to account for objects like 3C386. It is obvious in this context why there is little difference between the likelihoods of the GLFs. 2RR can fit the second lowest data point in Fig. 7 and Fig. 8, whereas 2RS and 2SS can both fit the lowest data point. None of the GLFs can fit the third lowest quasar fraction data point in either Fig. 7 or Fig. 8, which seems to arise from the presence of weak quasars at approximately $`0.1<\alpha <0.3`$ (see Fig. 4).

The source counts derived from the GLFs are compared with the 6C and 7C source count data (Sec. 2.4) and the source counts derived from the RLF from Willott et al. (2001) in Fig. 11. The RLFs from Willott et al. (2001) do not give a very good fit to the source count data (reduced $`\chi ^2=3.17`$), under-producing the 6C counts at $`0.73.0`$ Jy. However, the employment of the Bayesian ‘hyper-parameters’ in the definition of the likelihood function (Eqn. 18) leads to a sufficiently strong source count constraint, so that the source counts derived from our GLFs result in an improved fit to the data (reduced $`\chi ^2=1.68`$). The difference between the data and the RLFs can be accounted for by shot noise and by deviations of the number of objects in the 7CRS sample from the average number of sources in an area of sky due to large-scale structure, (e.g. Brand et al. 2003). Another consequence of the tighter source count constraint is that, despite the good agreement of the RLFs at most redshifts (Fig. 10), the low-$`\alpha `$ population $`1`$ provides less of the source counts than the low-luminosity population of Willott et al. (2001), and correspondingly the high-$`\alpha `$ population $`2`$ provides more of the source counts than the Willott et al. (2001) high-luminosity population, (which is also evident in Fig 12).

These subtle changes can have significant impact on the redshift distributions predicted for new radio source redshift surveys. For example, the peak of the model redshift distribution is shifted to lower redshifts for the GLFs compared to the Willott et al. (2001) RLFs at the flux density limit of a redshift survey at $`0.1\mathrm{Jy}`$ (Fig 12).

The simulations of the data for the two-population GLFs are shown in Fig. 13. All three GLFs show that radio quasars are more likely to be found at high values of $`L_{[\mathrm{OIII}]}`$ compared to radio galaxies, although the distribution of quasars in the high-$`\alpha `$ population is more evenly distributed in 2SS because of the constant opening angle of the torus. The most obvious problem with the simulations is that all of the GLFs under-produce the number of objects at the lowest radio and emission line luminosities. This is a consequence of the reduced normalization of population $`1`$ compared to the Willott et al. (2001) RLF at 3CRR flux densities (see Fig. 11 and Fig 12). This could be fixed by invoking a more complicated model.

## 6 Discussion

### 6.1 Quantitative Comparison of Unified Schemes

A preliminary investigation (Sec. 4) showed that a receding-torus scheme is strongly preferred over a simple non-luminosity-dependent unified scheme. However when a two-population model was investigated, there was little evidence to favour a receding torus in the high-$`\alpha `$ population. It is clear that, with enough scatter in the radio - optical relation, a small quasar fraction in the low-$`\alpha `$ population combined with even a constant value of $`\theta _{\mathrm{trans}}`$ in the high-$`\alpha `$ population will mimic the rise of quasar fraction with emission-line luminosity predicted by a luminosity-dependent unified scheme such as the receding-torus model. A simple unified scheme alone will not generate the observed increase in quasar fraction with $`\alpha `$ and the emission-line differences between radio quasars and radio galaxies but a receding-torus model or a two-population scheme with any sort of unified scheme in the high-$`\alpha `$ population will produce these effects. On the basis of quasar fractions and emission-line differences alone, there is not enough evidence to distinguish between possible luminosity dependences of the unified schemes, apart from ruling out schemes with no luminosity dependence, such as model 1S.

Simpson (2003) reviews the evidence for the receding-torus model highlighting the quasar fraction and emission-line difference arguments. However we have shown here quantitatively that the quasar fraction and the differences in emission-line luminosity in radio-selected samples are not very significant either individually or when considered together.

A major cause of the inability of the data to discriminate between the competing models is that the quasar fractions of the 6CE sample is very different ($`15.5\%`$) to the 3CRR ($`34.7\%`$) and 7CRS samples ($`33.8\%`$). Fig. 14 shows the quasar fractions as a function of the limiting flux density for the 3CRR, 6CE, 7CRS and 7CQ (Riley et al., 1999) samples and the curves predicted from the 2RS, 2SS and 2RR models. This illustrates the fact that although the 2SS model leads to very different predictions than the receding-torus model, they cannot yet be clearly distinguished by the data.

To re-investigate the issue of the emission-line differences between radio quasars and radio galaxies, we repeat the experiment of Jackson & Browne (1990). They compared each of $`12`$ 3C radio galaxies with radio quasars with $`0.19<z<0.85`$, matched in redshift to within $`0.05`$ and in luminosity to within $`30`$ per cent, and find that quasar $`[\mathrm{OIII}]`$ luminosities exceed those of radio galaxies by a factor of $`510`$. We plot histograms of the number of 3CRR radio quasars and radio galaxies with $`0.2z0.8`$ in $`L_{[\mathrm{OIII}]}`$ bins, Fig. 15. There is clearly an offset in the distributions, with the median value of $`L_{[\mathrm{OIII}]}`$ being $`36.10`$ for quasars and $`35.50`$ for radio galaxies, giving a shift of $`0.60`$, so that the median quasar has $`4`$ times more luminous \[OIII\] emission than the median radio galaxy. We found the distributions in $`L_{[\mathrm{OIII}]}`$ of 3CRR objects in this redshift range for the three two-population models. All three models reproduce a shift in the median value of $`L_{[\mathrm{OIII}]}`$ between quasars and radio galaxies. The 2RS model gave a shift of $`0.78`$ compared to $`0.56`$ for 2SS and $`0.64`$ for 2RR. This analysis shows that any of the two-population models considered here can reproduce the results of Jackson & Browne (1990) so that their test does not provide clear support for the receding-torus model. All models, with or without a receding torus, predict that quasars are brighter than radio galaxies because the probability of an object being classified as a quasar increases with $`\alpha `$ and $`L_{[\mathrm{OIII}]}`$.

### 6.2 The meaning of $`\alpha `$ and $`\beta `$

It is instructive to consider, if only briefly, possible interpretations for the axes found by the principal components analysis in the context of the receding-torus scheme and the two-population scheme. Since $`\alpha `$ accounts for the overwhelming majority of the scatter, it must be a manifestation of some important physical process. Higher values of $`\alpha `$ mean higher low-frequency radio luminosity and higher emission-line luminosity, i.e. more powerful objects. Emission-line luminosity is often assumed to be dominated by photo-ionization from a strong nucleus, though there is evidence that some part, specifically low-ionization lines like $`[\mathrm{OII}]`$, is powered by radiative shocks (Inskip et al., 2002). In this context an increase in accretion rate will cause an increase in both radio luminosity and emission-line luminosity (Rawlings & Saunders, 1991), although it is not immediately obvious in what proportions. We tentatively suggest that $`\alpha `$ is an indicator of accretion rate and $`\beta `$ is the scatter about this relation, which can probably be attributed to a number of factors: for example, environment and black hole mass. In this context and for the two-population model, it seems possible that $`\alpha _{\mathrm{cut},1}`$ and $`\alpha _{\mathrm{cut},2}`$ represent respectively a critical accretion rate above which AGN are shining at a significant fraction of their Eddington luminosity, and a critical rate below which they are shining at a tiny fraction.

In essence, it can be shown that, for the very similar values of scatter in the radio and emission-line distributions as found for this sample,

$$\alpha \mathrm{log}_{10}\left(\frac{L_{151}L_{[\mathrm{OIII}]}}{L_{151,\mathrm{av}}L_{[\mathrm{OIII}],\mathrm{av}}}\right),$$
(21)
and

$$\beta \mathrm{log}_{10}\left(\frac{L_{151}/L_{[\mathrm{OIII}]}}{L_{151,\mathrm{av}}/L_{[\mathrm{OIII}],\mathrm{av}}}\right).$$
(22)
$`L_{151}`$ gives an indication of the jet power which is in turn powered by the central engine and $`L_{[\mathrm{OIII}]}`$, arising from photo-ionization, is also a measure of the power of the central engine. $`\alpha `$ represents the product of the radio (mechanical) and optical (radiative) output, compared to the average product. $`\beta `$ represents the ratio of radio output to optical output compared to the average ratio.

As $`\beta `$ probably represents the scatter from a number of different sources, it would be interesting to see if $`\beta `$ correlates with any of the most likely contributions towards scatter. It has been shown by Magorrian et al. (1998) that the black hole mass correlates with the mass of the bulge in galaxies, and assuming a constant mass-to-light ratio, we expect the stellar luminosity to correlate with the black hole mass. Willott et al. (2003) find that the best-fit $`Kz`$ relation from the 3CRR, 6CE, 7CRS and 6C samples is given by

$$K=17.37+4.53\mathrm{log}_{10}z0.31(\mathrm{log}_{10}z)^2.$$
(23)
In Fig. 16 a measure of the stellar luminosity $`\mathrm{\Delta }K`$ (the excess $`K`$ luminosity over the $`Kz`$ relation of Eqn. 23 <sup>1</sup><sup>1</sup>1Willott et al. (2003) show that the empirical $`Kz`$ relation of Eqn. 23 lies very close to the locus expected for a giant elliptical of fixed mass whose stellar population formed at high redshift and is evolving passively. $`\mathrm{\Delta }K`$ represents (in magnitudes) the excess luminosity.) is plotted against $`\beta `$ for $`160`$ radio galaxies from the 3CRR, 6CE and 7CRS (I and II) samples (Willott et al., 2003). The Spearman rank correlation coefficient is $`0.064`$ with significance $`0.42`$ implying that $`\beta `$ is not significantly correlated with black hole mass. The correlation coefficient between $`\alpha `$ and K is $`0.131`$ with significance $`0.10`$.

Lacy et al. (2001) show that, for a sample of quasars from the FIRST Bright Quasar Survey (FBQS) and optically-selected quasars from the Palomar-Green survey, $`L_{5\mathrm{G}\mathrm{H}\mathrm{z}}M_{\mathrm{BH}}^{1.9\pm 0.2}(L/L_{\mathrm{Edd}})^{1.0}`$, where $`M_{\mathrm{BH}}`$ is the black hole mass. This agrees with the slope of $`1.05`$ found here in the $`\mathrm{log}_{10}L_{151}\mathrm{log}_{10}L_{[\mathrm{OIII}]}`$ relation if the narrow emission-line luminosity, $`L_{[\mathrm{OIII}]}`$ is proportional to the bolometric luminosity, $`L`$, and if the range in black-hole mass is small. We conclude that the scatter $`\sigma _\beta `$ is not dominated by the black-hole mass.

There are of course a variety of other factors which might correlate with $`\beta `$. For example, implicit in the work of Rawlings & Saunders (1991) is that $`\beta `$ might be expected to correlate positively with the richness of the radio source environment as, for a given jet power (assumed to scale with bolometric quasar luminosity), simple physical models predict that radio sources have higher radio luminosities in dense-gas environments (see also Barthel & Arnaud 1996), leading to more positive values of $`\beta `$. Sadly, despite the existence of XMM-Newton and Chandra, there are as yet no comprehensive studies of radio source environments that can be used to make such a test. We label on Fig. 16, several radio sources known to be in rich galaxy clusters from ROSAT observations (Hardcastle et al., 1988). The radio galaxy 3C295 ($`\beta =0.054`$) is one of the highest-$`\beta `$ radio galaxies in the sample, in agreement with the prediction, but 3C220.1 (which is not in Fig. 16 as it does not have a measured K magnitude) has $`\beta =0.015`$ and 3C280 has $`\beta =0.013`$, compared to the average value of $`\beta =0.0025`$ for radio galaxies in the redshift range $`0.4<z<1`$. Hardcastle et al. (1988) also finds that five 3CRR quasars with $`0.35<z<0.75`$ reside in rich cluster environments, and they have an average $`\beta =0.0024`$, compared to an average of $`0.0037`$ for quasars in this redshift range. So there are hints that this correlation may be present, but obviously the numbers involved are very small. We have also failed to find any significant correlation between $`\beta `$ and projected radio source linear size $`D`$. Naively one might expect a negative correlation because synchrotron and adiabatic losses should cause a systematic decline in radio luminosity over the lifetime of the radio source if the jet power and narrow-line luminosity remain fixed (Kaiser, Dennett-Thorpe & Alexander, 1997; Blundell, Rawlings & Willott, 1999). However, as discussed by Willott et al. (1999), such a correlation may be masked by a number of effects including changes with time of jet power and $`L_{\mathrm{phot}}`$, and the boosting of emission-line fluxes in small sources due to radiative bow shocks (Inskip et al., 2002).

## 7 Conclusions

A new approach to investigating unified schemes has been presented, based on a principal components analysis of the 3CRR, 6CE and 7CRS complete samples. Generalized luminosity functions have been derived based on these principal components: $`\alpha `$ which encodes the $`L_{151}L_{[\mathrm{OIII}]}`$ correlation and $`\beta `$ which encodes the scatter about this correlation. The main advantage of this new approach has been that the unified scheme parameters have been found by taking into account the intrinsic scatter in the $`L_{151}L_{[\mathrm{OIII}]}`$ correlation. The main conclusions to be drawn from this analysis are as follows.

1. A receding-torus model is strongly favoured over a simple non-luminosity-dependent unified scheme for GLFs with one population of radio sources.
2. With the extra constraint of 6C and 7C radio source counts and the normalization of the local RLF, two-population GLFs were derived. The GLFs give rise to RLFs which are very similar but smoother than those of Willott et al. (2001). There is very little difference in likelihoods between a GLF with a receding torus in both populations, a GLF with a torus opening angle that does not vary with ionizing luminosity in both populations, and a GLF with a receding-torus in the high-$`\alpha `$ population and a constant-opening-angle torus in the low-$`\alpha `$ population.
3. Two-population models reproduce the radio survey data well, and can provide a natural explanation for the rise in quasar fraction with emission-line luminosity and the emission-line differences between radio quasars and radio galaxies. The receding-torus may be a feature in both populations but this is not yet proved.

## Acknowledgements

JAG acknowledges the support of a graduate studentship from Oxford University and Magdalen College. SR acknowledges the support of a PPARC Senior Fellowship. This research has made use of the NASA/IPAC Extragalactic Database, which is operated by the Jet Propulsion Laboratory, Caltech, under contract with the National Aeronautics and Space Administration. We would like to thank Chris Simpson for a critical reading of an early version of the manuscript, and an anonymous referee for some very useful comments.


format: text

text:
# Hierarchical Clustering Using Mutual Information

## Abstract

We present a method for hierarchical clustering of data called mutual information clustering (MIC) algorithm. It uses mutual information (MI) as a similarity measure and exploits its grouping property: The MI between three objects $`X,Y,`$ and $`Z`$ is equal to the sum of the MI between $`X`$ and $`Y`$, plus the MI between $`Z`$ and the combined object $`(XY)`$. We use this both in the Shannon (probabilistic) version of information theory and in the Kolmogorov (algorithmic) version. We apply our method to the construction of phylogenetic trees from mitochondrial DNA sequences and to the output of independent components analysis (ICA) as illustrated with the ECG of a pregnant woman.

Classification or organizing of data is very important in all scientific disciplines and is fundamental for understanding and learning jain-dubes . Classification can be exclusive or overlapping, supervised or unsupervised. In the following we will be interested only in exclusive unsupervised classification, called clustering.

An instance of a clustering problem consist of a set of objects and a set of properties (called characteristic vector) for each object. The goal of clustering is separation of objects into groups using only the characteristic vectors. Cluster analysis organizes data either as a single grouping of individuals into non-overlapping clusters or as a hierarchy of nested partitions. The latter is called hierarchical clustering (HC). Because of wide spread of applications, there are a large variety of different clustering methods in usage, see e.g. jain-dubes  for an overview.

The crucial point of all clustering algorithms is the choice of a proximity measure. This is obtained from the characteristic vectors and can be either an indicator for similarity or dissimilarity. In the latter case it is convenient but not obligatory to satisfy the standard axioms of a metric (positivity, symmetry, and triangle inequality). Among HC methods one should distinguish between those where one uses the characteristic vectors only at the first level of the hierarchy and derives the proximities between clusters from the proximities of their constituents, and methods where the proximities are calculated each time from their characteristic vectors. The latter strategy (which is used also in the present paper) allows of course for more flexibility but might also be computationally more costly.

Quite generally, the “objects” to be clustered can be either single (finite) patterns (e.g. DNA sequences) or random variables, i.e. probability distributions. In the latter case the data are usually supplied in form of a statistical sample, and one of the simplest and most widely used similarity measures is the linear (Pearson) correlation coefficient. But this is not sensitive to nonlinear dependencies which do not manifest themselves in the covariance and can thus miss important features. This is in contrast to mutual information (MI) which is also singled out by its information theoretic background cover-thomas . Indeed, MI is zero only if the two random variables are strictly independent.

Another important feature of MI is that it has also an “algorithmic” cousin, defined within algorithmic (Kolmogorov) information theory li-vi  which measures the similarity between individual objects. For a thorough discussion of distance measures based on algorithmic MI and for their application to clustering, see li1 ; li2 .

Essential for the present application is the grouping property of MI,

$$I(X,Y,Z)=I(X,Y)+I((X,Y),Z).$$
(1)
Within Shannon information theory this is an exact theorem, while it is true in the algorithmic version up to the usual logarithmic correction terms li-vi . Since $`X,Y,`$ and $`Z`$ can be themselves composite, Eq.(1) can be used recursively for a cluster decomposition of MI. This motivates the main idea of our clustering method: instead of using e.g. centers of masses in order to treat clusters like individual objects in an approximative way only, we treat them exactly like individual objects when using MI as proximity measure and We thus propose the following scheme for clustering $`n`$ objects with MIC:
(1) Compute a proximity matrix based on pairwise mutual informations; assign $`n`$ clusters such that each cluster contains exactly one object;
(2) find the two closest clusters $`i`$ and $`j`$;
(3) create a new cluster $`(ij)`$ by combining $`i`$ and $`j`$;
(4) delete the lines/columns with indices $`i`$ and $`j`$ from the proximity matrix, and add one line/column containing the proximities between cluster $`(ij)`$ and all other clusters;
(5) if the number of clusters is still $`>2`$, goto (2); else join the two clusters and stop.

Shannon Theory: Here, $`XX_1,YX_2,\mathrm{}`$ are random variables. If they are discrete, entropies are defined as usual $`H(X)=_ip_i(X)\mathrm{log}p_i(X)`$ etc. The MI is defined as

$$I(X_1,\mathrm{},X_n)=\underset{k=1}{\overset{n}{}}H(X_k)H(X_1,\mathrm{},X_n).$$
(2)
Eq.(1) can be checked easily, together with its generalization to arbitrary groupings. It means that MI can be decomposed into hierarchical levels. By iterating it, one can decompose $`I(X_1\mathrm{}X_n)`$ for any $`n>2`$ and for any partitioning of the set $`(X_1\mathrm{}X_n)`$ into the MIs between elements within one cluster and MIs between clusters.

For continuous variables with densities $`\mu _X`$ etc., one first introduces some binning (‘coarse-graining’), and applies the above to the binned variables. If $`x`$ is a vector with dimension $`m`$ and each bin has Lebesgue measure $`\mathrm{\Delta }`$, then $`p_i(X)\mu _X(x)\mathrm{\Delta }^m`$ with $`x`$ chosen suitably in bin $`i`$, and

$$H_{\mathrm{bin}}(X)\stackrel{~}{H}(X)m\mathrm{log}\mathrm{\Delta }$$
(3)
where the differential entropy is given by

$$\stackrel{~}{H}(X)=𝑑x\mu _X(x)\mathrm{log}\mu _X(x).$$
(4)
Notice that $`H_{\mathrm{bin}}(X)`$ is a true (average) information and is thus non-negative, but $`\stackrel{~}{H}(X)`$ is not an information, can be negative, and is not invariant under homeomorphisms $`x\varphi (x)`$.

Joint entropies, conditional entropies, and MI are defined as above, with sums replaced by integrals. Like $`\stackrel{~}{H}(X)`$, joint and conditional entropies are neither positive (semi-)definite nor invariant. But MI, defined as

$$I(X,Y)=𝑑x𝑑y\mu _{XY}(x,y)\mathrm{log}\frac{\mu _{XY}(x,y)}{\mu _X(x)\mu _Y(y)},$$
(5)
is non-negative and invariant under $`x\varphi (x)`$ and $`y\psi (y)`$. It is (the limit of) a true information,

$$I(X,Y)=\underset{\mathrm{\Delta }0}{lim}[H_{\mathrm{bin}}(X)+H_{\mathrm{bin}}(Y)H_{\mathrm{bin}}(X,Y)].$$
(6)
In applications, one usually has the data available in form of $`N`$ sample points $`(x_i,y_i),i=1,\mathrm{}N`$ which are assumed to be i.i.d. realizations. There exist numerous algorithms to estimate $`I(X,Y)`$ and entropies. We use in the following the MI estimators proposed recently in Ref. mi , and we refer to this paper for a review of alternative methods.

Algorithmic Information Theory: In contrast to Shannon theory where the basic objects are random variables and entropies are average informations, algorithmic information theory deals with individual symbol strings and with the actual information needed to specify them. To “specify” a sequence $`X`$ means here to give the necessary input to a universal computer $`U`$, such that $`U`$ prints $`X`$ on its output and stops. The analogon to entropy, called here usually the complexity $`K(X)`$ of $`X`$, is the minimal length of an input which leads to the output $`X`$, for fixed $`U`$. It depends on $`U`$, but it can be shown that this dependence is weak and can be neglected in the limit when $`K(X)`$ is large li-vi .

Let us denote the concatenation of two strings $`X`$ and $`Y`$ as $`XY`$. Its complexity is $`K(XY)`$. It is intuitively clear that $`K(XY)`$ should be larger than $`K(X)`$ but cannot be larger than the sum $`K(X)+K(Y)`$. Finally, one expects that $`K(X|Y)`$, defined as the minimal length of a program printing $`X`$ when $`Y`$ is furnished as auxiliary input, is related to $`K(XY)K(Y)`$. Indeed, one can show li-vi  (again within correction terms which become irrelevant asymptotically) that

$$0K(X|Y)K(XY)K(Y)K(X).$$
(7)
Notice the close similarity with Shannon entropy. The algorithmic information in $`Y`$ about $`X`$ is finally

$$I_{\mathrm{alg}}(X,Y)=K(X)K(X|Y)K(X)+K(Y)K(XY),$$
(8)
and similarly for more than two strings. Within the same additive correction terms, one shows that it is symmetric, $`I_{\mathrm{alg}}(X,Y)=I_{\mathrm{alg}}(Y,X)`$, and can thus serve as an analogon to mutual information.

$`K(X)`$ is in general not computable. But one can easily give upper bounds: The length of any input which produces $`X`$ (e.g. by spelling it out verbatim) is an upper bound. Improved upper bounds are provided by any file compression algorithm.

MI-Based Distance Measures: When comparing objects with different marginal or joint informations, it seems intuitively clear that one should prefer relative distances over absolute ones, in order to minimize the dependence on the total information. We here use the quantity li1 ; cluster

$$D(X,Y)=1\frac{I(X,Y)}{H(X,Y)}$$
(9)
which is a metric , with $`D(X,X)=0`$ and $`D(X,Y)1`$ for all pairs $`(X,Y)`$. The algorithmic version is also universal: If $`XY`$ according to any non-trivial distance measure, then $`XY`$ also according to $`D`$.

A difficulty appears in the Shannon framework, if we deal with continuous random variables. As we mentioned above, $`\stackrel{~}{H}(X,Y)`$ is not invariant under homeomorphisms (including rescalings) and not even positive definite, while $`H_{\mathrm{bin}}`$ diverges when $`\mathrm{\Delta }0`$. We thus modified Eq.(9) by replacing $`H(X,Y)`$ by $`H_{\mathrm{bin}}(X,Y)`$ and replacing $`D(X,Y)`$ by the similarity measure

$$S(X,Y)=\underset{\mathrm{\Delta }0}{lim}(D(X,Y)1)\mathrm{log}\mathrm{\Delta }=\frac{I(X,Y)}{m_x+m_y}.$$
(10)
A Phylogenetic Tree for Mammals: We study the mitochondrial DNA of a group of 34 mammals (see Fig. 1). The same data Genebank  had previously been analyzed in li1 ; Reyes00 . This group includes among others some rodents, ferungulates, and primates.

Obviously we are here dealing with the algorithmic version of information theory, and informations are estimated by lossless data compression. For constructing the proximity matrix between individual taxa, we proceed essentially as in Ref. li1 , using the special compression program GenCompress GenComp .

In Ref. li1 , this proximity matrix was then used as the input to a standard HC algorithm (neighbour-joining and hypercleaning) to produce an evolutionary tree. Instead we use the MIC algorithm with distance $`D(X,Y)`$. The joining of two clusters is obtained by simply concatenating the DNA sequences. There is of course an arbitrariness in the order of concatenation sequences: $`XY`$ and $`YX`$ give in general compressed sequences of different lengths. But we found this to have negligible effect on the evolutionary tree.

The overall structure of this tree closely resembles the one shown in Ref. Reyes00 . All primates are correctly clustered and also the relative order of the ferungulates (blue whale to horse) is in accordance with Ref. Reyes00 . On the other hand, there are a number of connections which obviously do not reflect the true evolutionary tree, see for example the guinea pig with bat and elephant with platypus. But the latter two, inspite of being joined together, have a very large distance from each other, thus their clustering just reflects the fact that neither the platypus nor the elephant have other close relatives in the sample. All in all, however, already the results shown in Fig. 1 capture surprisingly well the overall structure shown in Ref. Reyes00 . Dividing MI by the total information is essential for this success. If we had used the non-normalized $`I_{\mathrm{alg}}(X,Y)`$ itself, the clustering algorithm used in li1  would not change much, since all 34 DNA sequences have roughly the same length. But our MIC algorithm would be completely screwed up: After the first cluster formation, we have DNA sequences of very different lengths, and longer sequences tend also to have larger MI, even if they are not closely related.

The concatenation of $`X`$ and $`Y`$ will of course not lead to a plausible sequence of the common ancestor, but it optimally represents the information about it. This information is essential to find the correct way through higher hierarchy levels of the evolutionary tree, and it is preserved in concatenating.

Clustering of Minimally Dependent Components in an Electrocardiogram: As our second application we choose a case where Shannon theory is the proper setting. We show in Fig. 2 an ECG recorded from the abdomen and thorax of a pregnant woman ECGdata . It is already seen from this graph that there are at least two important components in this ECG: the heartbeat of the mother and of the fetus. In addition there is noise from various sources (muscle activity, measurement noise, etc.). While it is easy to detect anomalies in the mother’s ECG from such a recording, it would be difficult to detect them in the fetal ECG.

As a first approximation we can assume that the total ECG is a linear superposition of several independent sources (mother, child, noise1, noise2,…). A standard method to disentangle such superpositions is independent component analysis (ICA) ICA . There, one tries to recover the sources by means of linear transformation $`s_i(t)=_{j=1}^nW_{ij}x_j(t)`$, where $`W_{ij}`$ is determined by minimizing the estimated MI between the $`s_i`$.

In reality things are not so simple. For instance, the sources might not be independent, the number of sources (including noise sources!) might be different from the number of channels, and the mixing might involve delays. For the present case this implies that the heartbeat of the mother is seen in several reconstructed components $`s_i`$, and that the “independent” components are not independent at all. In particular, all components $`s_i`$ which have large contributions from the mother form a cluster with large intra-cluster MIs and small inter-cluster MIs. The same is true for the fetal ECG, albeit less pronounced. To obtain clean recordings of the fetal and maternal ECGs, we proceeded as follows Harald .

Since we expect different delays in the different channels, we first used Takens delay embedding Takens80  with time delay 0.002 s and embedding dimension 3, resulting in $`24`$ channels. We then formed 24 linear combinations $`s_i(t)`$. We use the MI estimator mi , for details see cluster . Five of the resulting least dependent components contain strong contributions of the mother’s heartbeat, three are dominated by the fetus. The rest contains mostly noise cluster .

In plotting the actual dendrogram (Fig. 3) we used $`S(X,Y)`$ for the cluster analysis but used the MI of the clusters to determine the height at which the two branches join. The MI of the first five channels, e.g., is $`1.44`$ nats, while that of channels 6 to 8 is $`0.3`$ nats. For any two clusters (tuples) $`X=X_1\mathrm{}X_n`$ and $`Y=Y_1\mathrm{}Y_m`$ one has $`I(X,Y)I(X)+I(Y)`$. This guarantees, if the MI is estimated correctly, that the tree is drawn properly. The two slight glitches (when clusters (1 - 14) and (15 - 18) join, and when (21 - 22) is joined with 23) result from small errors in estimating MI. They do in no way effect our conclusions.

In Fig. 3 one can clearly see two big clusters corresponding to the mother and to the child. There are also some small clusters which should be considered as noise. For reconstructing the mother and child contributions to Fig. 2, we have to decide on one specific clustering from the entire hierarchy. We decided to make the cut at inter-cluster MI equal to 0.1, i.e. two clusters $`X`$ and $`Y`$ are joined whenever $`I((X),(Y))I(X,Y)I(X)I(Y)0.1`$. Reconstructing the first five traces of the original ECG from the child components only, we obtain Fig. 4.

In summary, we have shown that MI can not only be used as a proximity measure in clustering, but that it also suggests a conceptually very simple and natural hierarchical clustering algorithm. We do not claim that this algorithm, called mutual information clustering (MIC), is always superior to other algorithms. Indeed, MI is in general not easy to estimate. Obviously, when only crude estimates are possible, also MIC will not give very good results. But as MI estimates are becoming better, also the results of MIC should improve. The present paper was partly triggered by the development of a new class of MI estimators for continuous random variables which have very small bias and also rather small variances mi .

We have illustrated our method with two applications, one from genetics and one from cardiology. For neither application MIC might give optimal clustering, but it seems promising that one common method gives decent results for both, although they are very different.

The results of MIC should improve, if more data become available. This is trivial, if we mean by that longer time sequences in the application to ECG, and longer parts of the genome. It is less trivial that we expect MIC to make fewer mistakes in a phylogenetic tree, when more species are included. The reason is that close-by species will be correctly joined anyhow, and families – which now are represented only by single species and thus are poorly characterized – will be much better described by the concatenated genomes if more species are included.

We would like to thank Arndt von Haesseler, Walter Nadler and Volker Roth for many useful discussions.