© ESO, 2014

1. Introduction

Since hydrogen is the most common element in the universe, and makes up most of its mass, the distribution of this element is intimately connected with the dynamics of matter on scales from the cosmic down to the stellar. Although in its ionized or molecular forms hydrogen is problematic to detect, neutral atomic hydrogen (H i) emits (or absorbs, if it is cooler than the background), due to a hyperfine transition, radio waves at a rest frequency of 1420 MHz. The transition rate is very low, which makes the line faint and relatively difficult to detect, but this has two compensating advantages: firstly, clouds of hydrogen up to the galactic scale are usually optically thin, which makes it relatively straightforward to deduce the mass of the cloud from the intensity of the line; and secondly that the spectral line is intrinsically very narrow, making it a useful way to study the motions, whether thermal, internal or bulk, of the emitting cloud, via Doppler shift and broadening.

In recent years there have been a number of blind surveys for H i in the local universe, some complete (e.g. HIPASS, Barnes et al. 2001), others ongoing (ALFALFA, Giovanelli et al. 2005; EBHIS, Winkel et al. 2010; Kerp et al. 2011; CHILES, Fernández et al. 2013). Further, deeper blind searches for H i (LADUMA, Holwerda et al. 2012; WALLABY, Koribalski & Stavely-Smith 2009; DINGO, Meyer 2009; see also Duffy et al. 2012) are planned with the upcoming SKA precursor instruments MeerKAT (Jonas 2009) and ASKAP (Johnston et al. 2008). Such surveys complement those performed at other wavelengths and avoid potential biases from preferential selection of galaxies which are bright at these wavelengths.

The majority of hydrogen in the universe occurs in galaxies. The width of a galaxy’s H i spectral profile gives its bulk rotation speed, and the area under the profile is proportional to its H i mass – both quantities of importance in cosmology.

From simple geometrical considerations it is clear that in any survey of objects uniformly distributed in space, the number frequency of objects will increase as their angular size decreases. Since the detectability of objects takes a sharp downturn as their angular sizes become smaller than the instrument resolution, the most common survey object therefore can be expected to be one which is only just brighter than the survey sensitivity cut-off, and which has an angular size no larger than the beam of the instrument. Since detection of H i sources usually makes use also of spectral information, a galaxy at the limit of detection may be unresolved in any one channel map, but nevertheless kinematically resolved, such that the flux peak moves progressively across channels (see the EBHIS observation of DDO 154 described in Sect. 3.3.3 for an example of this).

The principal task therefore for post-detection processing of blind H i surveys is to measure as accurately as possible the width, area and other properties of low signal-to-noise (S/N) spectral lines of unresolved sources: and not only for spatially-integrated spectral lines, but also for the line at each spatial pixel across the extent of the source. This is the aim of the technique discussed in the present paper.

Fitting of a model to a spatially integrated spectral profile where the source is well resolved on the other hand should not be expected to yield better values of galaxy parameters than other methods, because with such sources there are problems determining which spatial pixels to include in the sum – significant numbers of pixels with small but non-zero contributions may be excluded, or pixels containing nothing but a chance spike in noise may be mistakenly included. We make use of the well-resolved THINGS observations here for practical, demonstrative reasons. It is unavoidable that pixel masking effects cause flux biases in such spectra, although careful treatment of the data can minimize this.

The question of spatial masking or selection is one that arises even if the source is unresolved, because its flux remains distributed across the imaging beam or point spread function (PSF) of the instrument. Any cut-off criterion will miss some contribution in the wings, but a too generous cut-off will include too much noise from off-source directions. Fitting a model which includes a spatial component shaped like the PSF seems like a simple way around this, but is unsatisfactory because the noise in adjacent spatial pixels is also convolved by the PSF, and is therefore not statistically independent. Statistically speaking, it is better to fit to data which is as unprocessed as practical – i.e., to include instrumental response in the model and fit to raw data, rather than to fit a simpler model to data which has been processed (with accompanying muddying of the statistical waters) in order to remove or at least systematize the instrumental response. But an exploration of such matters is beyond the scope of the present paper. Here we content ourselves with construction of global spectra by spatial masking and summing, and make the associated caveats about the reliability of resulting flux measurements.

The plan of the paper is as follows. The advantages of having a model of the H i spectral line is discussed in Sect. 2.1, and the model itself is described in Sect. 2.2. Various technical matters connected with the process of fitting, including the Bayesian methodology, estimation of uncertainties, and a suitable goodness-of-fit measure, are discussed in Sect. 2.4. The tests performed on the model are presented in Sect. 3. H i spectral lines were obtained under a variety of conditions of noise, spectral resolution and baseline, and were fitted by the model.

We posed three questions about the fits. Firstly, how good a fit is the model in the case of optimum frequency resolution and S/N – does it return good values of bulk properties of the galaxies? This question is addressed in Sect. 3.1, where we fit the model to a spatially integrated spectrum from each of the 34 THINGS galaxies observed by Walter et al. (2008).

Secondly, if we add noise to the spectra, and bin the channels more coarsely, does this affect the fitted parameter values? We address this question in Sect. 3.2.

Thirdly, how does the model perform fitting to a single-dish observation, where it is necessary also to fit a baseline? And can the model return useful kinetic information when the galaxy is at the limit of resolution? In Sect. 3.3, data from the EBHIS survey (Winkel et al. 2010; Kerp et al. 2011) are used to test this.

2. A model of the H I spectral line

2.1. Introduction

The aim of any model of a physical process is to approximate that process with a small, hence manageable, number of adjustable parameters. A good model will be a good match to the data; will have a small number of parameters; and it is also useful if the model is not just empirical but tied in some way to the underlying physics of the object. Since these desirables can be in opposition, a model is usually then also a compromise.

The shape of an H i spectral line represents the distribution and motions of neutral hydrogen within a galaxy. On top of the bulk rotation which prevents gravitational collapse there can be innumerable variations in the pattern of local H i flow from galaxy to galaxy, which will be reflected in differences between the shapes of the respective spectral lines. It might therefore seem difficult to formulate a simple model of the H i spectral line. However as shown in Sect. 3.1.2, in practice the 6-parameter model presented here works fairly well. Local deviations in flux density between the data and the fitted model don’t exceed about 10%, and tend not to affect either the total flux or the linewidth.

In all the H i surveys the authors are aware of, total flux and linewidth have been estimated from the data directly, without fitting a model to the line profile¹. So why use one? Firstly because it is more systematic: there is no need for either human intervention or ad-hoc prescriptions for the number of channels to consider. Secondly, as is shown in Sect. 3.2, the bias in linewidth measurement is much reduced through model fitting. Thirdly, as shown in Sect. 3.3, the modelling approach very naturally accommodates a modelling of the spectral background or baseline. It’s no longer necessary to decide where the line profile “ends” before fitting the baseline, since both can be considered together. Fourthly, fitting of a model arguably lends itself more easily to automated processing, which becomes an ever more pressing consideration as survey datasets grow in size.

Lastly, a parametrized model opens the door to the use of Bayesian techniques, which are becoming increasingly accepted as useful tools in astronomy. Advantages here fall under three main heads. The first of these is that a Bayesian formulation is the formally correct (and therefore optimum) procedure for estimating the parameters of interest and for incorporating prior knowledge. This method is applied in the present paper. The second head or category is the use of Bayes’ theorem to assess the relative suitability of differing models. This is also used in the present paper, to determine the best order of baseline model (Sect. B.6). Note that the same formulation can be used to estimate detection probability directly, although this topic is not explored in depth in the present paper: the Bayesian approach automatically takes into account both prior knowledge of the expected range of line shapes, as well as the total bandwidth and area of the survey. Such an approach to detection is more fundamental and rigorous than relying on either the 5-sigma rule or ad hoc prescriptions for calculating S/N. The third advantage of the Bayesian formalism is in the extraction of statistical results from large ensembles of low-quality data via hierarchical modelling (Loredo 2012). This technique is however not used here.

Roberts (1978) discussed the H i spectral line and explained why it often has a double-horned shape with a relatively steep rise and fall. Roberts explored several semi-physical models of the line profile, and in fact his model C is a special case of ours. Roberts was concerned rather to emphasize and explore the connection between galactic linewidth and luminosity now more usually associated with Tully & Fisher (1977), and his profile models don’t have the flexibility to fit a wide variety of H i line shapes.

In the second of two papers which presented a comprehensive simulation of gas in the local universe, Obreschkow et al. (2009b) described a profile model (which they apply to both H i and CO spectral lines) which uses 5 parameters. These parameters, labelled by the authors k₁ to k₅, are purely empirical in themselves, but can be derived from more fundamental properties of the line profile via a set of formulae (Eqs. (A.2) to (A.6), Appendix A of their paper). This second tier of parameters includes the flux density at profile centre, the maximum flux density, and the velocity widths at the 50% and 20% levels. In an earlier paper by some of the same authors (Obreschkow et al. 2009a), the second-tier parameters are derived by constructing each profile from an appropriate projection and integration of a comprehensive model of the distribution of H i gas density and velocity as a function of radius from the galaxy centre. This radial model in turn depends on several third-tier fundamental parameters, such as masses and characteristic radii for disk, halo and bulge components. It might be possible to find a reasonably simple way to connect these fundamental parameters with the eventual k values, but these authors did not attempt this themselves. Their intent was to generate the second-tier parameters for a large number of simulated galaxies and make these available in a database. The k-prescription was provided simply as a convenience for the user who wished to construct approximate profiles from these data.

One could make use of the k-model of Obreschkow et al. for fitting to real data, but there are a number of desirable improvements:

• 1.

  The k parameters provide no direct physical insight.

• 2.

  Many real H i spectral lines exhibit a noticeable asymmetry (see e.g. Richter & Sancisi 1994). The k model does not cater for this.

• 3.

  Since the prescribed functions of the k model are discontinuous at sharply defined velocity bounds, not aligned a priori with velocity channel boundaries, integrating the model over a finite channel width is a little fiddly.

• 4.

  There is no way to simulate any artefacts arising from the autocorrelation process which generates spectra from radio signals; nor can the k-model simulate applied filters, such as Hanning or Tukey filters.

The model described in the present paper addresses all these issues.

2.2. The model

In order to derive a model we need to arrive at a reasonable approximation to the motion of H i in a generic galaxy. We will begin by breaking the motions into two categories: bulk vs. random.

2.2.1. Bulk motions

We approximate bulk motions firstly by assuming that all gas in a (late-type, therefore gas-rich) galaxy rotates about a common centre and in a common plane. As shown for example in de Blok et al. (2008), this is a reasonable assumption for the majority of spiral galaxies. We don’t assume that the density of gas is azimuthally symmetric: departures from same are catered for, albeit in a crude manner, via the asymmetry parameter of our model.

As is well known, for most galaxies with a total H i mass larger than about 10⁹ solar masses, the curve of rotation speed as a function of radius from the centre is seen to be remarkably flat outside the core. We assume here that the rotation curve is always exactly flat outside a certain radius. Within the core itself, the velocity usually rises steadily from (nominally) zero at the galactic centre, then turns over smoothly when it reaches the “flat” value of velocity. In the present model we approximate this rise by a straight line, such as would be observed for example in a rotating solid body; we also assume that the density of gas in the core is uniform. This simplified rotation curve is similar in shape to that labelled “C” in Fig. 3 of Roberts (1978).

Understanding of the way a profile model represents the bulk gas motions is assisted by mapping the gas motions and densities in phase space – that is, presenting the gas density not as a function of spatial location but of velocity components. The line profile may then be obtained by projecting the phase space gas distribution onto a line directed towards the observer. This is diagrammed in Fig. 1. In phase space, the outer mass of gas with a constant rotation speed appears as a ring of infinitesimal thickness, whereas the “solid-rotating” inner part appears as a uniform disk.

The ring gives in projection the characteristic double-horned profile so often observed in H i spectral lines: ${\mathit{s}}_{\mathrm{outer}}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{S}}{\mathit{\pi }\mathrm{\Delta }\mathit{v}}{\mathit{\rho }}^{-1}\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}}{\mathrm{\Delta }\mathit{v}\mathit{/}\mathrm{2}}\right)$(1)where $\begin{array}{ccc}{\mathit{\rho }}^{\mathrm{2}}\mathrm{\left(}\mathit{u}\mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{-}{\mathit{u}}^{\mathrm{2}}\mathrm{for}\mathrm{\left|}\mathit{u}\mathrm{\right|}\mathit{<}\mathrm{1}\mathit{,}\\ & \mathrm{=}& \mathrm{0}\mathrm{else}\mathrm{.}\end{array}$Here S is the total flux from the H i in this section, Δv is the range between maximum and minimum gas velocities in the line of sight, and v_ctr is the mean line-of-sight velocity. The disk in projection yields a half ellipse:

${\mathit{s}}_{\mathrm{inner}}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\frac{\mathrm{4}\mathit{S}}{\mathit{\pi }\mathrm{\Delta }\mathit{v}}\hspace{0.17em}\mathit{\rho }\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}}{\mathrm{\Delta }\mathit{v}\mathit{/}\mathrm{2}}\right)\mathrm{·}$These are the two fundamental components of our model. Asymmetry in the line profile is accommodated by multiplying both components by

$\mathrm{1}\mathrm{+}\mathrm{2}\mathit{\alpha }\mathrm{(}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}\mathrm{)}\mathit{/}\mathrm{\Delta }\mathit{v}\mathit{.}$To describe the model so far we require 5 parameters: the line centre v_ctr; the so-called intrinsic line width Δv; the total flux S; the fraction f of the gas which is found in the “solid-rotating” part which is associated with the inner part of the galaxy; and the asymmetry parameter α. The model so far is thus represented by $\begin{array}{ccc}{\mathit{s}}_{\mathrm{intrinsic}}\mathrm{\left(}\mathit{v}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}\mathit{S}}{\mathit{\pi }\mathrm{\Delta }\mathit{v}}\hspace{0.17em}\left[\mathrm{1}\mathrm{+}\frac{\mathrm{2}\mathit{\alpha }\mathrm{(}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}\mathrm{)}}{\mathrm{\Delta }\mathit{v}}\right]\\ & & \mathrm{\times }\hspace{0.17em}\left[\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}{\mathit{\rho }}^{-1}\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}}{\mathrm{\Delta }\mathit{v}\mathit{/}\mathrm{2}}\right)\mathrm{+}\mathrm{2}\mathit{f\rho }\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}}{\mathrm{\Delta }\mathit{v}\mathit{/}\mathrm{2}}\right)\right]\mathrm{·}\end{array}$(2)

2.2.2. Random motions

Random motions include both thermal motions of the atoms and turbulent motions. The latter may occur on length scales much larger than the atomic, but provided they are much smaller than the scale of bulk, i.e. galaxy-scale motions, we can lump them together with thermal motions. In the present study we assume that the distribution of random velocities is the same throughout any galaxy. Maps of the local linewidth in Walter et al. (2008) indicate that this is a reasonable assumption for the THINGS galaxies. In this approximation, random motions may be treated mathematically as a 3-dimensional convolution of the bulk-motion line profile. Further, if we approximate the distribution of random velocities by a 3D Gaussian, then the 1D projection of this in any direction is also a Gaussian: $\mathit{g}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }{\mathit{v}}_{\mathrm{rand}}\sqrt{\mathrm{2}\mathit{\pi }}}\exp \left(\frac{\mathrm{-}{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}\mathrm{[}\mathrm{\Delta }{\mathit{v}}_{\mathrm{rand}}{\mathrm{]}}^{\mathrm{2}}}\right)\mathrm{·}$(3)The characteristic width Δv_rand is the sixth and final model parameter.

It should be emphasized that representation of random gas motions by a single Gaussian is only an approximation. In reality there may be several co-located components exhibiting different degrees of dispersion (see for example Braun 1997).

The full equation for the profile model is $\mathit{s}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}{\mathit{s}}_{\mathrm{intrinsic}}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathit{\star }\mathit{g}\mathrm{\left(}\mathit{v}\mathrm{\right)}$(4)where ⋆ indicates convolution. The full list of 6 model parameters is v_ctr, Δv, S, Δv_rand, f and α. These have natural ranges, that is, ranges imposed by physical reasonableness, as follows:

• v_ctr: constrained in practice by the ends of the spectrum.

• Δv> = 0

• S> 0

• Δv_rand> 0

• 0 < = f< = 1

• − 1 < = α< = 1.

Note however that, for some of the fits reported in the present paper, up to 6 additional parameters were used for fitting the baseline.

The asymmetry and fraction-solid parameters can be expected to be correlated respectively with the third and fourth moments of the spectral line, as described by Andersen & Bershady (2009). These authors make a comprehensive study of the statistics of H i and H ii spectral line shapes and their relation to galaxy morphology. We don’t here retrace that analysis, but only observe that calculation of the appropriate moments of our spectral line model is easy and avoids the necessity, mentioned by Andersen & Bershady, of a priori human selection of a velocity range, when extracting moments from data. Other measures of asymmetry, such as that of Tifft & Cocke (1988), are equally easy to calculate from our model.

2.3. Fourier-space formulation

Although the model can be calculated directly using Eq. (4) and the preceding formalism, there are several advantages to calculating the profile first in Fourier space, then transforming to velocity space. Firstly, the convolution in Eq. (4) turns into a product; secondly, the singularities in Eq. (1) are avoided; and thirdly, the process mimics the processing of real signals in an XF-type correlator, and thus allows inclusion of some of the artefacts which result from same.

This Fourier-space model construction was followed in all the calculations described in the present paper. To make it easy for readers to do this themselves, the Fourier transform of Eq. (4) is given in Appendix A.

2.4. Fitting considerations

A frequent use for such a model is to fit it to data. For this purpose one needs an objective function describing the goodness of the fit, and one must choose an algorithm for minimizing this function. These considerations are discussed in Sects. 2.4.1 through 2.4.4.

2.4.1. Bayesian formulation

According to Bayes’ theorem, given a set of measurements of flux density y, the posterior probability density function of the model parameters p(q | y), where we use q as shorthand for the six model parameters, is given by (see e.g. D’Agostini 2003) $\mathit{p}\mathrm{\left(}{q}\mathrm{\right|}{y}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathit{E}}\hspace{0.17em}\mathit{p}\mathrm{(}{q}\mathrm{)}\hspace{0.17em}\mathit{p}\mathrm{(}{y}\mathrm{\left|}{q}\mathrm{\right)}\mathit{,}$(5)where E, known as the evidence, is just a normalizing constant: $\mathit{E}\mathrm{=}\mathrm{\int }\mathit{d}{q}\hspace{0.17em}\mathit{p}\mathrm{\left(}{q}\mathrm{\right)}\hspace{0.17em}\mathit{p}\mathrm{\left(}{y}\mathrm{\right|}{q}\mathrm{)}\mathit{.}$(6)The first function in the integrand is the prior probability distribution which represents our prior knowledge of the parameter values; the second is the likelihood. For N data values y_j which include Gaussian-distributed noise, the likelihood is given by $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{y}\mathrm{\right|}{q}\mathrm{)}& \mathrm{=}& \underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{j}}\sqrt{\mathrm{2}\mathit{\pi }}}\exp \begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{-}\mathrm{[}{\mathit{y}}_{\mathit{j}}\mathrm{-}\mathit{s}\mathrm{(}{\mathit{v}}_{\mathit{j}}\mathit{,}{q}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{j}}^{\mathrm{2}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\\ & \mathrm{=}& \mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{-}\mathit{N}\mathit{/}\mathrm{2}}\exp \left(\frac{\mathrm{-}{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{j}}}\end{array}$(7)where s(v_j,q) is the profile model evaluated for velocity channel j, σ_j is the standard deviation of the noise in channel j, and χ² has its usual formulation.

A Bayesian fit (that is, a fit procedure which optimizes the Bayesian posterior probability) may tend towards being data-dominated, or it may be prior-dominated. The first case occurs if the data has high S/N and if there has not been much earlier fitting experience. The original THINGS dataset observed with the VLA matches this criterion for us. The EBHIS observations of the THINGS galaxies however generally have much lower S/N values, as do the semi-simulated profiles described in Sect. 3.2. For the latter two cases we therefore thought it appropriate to set some cautious priors, derived from the fits to the VLA THINGS profiles. Some of the model parameters are poorly constrained by the data in these lower-S/N fits and these are thus prior-dominated.

The low-S/N priors are described in detail in Appendix B.2.

2.4.2. Fitting algorithms

Three fitting techniques have been used in the present study: the Levenberg-Marquardt (LM) method (Press et al. 1992, chapter 15.5), simplex optimization (Nelder & Mead 1965, see also Press et al. 1992 Chap. 10.4) and Markov chain Monte Carlo (MCMC), specifically using the Metropolis algorithm (Metropolis et al. 1953, see also Bhanot 1988 for other useful references). For fitting to the THINGS profiles, both from the original VLA observations (Sect. 3.1) and as observed as part of EBHIS (Sect. 3.3), the LM procedure was used to obtain initial values of the approximate centre and width of the posterior, then the exact form of the posterior was explored via MCMC. For the semi-artificial profiles described in Sect. 3.2, the simplex procedure alone was used.

The simplex algorithm is robust but slow. In practice it was “fast enough”, requiring only about 10 to 20 s on a standard laptop to fit a 6-parameter model to 450 data points.

Technical details of the MCMC are given in Appendix B.4.

2.4.3. Uncertainties

The parameter uncertainties quoted in Table 1 are the square roots of the diagonal elements of a covariance matrix estimated from the ensemble of MCMC points. For the LM and simplex fits, the covariance matrix was obtained where necessary by inverting the matrix of second derivatives (known as the Hessian) of the posterior with respect to the parameters. As can be seen in the figures in Appendix C, the MCMC values are often about 50% larger than the LM ones. The reason for this is not known.

All flavours of uncertainty scale with the uncertainties in the original values of flux density. The calculation of these is described in Appendix B.5.

2.4.4. Goodness of fit

The value of χ² is commonly used to assess how well a model fits the data. In effect, reduced χ² is a measure of the ratio between the data-minus-model residuals and the amplitude of the measurement noise. Noise plays the vital role in this because it affects the probability that the observed residuals would occur by chance with a perfect fit. In the present case however this is not quite what we want. We expect from the start that the model will not be a perfect fit to the data, so in a sense this question is already decided in the negative: what we want is rather to measure the size of the imperfections. Noise plays no role in this, provided only that it is not so large as to swamp the deviations; hence the bare value of χ² will not return the information we want.

The line shape of NGC 925 in Fig. 2 provides a good example of the issue. Clearly seen are many local “bumps” or “wiggles” in the H i distribution which the model is not able to track. It is the relative size of these bumps and dips compared to the average height of the profile which it would be most useful to know. For example, from Fig. 2 one can clearly see that the model is a cleaner fit to NGC 3184 than to NGC 925: what we want is some formula to quantify this difference.

The formula we devised to estimate the wiggle fraction J is ${\mathit{J}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{S}}^{\mathrm{3}}}{\left(\frac{\mathrm{\Delta }\mathit{v}}{\mathrm{\Delta }{\mathit{v}}_{\mathrm{chan}}}\right)}^{\mathrm{2}}\sum _{\mathit{j}}{\mathit{y}}_{\mathit{j}}\left[{\mathrm{(}{\mathit{y}}_{\mathit{j}}\mathrm{-}{\mathit{s}}_{\mathit{j}}\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathit{j}}^{\mathrm{2}}\right]\mathit{.}$(8)where S and Δv are respectively the total flux and line width model parameters as described in Sect. 2.2, and Δv_chan is the width of the spectral channels. The reasoning behind this formula is as follows. Firstly, we want the square of the residual in each channel. We then want to subtract the noise from this in quadrature. The resulting term should be weighted by the flux density y_j of the line profile. This is normalized by dividing by the sum S of the model profile values. The result so far is the square of the average residual, within the line profile, due solely to model/data mismatch. Finally the square root is taken and that result divided by the average flux density of the model, which is

$\mathrm{⟨}\mathit{s}\mathrm{⟩}\mathrm{=}\frac{\mathrm{\Delta }{\mathit{v}}_{\mathrm{chan}}}{\mathrm{\Delta }\mathit{v}}\mathit{S}\mathit{.}$The result is the average fractional residual due to “wiggles”. This value is given in Col. 8 of Table 1.

Local fluctuations in H i can be visually deceptive. A good example is the global spectrum of NGC 4214, which is shown together with a model fit in Fig. 2. At first inspection it is not clear why the model has not better fitted the seemingly regular double-horned profile. However, the linewidth of this galaxy is, at 56 km s^-1, relatively small: less than 5 times the “sigma” width of the turbulent broadening, which here is 12 km s^-1. Further consultation of Table 1 shows that the fraction of solid rotation fitted (which tends to reduce the depth of the valley between the horns) is small – consistent with zero. In fact the fitting procedure has chosen the sharpest and deepest possible double-horned profile which, after smoothing by the turbulent-broadening Gaussian, is consistent with the line slopes. With this amount of smoothing it is not possible for the model to follow deviations over velocity scales less than 10 km s^-1, as occurs with NGC 4214. These fluctuations only fool the eye into assuming a regularity which isn’t in fact there, because NGC 4214 is too narrow to make obvious the actually random nature of its wiggles; in contrast to the spectrum of NGC 925, for example (shown in the same figure).

It is worth observing too that, although the wiggles in NGC 4214 do rather offend the eye, the J factor for this galaxy is only 2.3%, well below the average seen in Fig. 3.

3. Tests of the model

3.1. Spectra from the original THINGS observations

3.1.1. Introduction

Walter et al. (2008) used the VLA to observe the H i distribution in 34 nearby late-type galaxies. Spatial resolution, frequency resolution and S/N were all relatively high, certainly when compared to the most common objects detected in blind H i surveys. The set of observations is known as THINGS (The H i Nearby Galaxy Survey). We fitted our profile model to a spectrum extracted from THINGS observations of each of the 34 targets.

Walter et al. (2008) provide these spectra in the online version of their paper. According to their description, the global spectral profile for each galaxy is obtained by adding together a subset of pixels for each channel of the data cube. Pixels where there was no measurable emission from H i were not included in the sum. This has the effect of making the noise in a channel proportional to the square root of the number of unmasked pixels. The mask cubes are not available on the THINGS web site but were kindly provided on request by F. Walter. This allowed us to estimate the noise per channel for each global profile. The exact procedure for doing so is described in Appendix B.5.

Some of the THINGS galaxies have a velocity range which straddles zero. For some of these, interference from Milky Way H i is evident. A few channels where this effect was obvious have been omitted when fitting to the spectra for galaxies DDO 53, NGC 2976, NGC 3077 and NGC 6946. Note that the total flux values in Col. 4 of Table 1 are the values of the respective parameter of the fitted model: thus these represent an interpolation over missing channels for the four galaxies mentioned. One can however also obtain the total flux under the fitted profile by adding together the samples of the profile model obtained at each channel. For Fig. 4, in which the flux under the data spectrum is compared to the model, the relevant channels have been omitted from these sums for both the data and the model for the four galaxies mentioned. The percentage differences in model flux from the values in the Table are respectively 7, 13, 4 and 3.

The question discussed in the remainder of Sect. 3.1 is: how good a fit is the model in the case of optimum frequency resolution and S/N – does it return good values of bulk properties of the galaxies?

3.1.2. Graphical and tabular display of fit results

Five examples of profiles fitted to THINGS spectra are shown in Fig. 2. These are the five galaxies chosen as “simulation inputs” in Sect. 3.2. Figures showing all 34 individual fit results are given in Appendix C.

The results of the fits are given in Table 1. Shown in Cols. 2 to 7 are the mean values with uncertainties for each of the six model parameters. Also given in Col. 8 is the “wiggle fraction” J calculated from Eq. (8). A low value indicates a good match between model and data. Note that it is possible for J² to be negative. In this case, formally speaking, the root is imaginary, as listed in the table. This has no physical meaning, it just indicates that any local deviations in H i from the model are insignificant compared to the measurement noise.

A histogram showing the distribution of wiggle fraction is given in Fig. 3.

Column 10 of the table gives the velocity separation between points on the fitted profile where the flux density decreases to 20% of its maximum value. The remaining Cols. 9 and 11, are described in the following section.

3.1.3. Numerical comparisons between model and data

The reason for doing these fits is to see how well the model matches a variety of real spectra. One can gain a qualitative impression from looking at the profiles but it is more informative to perform some numerical comparisons between the fitted model parameters and equivalent values from other sources. Some examples are shown in Figs. 4 to 7, which are described individually below.

Figure 4 compares the total flux S_fit from summing valid values of the fitted flux density to a similar sum S_data over the raw data values. What is displayed in the figure is the fractional difference between these “fit” and “data” flux values for each galaxy.

For about half of the galaxies, the difference between fluxes S_data and S_fit as displayed in Fig. 4 is consistent with zero. Almost all of the galaxies have flux differences less than about 1%. Where there is a perceptible difference, the total flux values for the data (S_data) are always larger than for the fitted profile.

No detailed explanation for this flux anomaly is known at present, but any explanation ought to start with the observation that the model is nonlinear in some of its parameters, and that it doesn’t have support across the whole spectrum. Worth particular attention are the line wings, which in the model become closer to Gaussian in shape the further away from line centre one goes. This Gaussian component is intended to model the distribution of random gas velocities. However, it is known that the true distribution of velocities is better described by a sum of at least 2 Gaussians (Ianjamasimanana et al. 2012). It is not hard to see that a Gaussian which is the best fit to the steep line edges might nevertheless miss flux present in broad, non-Gaussian wings for example.

As a toy model, with a non-linear parameter and limited support, to demonstrate how a pedestal can lead to flux underestimation in such circumstances, consider a simplified situation in which we wish to fit the two-step profile shown in Fig. 5 with a simple top-hat model with three free parameters: its left and right edges, and its height s.

Obviously there is nothing to be gained by moving the right edge of the model away from v₂; and given any two values for the edge locations, determining the best-fit height is trivial. The only open question is where within the range [ v₀,v₁ ] should we place the left edge v of the profile such that χ² is minimized. We define a kind of non-discrete, noiseless analog of χ² as

${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{v}}_{\mathrm{1}}\mathrm{-}\mathit{v}\mathrm{)}\mathrm{(}\mathit{s}\mathrm{-}{\mathit{s}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{1}}\mathrm{\left)}\mathrm{\right(}{\mathit{s}}_{\mathrm{1}}\mathrm{-}\mathit{s}{\mathrm{)}}^{\mathrm{2}}\mathit{.}$After optimizing the model height s this becomes

${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{s}}_{\mathrm{1}}\mathrm{-}{\mathit{s}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{1}}\mathrm{)}\frac{{\mathit{v}}_{\mathrm{1}}\mathrm{-}\mathit{v}}{{\mathit{v}}_{\mathrm{2}}\mathrm{-}\mathit{v}}\mathrm{·}$Clearly this has its lowest allowed value when v = v₁. But this means that the profile is not fitting the pedestal at all, and the data then has more flux than the model, which is consistent with the characteristic upward deviations seen in Fig. 4. Pedestals are noticeable in 5 of the galaxy profiles shown in Figs. C.1 to C.34, namely for Ho II, NGC 628, NGC 5194, NGC 5236 and NGC 5457. Of these, all but Ho II are among the group of 6 worst cases of flux underestimation.

Figure 6 compares the velocity widths at the 20% height of the fitted profiles (W₂₀) to the same figure calculated from the raw data values. For both model and data, the 20% flux density level was calculated with respect to the maximum value for that spectrum; the velocity at this level on either wing of the profile was calculated by linear interpolation between the most distal pair of velocities which straddled the 20% level. This simple technique is only accurate for the data because of the relatively high S/N of the THINGS observations. Some discussion of difficulties which arise in linewidth estimation when the spectrum is noisy is given in Sect. 3.2.2.

The agreement between the fitted width and the widths from the data is seen to be very good. In few cases is the difference larger than the 1% level. The anomalously high value (W_20,fit>W_20,data) belongs to the M81 dwarf B; the three low values (W_20,data>W_20,fit) to DDO 53, NGC 4449 and NGC 5236. For M81 dwB and DDO 53, the width of the spectrum channels is several percent of the W₂₀ width, and therefore entirely accounts for the error. In NGC 4449 and 5236 the profile model is clearly seen not to be a good fit. However, we note that NGC 4449 is an interacting galaxy with a significant amount of H i outside the disk. And for NGC 5236 the THINGS observations have many missing spacings, which may explain the curious steps in its line wings, which are probably the cause of the inflated value of W_20,data.

De Blok et al. (2008) fitted a tilted-ring model to a 19 member subset of the THINGS galaxies and obtained high-precision H i rotation curves for these. In Fig. 7 we compare their systemic velocities with two measures of line centre from the fitted profile, namely the fitted line centre parameter v_ctr, and the velocity ⟨ V₂₀ ⟩ obtained from the mean of the 20% velocities. The systemic velocities V_sys are taken from Col. 6 of Table 2 of de Blok et al. With exception of the outlier at about −14 km s^-1, which is NGC 3627, the spread in differences is centred on zero with a standard deviation on the same order as the channel width for these galaxies, which was about 5.2 km s^-1 for 12 out of the 19, half that for the rest.

De Blok et al. found that a V_sys calculated from the global line wings of NGC 3627 returned values in the range 717 to 720 km s^-1, which fall much closer to our value of 720.6. NGC 3627 appears to have kinematic asymmetries, such as have been described and discussed by Swaters et al. (1999). However, in the two galaxies studied by Swaters et al. the disturbances seem to be confined to the inner velocity regions. Disturbances in the inner-disk kinematics don’t affect the rise and fall of the global profile, and thus should have only a small effect on the V_sys derived via fitting the present model. In contrast to this, for NGC 3627, Fig. 79 of de Blok et al. shows a significant excess of gas at velocities 10 km s^-1 or more greater than their tilted-ring fit, right at the trailing edge of the velocity distribution. The present model is sensitive to such distortions, so this is arguably the reason for the inconsistency between our value of V_sys for this galaxy and that of de Blok et al. (2008).

Andersen & Bershady (2009) found that line centres obtained via a moment analysis were biased in proportion to the asymmetry as measured by the 3rd moment. Since asymmetry is built in to our model from the start, we would not expect it to suffer from a similar problem. Figure 8 confirms that this is the case. This figure is analogous to Fig. 8 of Andersen & Bershady. Despite that both our velocity scatter and uncertainties are smaller, no correlation is apparent in our Fig. 8.

3.2. Fitting to coarsened and noisified data

3.2.1. Introduction

We argue in the introduction that the main use for this profile model is to extract information in a systematic way from the kinds of H i spectra most commonly encountered in blind surveys, namely with low S/N, and a channel width of about 10 km s^-1. In order to test how well the model serves this purpose, data is needed with a priori known values of line width, total flux etc. to allow comparison with the values obtained from fitting the model. Also a large number of test spectra is desirable to reduce or at least ascertain the uncertainties in the results. The raw THINGS spectra will not serve for this: they are too few in number, have unknown a priori values, have too high S/N, and significantly higher spectral resolution than we normally expect to encounter in blind H i surveys.

An alternative would be simply to simulate the required numbers of noisy, coarsely binned profiles. However such a simulation always carries with it some uncertainty about how applicable its results are to real measurements. Simulations are also vulnerable to the criticism that one only gets out only what one puts in.

We can however make use of the THINGS profiles to generate realistic data simply by binning them into wider velocity channels and adding a lot more noise. Lewis (1983) made use of a similar method in investigating line width measurement biases. This is also the approach which is taken in the present section.

What we do here is run several Monte Carlos simulations, in each of which an ensemble of test spectra is generated. The line-profile model is fitted to each spectrum in an ensemble, and the width, centre and total flux of the line are also estimated using traditional direct methods. For each ensemble, and for each property of interest (line width etc.), an average is formed from the values obtained from the fitted profiles on the one hand and the direct measurement on the other. The input values of the properties are known and so are available for comparison. The aim is to demonstrate that the average value of each property is better estimated via model fitting than by direct measurements from the spectrum.

3.2.2. The simulation Monte Carlo

Five of the THINGS galaxies were selected, the sample being chosen so as to cover a range of shapes and sizes. These five are the ones shown in Fig. 2. They all had channel widths of about 2.6 km s^-1, except for NGC 4214, for which the width is half this. For each of the five, a Monte Carlo ensemble of 100 spectra was generated for each of a set of twelve S/N values, the S/N figure being calculated according to the ALFALFA formula, which is discussed in Sect. 3.2.3. The twelve values were evenly spaced in a logarithmic sense over a range between about 2 and 100. Different amounts of S/N offset were applied to the five galaxies to prevent overlay of graph points when results from several galaxies appear on the same figure.

As already mentioned, the desired channel width in this Monte Carlo was about 10 km s^-1. This was most simply achieved by averaging an integer number of the channels of the original spectra, that number being 4 for all but NGC 4214, for which a factor of 8 was used. This scheme also permitted dithering of individual spectra by random offsets in the range from zero to one less than the rebinning factor. This dithering was to avoid systematic effects from the centre velocity of the input spectrum always having the same relation to the channel boundaries in the rebinned scheme.

3.2.3. Calculating the S/N

Signal-to-noise ratio can be defined most simply as the ratio between the maximum or peak height of the input spectrum and the noise standard deviation or root mean square (rms). However this does not correlate well with the detectability of H i lines, because a broader line is usually easier to detect than a narrow one having the same peak-to-rms S/N. For this reason it is usual to define some alternative measure which takes account of the line width. The ALFALFA survey makes use of the following formula (Eq. (2) of Haynes et al. 2011):

$\mathit{S}\mathit{/}{\mathit{N}}_{\mathrm{ALFALFA}}\mathrm{=}\frac{\mathit{S}}{\mathit{\sigma }{\mathit{W}}_{\mathrm{50}}}\sqrt{\frac{\min \mathrm{\left(}\mathrm{400}\mathit{,}{\mathit{W}}_{\mathrm{50}}\mathrm{\right)}}{\mathrm{2}\mathrm{\Delta }{\mathit{v}}_{\mathrm{chan}}}}\mathrm{·}$We used the same to calculate the S/N values for our Monte Carlo.

Note firstly that, for any given profile shape, S/N_peak ∝ S/N_ALFALFA, but the proportionality constant differs from profile to profile. Note secondly that ALFALFA take S/N_ALFALFA> 6 as their detection criterion (Haynes et al. 2011).

A sample of a spectrum produced in the Monte Carlo is shown in Fig. 9. This has a S/N equal to 3.9, which is a little below the ALFALFA detection criterion. As such it represents a very typical example of the quality of profile commonly seen in such surveys.

3.2.4. Fitting to the data and analysing the results

The simplex algorithm was used to fit each spectrum. Several fitting runs were done with increasing values of maximum permitted number of iterations, and decreasing values of the convergence criterion, in order to be sure that the fits were properly converging. The posterior described in Sect. 2.4.1 was the object function fitted, but this time some attention was paid to priors. Details of these and their justification can be found in Appendix B.2.

There are two ways one can analyse such data: one can compare the fitted values of the model parameters to those values obtained from fits to the original THINGS template spectra; and one can compare measurements derived from fitting a model to those obtained via more traditional, non-parametric methods. Both sorts of analysis have been done for the “coarse” data. The corresponding results are described respectively in Sects. 3.2.5 and 3.2.6.

It is always going to be difficult to obtain reliable estimates of line parameters from data at such low S/N values as in Fig. 9. There are three quantities which are of central importance: the total flux under the spectral line, its central or systemic velocity, and its width.

In the case of total flux, a typical non-parametric estimation method is that described in Haynes & Giovanelli (1984): i.e. simply to “integrate over the observed signal”. Implicit in this however is a knowledge of where the “observed signal” begins and ends. For high S/N signals this is often done by visual inspection. In this case exactitude is not too important – one can afford to be generous, since the added noise per extra channel is small. Judging the line boundaries is of course much more problematic in the weak-signal limit. There are also practical difficulties involved in human inspection of large numbers of spectra.

Line width and line centre are often derived from the same source, namely two measurements of the low- and high-velocity edges of a line. Line width is taken from their difference and line centre from their average. Bicay & Giovanelli (1986) present several typical methods of estimating the edge velocities. These all start by deciding a flux density value, then interpolate linearly between adjacent channels which straddle this value. The flux density chosen is some percentage of a characteristic value for the line as a whole, which may be the mean flux density, or the maximum within the line profile, or some function of the peaks at the horns of the profile, in cases where these can be measured. In all cases such values for the characteristic flux density depend upon either the channel width or on implicit assumptions about where the line begins and ends.

The largest blind H i survey (still on-going) is ALFALFA (Giovanelli et al. 2005); their most recent catalog release paper is Haynes et al. (2011). The latter authors describe their edge-estimation procedure as similar to that of Springob et al. (2005), in which polynomials (in practice nearly always “polynomials of order 1”, i.e. straight lines) are fitted to several channels at the line rolloff on each side, the interpolation flux density being 50% of the horn height on that side. In practice this is not much different to schemes already mentioned because, unlike the profiles shown in Fig. 2 of Springob et al. (2005), in ALFALFA the effective channel width of 10 km s^-1 would not provide many channels in the rolloff.

To settle on a point of comparison therefore we decided to employ, as our “traditional, non-parametric method” of estimating the line edges, the simple scheme as follows:

• 1.

  Determine s_max, the highest flux density within the spectrum.

• 2.

  Identify the number j_lo of the lowest-velocity channel whose s value exceeds 0.5s_max.

• 3.

  Identify the number j_hi of the highest-velocity channel whose s value exceeds 0.5s_max.

• 4.

  Interpolate between j_lo − 1 and j_lo to obtain v_lo, and between j_hi and j_hi + 1 to obtain v_hi.

It is thus a maximizing algorithm.

3.2.5. Monte Carlo results I: comparison with high-S/N fits

Some of the more important results of the Monte Carlo described in Sects. 3.2.2 to 3.2.4 are shown in Figs. 10 to 15.

In Fig. 10, the variation of the mean value of total flux fitted to the spectra is shown against the S/N, for all five galaxies in the trial. Some upward bias in the values becomes visible at low S/N values, but this always remains within 1 standard deviation of nominal.

Figure 11 is a similar plot, this time for the intrinsic line width. DDO 53 is not shown, because line width is small and poorly constrained for this galaxy, even in fitting to the raw THINGS spectrum. Similar slight biases are observed in three of the four galaxies; it is not clear why NGC 4214 seems to be so badly affected at low S/N.

It is important to note that, at the highest S/N levels shown, the linewidth clearly asymptotes to the value measured in Sect. 3.1 and tabulated in Col. 3 of Table 1, despite the fact that the spectral resolution of the present data is several times poorer than the original THINGS observations. This emphasises that no “instrumental broadening” correction to linewidth measures is necessary when fitting a line-profile model.

The curious under-estimation by up to 2σ of the linewidth for DDO 154 at quite large values of S/N remains unexplained, although it is worth noting that it is never more than about 8 km s^-1, which is less than the channel width of these spectra. In the next section it is demonstrated however that for DDO 154 the W₅₀ linewidth remains correct in this range – so clearly the decrease in the fitted “intrinsic” value is compensated by an increase in the “turbulent broadening” value.

Note that the mean fitted linewidth for all five galaxies remains accurate down to the ALFALFA limit of detectability.

The mean and standard deviation of the line centre offset (that is, the difference between the fitted line centre and the input value) are displayed in Fig. 12 just for DDO 154. All the galaxies show similar results. No significant bias is seen, and the scatter remains within the channel width for ALFALFA S/N values greater than 5.

From these results we can conclude that, so far as the important parameters of the fitted profile model go, there is little bias evident when the spectral resolution and S/N are lowered to typical levels expected in blind surveys.

3.2.6. Monte Carlo results II: comparison with non-parametric methods

Figures 13 through 15 all make the same comparison, but a single galaxy is plotted per figure for clarity. The other two galaxies have not been shown because their results are similar. For each value of S/N in each plot, two values of W50, the line width at 50% of maximum height, are shown: the error bar centred on a triangle comes from measuring the fitted profile; that centred on the square, from the raw spectrum, using the simple algorithm discussed in Sect. 3.2.4. It is obvious from all plots that finding W50 from the fitted profile yields a more accurate result. The effect is most stark for broad profiles, which for the same ALFALFA-formula S/N have a much lower peak-to-rms value.

3.2.7. Discussion

Lewis (1983) made a similar study of the effect of noise on measurements of H i line position and width. Although Lewis’s study was in many respects more comprehensive than the present one, this author applied traditional direct methods to measure the positions of line edges, and deprecated the use of functions fitted to the line shoulders over several channels. Lewis argued that the number of channels spanning the line shoulder was of the same order as the minimum number of parameters required by any useful fitting function, and therefore that fitting was equivalent merely to interpolation, and could therefore not be expected to yield any improvement in precision.

In opposition to this view is the Bayesian probability theory, which says that the most complete knowledge possible of the probability distribution of linewidth values is obtained by use of Bayes’ theorem in conjunction with a model of the spectral line. A seeming oversupply of model parameters is no objection provided the unwanted ones are marginalized out. The only necessary criterion is that the model reproduces accurately the widths of real spectral lines. For the present model, this is well demonstrated in Sect. 3.1.3.

Nevertheless it seems reasonable that a model with 6 parameters will be poorly constrained if the spectral line extends over no more than about the same number of velocity channels. It is of interest therefore to consider the likely shape of the Bayesian posterior distribution in the extreme case that the spectral line has most of its flux in a single channel. In this case it seems clear that the line centre and total flux parameters will remain well constrained. On the other hand, the fraction-solid and asymmetry parameters can be expected to be completely unconstrained, with the posterior density remaining about the same if they are varied through their respective ranges (as for example in Fig. B.5). Since the Gaussian spread of the spectral line is modelled by the Δv_rand parameter, its value will be most tightly constrained by the flux densities in the channels adjacent to the central one. If these are insignificant, then this parameter will be unconstrained within a range between zero and about the width of a channel. The final parameter, the intrinsic linewidth, will be expected to have the same character (this again is similar to what is seen in Fig. B.5).

In fact none of these indeterminacies and degeneracies matter, because in practice one integrates the posterior over all the dimensions (i.e. parameters) of no present interest. For example if one wants a probability distribution of W₅₀ for a given spectrum, this is very easy to obtain from the Bayesian formulation for the respective data. The number of parameters in the model should be viewed merely as a detail of the way the problem is worked out, having no importance for the final result.

Another issue which deserves some discussion concerns the biases in fit parameters observed in Figs. 10 to 15. If the Bayesian formulation is optimum, one might ask, why does it still lead to biased results? The answer is that the Monte Carlo procedure adopted in the present section is not the way one ought to proceed with real data. In no case is it optimum to fit separate spectra, then average the resulting parameter values. This is done here only because we have no option if we wish to make comparison with standard methods of estimating line parameters directly from the data. If one did not have the necessity of comparison with non-Bayesian methods, the correct approach depends on the nature of the quantity to be estimated. Either we have an ensemble of observations of the same object, or an ensemble of observations of different objects. In the former case, the correct Bayesian approach is to fit once to the entire data set. In this case we are back (given sufficient data) in the high S/N regime and any biases will be correspondingly insignificant. If the objects are all different, then the correct approach is hierarchical modelling. An interesting recent example of this treatment is Brewer et al. (2014). Essentially one defines a hyper-model which describes the distribution of values of properties across the ensemble of objects, and tries to constrain its hyper-parameters. Although any individual spectrum is low S/N, again a sufficiently large data set will eventually reach the high-S/N, thus low-bias regime as regards the hyperparameters.

3.3. The THINGS galaxies in EBHIS

3.3.1. Overview

The semi-simulated spectra generated in Sect. 3.2 still fall some way short of real life and therefore cannot be expected to present the full range of difficulties one encounters in trying to analyse real observations. We wanted also to test the profile model on single-dish data. The EBHIS survey (Effelsberg Bonn H i Survey, Winkel et al. 2010; Kerp et al. 2011) is convenient for this, since its survey coverage has included most of the THINGS galaxies, yet as a single-dish survey, it necessarily has much reduced spatial resolution when compared to the interferometer used by Walter et al. (2008) to make their THINGS observations. An extensive comparison of the EBHIS observations of the THINGS galaxies might be interesting and useful, but is beyond the scope of the present paper. Here we just wish to show that profile-fitting works well also under less than ideal conditions. Specifically we are interested in the following questions:

• Can we cope with non-flat baselines and interference?

• Can profile-fitting at separate pixels across the source return a map of the source flux which is better (less noisy) than one obtained simply by adding flux densities over a range of channels?

• Can the variation in the asymmetry parameter across the source give clues to the galaxy orientation and inclination, even with poorly resolved sources?

3.3.2. Baseline fitting

The baseline problem differs between single-dish vs. interferometer spectra. Single-dish measurements tend to be worse affected by Fabry-Perot type effects which can produce waves in spectral baselines. Spectra made with an interferometer on the other hand tend to have flat baselines, although a spatially superimposed continuum source can raise and otherwise perturb the baseline. Future, deep interferometric surveys for H i may, however, begin to be affected by source confusion, which may necessitate wider use of baseline fitting.

Traditionally, non-flat baselines have been dealt with by fitting a function, often a simple polynomial, to stretches of channels on either side of the line of interest. The fit is interpolated across the line profile but channels within the profile make no contribution to the fit. This is less than ideal for three reasons: firstly because of the restricted number of channels which contribute information to the fit; secondly because it relies on a decision as to where the profile begins and ends; thirdly, selection of the degree of the fitting polynomial is usually done “by eye”.

The approach taken in the present paper is simply to add a baseline function to the existing profile model, and fit the combined model to all selected channels. This also lends itself to a Bayesian approach to determining the best order of the baseline function. This technique is more fully described in Appendix B.6.

Chebyshev polynomials have been preferred for the baseline function in the present paper, because they are easy to generate, and more orthogonal and numerically tractable than simple polynomials. Note however that a Chebyshev polynomial is still just a polynomial, and an nth-order Chebyshev fitted to some data will result in exactly the same function as fitting a simple nth-order polynomial would (B. Winkel, priv. comm.). The Chebyshev is just better behaved from a computational point of view.

The EBHIS cubes have 3000 spectral channels spanning nearly 3900 km s^-1, which is far wider than any of the lines fitted. It isn’t necessary to fit to the full velocity range, so in all cases a small section only, about 3 or 4 times wider than the spectral line, was selected for fitting purposes.

3.3.3. Results

Results are presented here for only a single galaxy, DDO 154.

Figure 16 shows the spectrum of DDO 154 obtained by integrating a 9 × 9 pixel area of the relevant data cube. The profile shown is the mean of 10⁴ iterations of a converged MCMC, but it is almost indistinguishable to the eye from the best-fit profile from a Levenberg-Marquardt optimisation. The baseline is fitted by a sum of Chebyshev polynomials truncated at order 5 (i.e. giving 6 orders in total), which number gave the maximum Bayesian evidence (see Appendix B.6).

A spectrum summed across the entire breadth of this cube showed no evidence for narrow-band RFI in any channel, so no channels were excised before fitting.

Because there is a significant baseline contribution in this spectrum, careful attention was paid to baseline priors, which were estimated from the non-source portions of the cube as described in Sect. B.7. As can be seen from the figure, the prior accounts for quite a lot of the baseline over the small spatial extent of the source. This gives one confidence that the baseline subtraction is accurate, and thus also the fitted total flux of the spectral line.

A Levenberg-Marquardt optimization returned a best-fit value of 102 ± 4 Jy km s^-1 for the total flux of DDO 154, whereas the MCMC returned a mean value of 100 ± 5 Jy km s^-1. These values are significantly higher than the value of 82 fitted to the data of Walter et al. (2008). However, it is known that interferometric measurements can miss flux. Our value agrees very well with the value of 105 ± 5 Jy km s^-1 reported by Carignan & Purton (1998) from a careful combination of interferometer and single-dish data.

Next, the baseline-plus-line model was fitted to spectra extracted at single spatial pixels across the same 9 × 9 area integrated to produce the spectrum in Fig. 16. 9 × 9 maps of the fitted parameter values were made. Figures 17 and 18 show two of these, respectively for the total-flux parameter and the asymmetry parameter. (The pixels of Fig. 18 were set to null where the total flux fell below 15% of maximum.) It is of interest to compare these to the much higher-resolution versions in Fig. 53 in Walter et al. As one might expect, the single-dish EBHIS survey seems to detect more extended structure than Walter et al’s interferometric measurements.

Clearly visible in the asymmetry map is the orientation of the galaxy’s major and minor axes, even though DDO 154 is barely resolved by the 9 arcmin FWHM beam of the Effelsberg telescope at L band (Kerp et al. 2011).

The remaining parameter maps are of less interest and are not shown. The map of the fitted line centre also (as one might expect) shows indications of the major and minor axes, but the signal is not so unequivocal as with the asymmetry map.

As far as the baseline fits go, no obvious trends are visible either in the 6 maps of the individual Chebyshev coefficients, or a map of their rms value (Fig. 19). There is no bright continuum radio source to disturb the baseline at the location of DDO 154: NVSS (Condon et al. 1998) shows that there is no source with an average flux density at L band greater than 90 mJy/beam within the approximately 30′′ by 30′′ square shown in Figs. 17 to 19. The first Chebyshev coefficient (which is just a DC offset) does appear slightly elevated at the centre of the source, but only at the 2-sigma level.

4. Conclusion

The two most desirable measurements of the global H i spectral profile of a galaxy are its width and area. Combined with optical measurements of the inclination and brightness of the galaxy, these allow the galaxy’s H i mass and intrinsic optical luminosity, and thence its distance, to be estimated via the baryonic Tully-Fisher relation. The systemic velocity of the galaxy comes next in importance. This quantity, when compared to the Hubble velocity appropriate to that distance, additionally allows one to calculate the peculiar or local velocity of the galaxy in the line of sight to the observer. A profile model is described here which permits one to make accurate estimates of these quantities.

In Sect. 3.1 the model was fitted to a variety of H i spectra from the THINGS survey (Walter et al. 2008). These galaxies were observed at high velocity resolution and under very low-noise conditions, and exhibited negligible baselines.

The fits are nearly always excellent. Local densities and voids in the H i distribution give rise to noisiness in the profiles, but residuals from these wiggles typically have amplitudes less than about 10% of the profile height. Such local deviations seem to have little effect on the bulk fitted properties.

Particularly important to note is the close fit to the slopes of the lines, which means that linewidth parameters of the profile such as W₅₀ or W₂₀ are accurately (within 2%) reproduced by the best-fit model. There is a great deal of existing lore on how to correct such linewidth measurements so as to obtain estimates of the maximum rotation speed of the galaxy, which is the quantity desired for Tully-Fisher calculations. The availability of a good model proxy for the linewidth allows one to make use of such formulae without radical modification, even in cases of low spectral resolution or S/N.

These correction formulae include a correction for instrumental broadening. Strictly speaking, this is no longer necessary if the linewidth is taken from the model profile, since this profile is what would actually be seen with no instrumental broadening (i.e., infinite resolution). Trials in Sect. 3.2 in which the model is fitted to spectra with coarser velocity binning indeed show no significant “instrumental” bias in linewidth in the high-S/N limit. Biases do begin to show themselves as S/N decreases toward the limit of detectability, but they are always significantly less than those observed with linewidths measured directly from the noisy, coarsely-resolved profiles. As described in Sect. 3.2.7, correct Bayesian treatment of ensemble observations can circumvent ill effects from such low-S/N biases.

Comparison of the total flux parameter of the fitted model to that obtained by summing the observed flux densities indicates again that the model value of flux is an accurate proxy to the real one (within 1%). For about half the 34 galaxies fitted, the difference was less than the measurement uncertainty, which was nearly always significantly less than 1%. There is however a significant minority of galaxies for which the total flux is slightly underestimated. It is thought likely that this is because of the frequent occurrence of extended tails in the profile wings, which are underfitted by the Gaussian wings of the model.

Systemic velocities are returned within measurement errors and without bias by the model fits, even for asymmetric profiles.

The use of the remaining three model parameters (turbulent broadening Δv_rand, “solid rotating” fraction f and asymmetry α) has been little explored. The Δv_rand value fitted is typically about 12 km s^-1 (see Appendix B.2), which is 2 or 3 km s^-1 larger than typical values seen in spatially resolved maps of velocity distribution. As explained in Sect. 3.1.3, this is likely an artefact of the simplistic rotation curve implied by the model. Deviations of the rotation curves of real galaxies from this simple framework are readily absorbed into the Δv_rand parameter, tending therefore to inflate its value.

It might be of interest to see whether the f and α parameters can be useful as proxies for features of galaxy morphology, in analogy with studies such as those of Richter & Sancisi (1994) and Andersen & Bershady (2009), but this question is beyond the scope of the present paper.

In Sect. 3.2 it was shown that the model remains a useful way to extract parameters of interest even when the spectral resolution and S/N more nearly approximate those of the bulk of detected sources in surveys such as ALFALFA.

The final test performed on the model was to fit to a representative H i profile from the EBHIS survey (Winkel et al. 2010; Kerp et al. 2011). In contrast to the THINGS observations, this survey was done using a single dish antenna, and thus comes with a significant amount of baseline. This was dealt with in a Bayesian fashion by using the non-source areas of the relevant data cube to constrain the baseline priors. The success of this procedure is shown by the accurate agreement between the resulting fitted value of total flux and previous careful measurements for this galaxy (DDO 154).

Fitting the profile to individual pixels instead of to a spatial sum across several pixels allowed something of the spatial distribution of neutral hydrogen in this galaxy to be mapped in a low-noise fashion, despite the intrinsically poor spatial resolution of the survey. The orientation of the galaxy’s rotation axis could also be approximately determined via mapping the fitted value of the asymmetry parameter. The success of these fits shows that the model is not only suitable for H i profiles which result from a complete spatial integration across the extent of the galaxy, but also for partially integrated profiles.

---

Acknowledgments

I.M.S. wishes to thank F. Walter and J. Kerp for generous provisions of non-public data, and is grateful to L. Flöer, A Schröder and B. Winkel for help and advice. We are grateful to the anonymous referee for helpful and informative comments.

References

Appendix A: The model in Fourier space

Appendix A.1: The transform components

As mentioned in Sect. 2.3, the model is generated in Fourier space partly to avoid singularities, partly to facilitate convolution, and partly to mimic the process by which real radio spectra are generated in an XF-type correlator. To follow the latter process requires the accurate construction of a synthetic discrete autocorrelation function Z_j for j = −N to N − 1. This is done by sampling the Fourier transform Z(τ) of the profile model at appropriate intervals in τ, which is proportional to the autocorrelation lag time. The correct autocorrelation values result, provided only that the spectral line to be modelled has no significant power outside the chosen velocity interval [ v_lo,v_hi ]. This interval needs to be decided on before calculating the transforms. Together with the number of channels N, it defines the channel width Δv_chan.

The Fourier transform Z(τ) of the entire profile model is obtained by multiplying together the separate transforms Z_intrinsic(τ) and G(τ) of the expressions respectively for s_intrinsic(v) and the dispersive convolver g(v) given in Eqs. (2) and (3). The transform of the intrinsic or undispersed profile evaluates to $\begin{array}{ccc}{\mathit{Z}}_{\mathrm{intrinsic}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}& \mathrm{=}& \mathrm{2}\mathit{S}\exp \mathrm{\left(}\mathrm{i}\mathit{\psi \tau }\mathrm{\right)}[\left\{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}\mathit{f}}{\mathit{\tau }}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right\}\\ & & \mathrm{+}\hspace{0.17em}\mathrm{i}\mathit{\alpha }\left\{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}\mathit{f}}{\mathit{\tau }}\left[\frac{\mathrm{2}}{\mathit{\tau }}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{-}{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right]\right\}]\end{array}$(A.1)where J₀ and J₁ are respectively the zeroth and first order Bessel functions of the first kind, and

$\mathit{\psi }\mathrm{=}\frac{\mathrm{2}\mathrm{(}{\mathit{v}}_{\mathrm{hi}}\mathrm{-}{\mathit{v}}_{\mathrm{ctr}}\mathrm{)}}{\mathrm{\Delta }\mathit{v}}\mathrm{·}$The transform of the dispersive-motion convolver gives

$\mathit{G}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\exp \left[\mathrm{-}\mathrm{2}{\left(\frac{\mathrm{\Delta }{\mathit{v}}_{\mathrm{rand}}}{\mathrm{\Delta }\mathit{v}}\mathit{\tau }\right)}^{\mathrm{2}}\right]\mathrm{·}$The correct sampling is at τ = jΔτ for integer j ∈ [ − N,N − 1 ], where

$\mathrm{\Delta }\mathit{\tau }\mathrm{=}\frac{\mathit{\pi }\mathrm{\Delta }\mathit{v}}{\mathrm{2}\mathit{N}\mathrm{\Delta }{\mathit{v}}_{\mathrm{chan}}}\mathrm{·}$Recall that Δv in all these expressions is the linewidth parameter of the model, as described in Sect. 2.2.

Appendix A.2: Computing close to singularities

Note that we did not group all the J₀ and J₁ terms together in Eq. (A.1). This is because terms with τ in the denominator need special treatment as τ → 0.

The expression J₁(τ) /τ in Eq. (A.1) is analogous to the sinc function. For values of τ close to zero it should be approximated by its power series

$\frac{{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}{\mathit{\tau }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{8}}\mathrm{+}\frac{{\mathit{\tau }}^{\mathrm{4}}}{\mathrm{192}}\mathrm{-}\mathit{...}\hspace{0.17em}\right)\mathit{.}$Let us denote the last set of terms in Eq. (A.1) as

$\mathit{E}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\tau }}\left[\frac{\mathrm{2}}{\mathit{\tau }}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{-}{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right]\mathit{.}$The components of this function also have singularities at τ = 0 which cause computational difficulties as τ → 0. The power series for E evaluates to

$\mathit{E}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathit{\tau }}{\mathrm{8}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{12}}\mathrm{+}\frac{{\mathit{\tau }}^{\mathrm{4}}}{\mathrm{384}}\mathrm{-}\mathit{...}\hspace{0.17em}\right)\mathit{.}$

Appendix A.3: Derivatives

It can be convenient (for example when doing Levenberg-Marquardt optimization) to have expressions for the derivatives of the model with respect to each of the six parameters. We give the Fourier transforms of these here.

For convenience we define C(τ) such that the whole-model transform can be expressed as

$\mathit{Z}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathit{Z}}_{\mathrm{intrinsic}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{\times }\mathit{G}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{SG}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\exp \mathrm{\left(}\mathrm{i}\mathit{\psi \tau }\mathrm{\right)}\mathit{C}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathit{.}$C thus contains all terms within the outer square bracket of Eq. (A.1). The derivatives can then be given as $\begin{array}{ccc}& & \frac{\mathit{\partial Z}}{\mathit{\partial }{\mathit{v}}_{\mathrm{ctr}}}\mathrm{=}\frac{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{\tau }}{\mathrm{\Delta }\mathit{v}}\mathit{Z}\\ & & \frac{\mathit{\partial Z}}{\mathit{\partial }\mathrm{\Delta }\mathit{v}}\mathrm{=}\mathrm{2}\mathit{SG}\exp \mathrm{\left(}\mathrm{i}\mathit{\psi \tau }\mathrm{\right)}\frac{\mathit{\partial C}}{\mathit{\partial }\mathrm{\Delta }\mathit{v}}\\ & & \frac{\mathit{\partial Z}}{\mathit{\partial S}}\mathrm{=}\mathit{Z}\mathit{/}\mathit{S}\\ & & \frac{\mathit{\partial Z}}{\mathit{\partial }\mathrm{\Delta }{\mathit{v}}_{\mathrm{rand}}}\mathrm{=}\mathrm{-}\mathrm{4}\mathrm{\Delta }{\mathit{v}}_{\mathrm{rand}}{\left(\frac{\mathit{\tau }}{\mathrm{\Delta }\mathit{v}}\right)}^{\mathrm{2}}\mathit{Z}\\ & & \frac{\mathit{\partial Z}}{\mathit{\partial f}}\mathrm{=}\mathrm{2}\mathit{SG}\exp \mathrm{\left(}\mathrm{i}\mathit{\psi \tau }\mathrm{\right)}\frac{\mathit{\partial C}}{\mathit{\partial f}}\\ & & \frac{\mathit{\partial Z}}{\mathit{\partial \alpha }}\mathrm{=}\mathrm{2}\mathit{SG}\exp \mathrm{\left(}\mathrm{i}\mathit{\psi \tau }\mathrm{\right)}\frac{\mathit{\partial C}}{\mathit{\partial \alpha }}\mathrm{·}\end{array}$The C derivatives are as follows: $\begin{array}{ccc}\frac{\mathit{\partial C}}{\mathit{\partial }\mathrm{\Delta }\mathit{v}}& \mathrm{=}& \frac{\mathrm{-}\mathit{\tau }}{\mathrm{\Delta }\mathit{v}}\{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\mathrm{2}\mathit{fE}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\\ & & \mathrm{-}\hspace{0.17em}\mathrm{i}\mathit{\alpha }\left[\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}\left({\mathit{J}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{-}\frac{{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}}{\mathit{\tau }}\right)\mathrm{+}\frac{\mathrm{2}\mathit{f}}{\mathit{\tau }}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{-}\frac{\mathrm{6}\mathit{f}}{\mathit{\tau }}\mathit{E}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right]\}\\ \frac{\mathit{\partial C}}{\mathit{\partial f}}& \mathrm{=}& \left[\mathrm{-}{\mathit{J}}_{\mathrm{0}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\frac{\mathrm{2}}{\mathit{\tau }}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right]\mathrm{+}\mathrm{i}\mathit{\alpha }\left[\mathrm{-}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\mathrm{2}\mathit{E}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right]\\ \frac{\mathit{\partial C}}{\mathit{\partial \alpha }}& \mathrm{=}& \mathrm{i}\left[\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}\mathrm{)}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\mathrm{2}\mathit{fE}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\right]\mathit{.}\end{array}$

Appendix B: Model fitting – technical issues

Appendix B.1: Priors in general

Bayes’ theorem (Eq. (5)) describes the posterior probability distribution of the model parameter values for a given set of data. This represents all that can be known about the parameters. However usually it suffices to settle for less than this. Typically one wishes to approximate the posterior via a single set of optimum parameter values coupled with widths of the posterior function in each dimension of the parameter space. This represents a full description of the posterior only if this function is an unskewed Gaussian; but it is often a good enough approximation. “Optimum” can mean the values of the parameters at the maximum of the posterior, but mean values are often an acceptable proxy for these. Only if the posterior is very asymmetric will these two sets of values be significantly different.

The posterior is a product of the likelihood and a function which represents prior knowledge of the parameter distribution. It is often easier to specify the likelihood than the priors. A typical example is where nothing definite at all is known about the likely distribution of the parameter values. What form of prior best encodes such a state of ignorance? There has been considerable debate about the correct forms for “ignorance” priors (see for example references in D’Agostini 2003). The question depends however to some degree on what one wants to do with the posterior. If a value for the Bayesian evidence is what is needed, then priors have to be more carefully chosen: in particular, they must be integrable. Priors which are not integrable are named “improper”. However in most cases discussed in the present paper, the location of the maximum in the posterior is all that is desired. If in addition the S/N of the data is high, the shape of the posterior will be dominated by a sharply-peaked likelihood, and the exact form of the priors becomes unimportant. One may even choose an improper prior without disadvantage.

Unbounded, thus improper, priors were used for most of the parameters in fitting the profile model to the original, high-S/N THINGS data. The fraction-solid and asymmetry parameters f and α have natural high and low bounds, but all the other parameters are unbounded on one or both sides.

Appendix B.2: Priors chosen for the low S/N spectra

The distributions of fitted parameters from the high-S/N fits were used to generate some priors for the low-S/N fits. Identical priors were applied in fitting both the semi-simulated spectra in Sect. 3.2 and the EBHIS spectra in Sect. 3.3.

The spectral-line fits performed on the 34 THINGS galaxies in Sect. 3.1 don’t provide a very numerous sample, but some trends do become apparent. Using this information, non-flat priors were decided for four of the six model parameters, namely the linewidth Δv, the dispersion velocity width Δv_rand, the H i mass fraction in the core f, and the asymmetry parameter α. Histograms of the distribution of these four across the 34 THINGS fits are shown in Figs. B.1–B.4.

A logarithmic histogram of the fitted linewidth parameter Δv is shown in Fig. B.1. Despite the single galaxy (DDO 53) with linewidth close to 10 km s^-1, a log-normal prior centred on 2.2 (i.e., at 10^2.2 km s^-1) and of width 0.3 seemed a reasonable fit to the THINGS values.

Figure B.2 shows the dispersion width. As already discussed, this is typically wider than would be expected from turbulence and thermal velocities alone. For this parameter we chose a log-normal prior function with a center in log space of 1.17 and a width of 0.14. It isn’t a particularly good fit to the histogram however: something with a sharper low-end cut-off and a longer tail would be better.

In Fig. B.3 a similar histogram for the core-fraction parameter is presented. The decrease in the incidence of this parameter as it increases in value is consistent with an exponentially decreasing prior, although a value of 0.3 was chosen for the characteristic constant so as to be consistent with the presence of the outlier in the histogram at the high end of the scale.

The final histogram in Fig. B.4 shows the asymmetry parameter. Values which depart from the midpoint by more than about 0.2 are seen to be unusual. The outliers here occur for spectra for which the model is not a good fit to the profile. A Gaussian with sigma equal to 0.2 was chosen for this prior.

In addition to these four priors, and purely in order to exclude nonsense results due to fitting the profile to noise away from the true spectral line, we also chose a prior for the line centre parameter, this being a Gaussian centred on the true line centre (as determined from the high-S/N fits), and with a width somewhat arbitrarily set equal to 20 km s^-1. This is considered acceptable because we assume that in a real survey, the rough location of the line would already been found by other means. We remind the reader that we are not treating source detection here, but source characterisation.

For the remaining total-flux parameter, a flat prior was retained, with non-physical, negative-valued results excluded as before.

To summarize, the priors adopted in the present section were:

$\begin{array}{c}\\ \mathit{x}& \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ {\mathit{v}}_{\mathrm{ctr}}& \exp \mathrm{(}\mathrm{-}\mathrm{0.5}\mathrm{[}\mathrm{\{}\mathit{x}\mathrm{-}{\mathit{v}}_{\mathrm{true}}\mathrm{\}}\mathit{/}\mathrm{20}{\mathrm{]}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{20}\sqrt{\mathrm{2}\mathit{\pi }}\mathrm{)}\\ \mathrm{\Delta }\mathit{v}& \exp \mathrm{(}\mathrm{-}\mathrm{0.5}\mathrm{[}\mathrm{\left\{}\log \mathrm{10}\mathrm{\right(}\mathit{x}\mathrm{)}\mathrm{-}\mathrm{2.2}\mathrm{\}}\mathit{/}\mathrm{0.3}{\mathrm{]}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{0.3}\ln \mathrm{10}\sqrt{\mathrm{2}\mathit{\pi }}\mathrm{)}\\ \mathit{S}& \mathit{x}\mathit{>}\mathrm{0}\\ \mathrm{\Delta }{\mathit{v}}_{\mathrm{rand}}& \exp \mathrm{(}\mathrm{-}\mathrm{0.5}\mathrm{[}\mathrm{\left\{}\log \mathrm{10}\mathrm{\right(}\mathit{x}\mathrm{)}\mathrm{-}\mathrm{1.17}\mathrm{\}}\mathit{/}\mathrm{0.14}{\mathrm{]}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{0.14}\ln \mathrm{10}\sqrt{\mathrm{2}\mathit{\pi }}\mathrm{)}\\ \mathit{f}& \mathrm{0.3}\exp \mathrm{(}\mathrm{-}\mathrm{0.3}\mathit{x}\mathrm{)}\\ \mathit{\alpha }& \exp \mathrm{(}\mathrm{-}\mathrm{0.5}\mathrm{[}\mathit{x}\mathit{/}\mathrm{0.2}{\mathrm{]}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{0.2}\sqrt{\mathrm{2}\mathit{\pi }}\mathrm{\left)}\mathrm{for}\mathrm{\right|}\mathit{x}\mathrm{|}\mathrm{\le }\mathrm{1}\mathit{,}\mathrm{=}\mathrm{0}\mathrm{else}\mathrm{.}\end{array}$

Appendix B.3: Possible fitting problems

Appendix B.3.1: Degeneracy

The chosen set of parameters are convenient for many reasons but have the problem of not being orthogonal under all circumstances. For example, if the intrinsic linewidth Δv is small compared to the dispersion width Δv_rand, the f and α parameters begin to be poorly constrained. Figure B.5 serves to illustrate this point. And for at least one of the fits in Sect. 3.1, the best-fit value of α is close to unity, which causes nearly all the other parameters to be poorly constrained. These degeneracies can occur for both scenarios described in Sect. 2.4.1 and can cause convergence problems when trying to fit the profile. An MCMC is again the preferred way to deal with such cases.

Appendix B.3.2: Lag-space problems

Filtering in lag space in general causes the noise values in adjacent velocity channels to be correlated. The only common filtering scenario which avoids this is Hanning filtering followed by discarding every second channel. This was the scheme employed by Walter et al. (2008), so we have not had to address this problem. Other data sets are not so favourable though: for example, HIPASS made use of a 25% Tukey filter (Barnes et al. 2001) which introduces correlations between adjacent channels. The way to deal with this is called generalized least squares and involves reformulating the expression for χ² to have the form ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{s}\mathrm{-}{m}{\mathrm{)}}^{\mathrm{T}}{{C}}^{-1}\mathrm{(}{s}\mathrm{-}{m}\mathrm{)}\mathit{.}$(B.1)Here the T superscript indicates transpose and C is the covariance matrix of the data channels. For uncorrelated data C is of course diagonal and Eq. (B.1) reverts to the standard expression. Equation (B.1) seems straightforward but in fact C is often singular for correlated data. This is effectively because there are now fewer than N degrees of freedom. Hanning smoothing for example reduces the number of degrees of freedom by a factor of 2, which is why one usually discards half the channels afterwards – they don’t contain any new information. In general the solution to this is not to try to invert C but to calculate instead its pseudo-inverse (Penrose 1955), which reduces the rank of the problem to the proper number. But this is not further elaborated here.

An alternative approach to the case of correlated noise would be to back-transform the data and perform the fitting in lag space, where there is no correlation between lag samples. This is a little complicated if the velocity channels have different standard deviations, which is the case in Sect. 3.1, so it has not been attempted here.

When fitting to a real spectral line, in principle both the velocity range and the number of channels should match that in the original observation. This is because the end points where the autocorrelation function is truncated (particularly the low-frequency, i.e. high-velocity end) generate some ringing at places in the spectrum where its slope is steep (see e.g. Willis & Bregman 1993). If there is such ringing in the data, it’s obviously advantageous if the model reproduces this artefact.

There were two reasons why, for many of the THINGS galaxies, this was difficult to do. Firstly, the cubes published by Walter et al. (2008) don’t include all the observed channels. The starting channels had to be estimated from the tabulated description of the observations. Secondly, for some of the targets, Walter et al spliced together several observations taken with different velocity ranges. These cannot be accurately represented via a single simulated autocorrelation. These inaccuracies were however not regarded as significant.

Appendix B.4: Markov chain Monte Carlo

A Markov chain Monte Carlo (MCMC) is a sequence of steps within the parameter space. According to the Metropolis-Hastings algorithm (Metropolis et al. 1953), at each step a random displacement from the last good position is generated, and a test is applied to the value of the posterior at the resulting trial position. If the test is passed, the trial position is accepted as the next position in the chain; if not, a new trial is generated. It can be shown that the distribution of accepted positions will in time converge to that of the posterior.

This procedure is simple and robust to apply, but practical issues can arise, principally two: the speed of convergence, and how to know when the MCMC has converged. If the MCMC starts a long way from the centre of the posterior, and if the dimensions of the step size distribution are very different from those of the posterior, the rate of acceptances can be very low. This is the reason for using a LM fit in Sect. 3.1 to find the location and dimensions of the posterior maximum before starting the MCMC. After starting the MCMC, convergence speed was further improved via an iterative procedure, described as follows. The MCMC was divided into a sequence of separate chains. To begin with, the chains were short - only tens to hundreds of accepted steps were required. After finishing each chain, the covariance of the parameter points was calculated and used to specify the shape of an N-dimensional Gaussian step-size function for use in the next chain (N being the number of model parameters). The chains thus gradually converged at a faster rate, and more acceptances could be demanded. The final sequence of chains used had

${\mathit{N}}_{\mathrm{A}}\mathrm{=}\mathrm{\left\{}\mathrm{30}\mathit{,}\mathrm{100}\mathit{,}\mathrm{200}\mathit{,}\mathrm{1000}\mathit{,}\mathrm{1000}\mathit{,}\mathrm{10000}\mathrm{\right\}}\mathit{,}$where N_A is the number of acceptances required.

For a Gaussian posterior in a 6-dimensional parameter space, the fraction of accepted steps is expected to be a little less than 0.3 (estimated from Figs. 3 and 4 of Hanson & Cunningham 1998). This figure may be compared with the distribution of this fraction among the THINGS galaxies shown in Fig. B.6.

The parameter values and uncertainties listed in Table 1 were obtained by calculating the means and covariances of the final chain of 10⁴ acceptances.

Others who employ the MCMC method commonly use far more points in the final chain than 10⁴; typically also convergence tests are applied; sometimes multiple chains are used, with sparse sampling of the final chains, in order to avoid the obvious correlation between successive points. We have not felt any of this to be necessary. Convergence was judged by monitoring the acceptance rate from chain to chain in the sequence, also by examining some plots of the successive parameter values. When no drift in the mean could be detected, it was judged that these chains have converged enough. The sequence given above was found to be adequate for all 34 profiles. Correlation was thought not to be a problem, since it is the ensemble of points which matters to us, and if there are enough points, then any sparse subset of them ought to have the same distribution as any other, and thus of the whole.

Appendix B.5: Calculation of the spectrum uncertainties

The THINGS spectra were formed by summing, for each velocity channel, the values of a subset of pixels in the corresponding plane of the data cube. The subset was chosen by Walter et al. (2008) to include only those pixels where some H i emission was detectable above the background. The uncertainty in each summed value can be expected to scale with the square root of the number of pixels summed, but calculating an exact value has some complications, because the noise values in adjacent spatial pixels are already correlated via convolution by the interferometer beam.

It is easiest to examine this in the Fourier domain. Given input uncorrelated noise of standard deviation σ which is then convolved by a beam b_j,k, the standard deviation or rms σ_I of the output is given by$\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{I}}^{\mathrm{2}}& \mathrm{=}& {\mathit{\sigma }}^{\mathrm{2}}\mathrm{⟨}{B}{{B}}^{\mathrm{\ast }}\mathrm{⟩}\\ & \mathrm{=}& {\mathit{\sigma }}^{\mathrm{2}}\mathrm{\Sigma }{\mathit{b}}^{\mathrm{2}}\mathrm{or}{\mathit{\sigma }}^{\mathrm{2}}{b}\mathrm{·}{b}\end{array}$(B.2)where B is the discrete Fourier transform of the beam b, ∗ signifies complex conjugation, and the angle brackets represent the mean value.

This is the state of the data in each image plane of the cube as we receive it. Note that in interferometry practice it is standard to normalize the beam such that b_0,0 = 1. This allows one to read the flux density of point sources from the point image value I_j,k at their locations, and gives rise to the expression “janskys per beam” for the units of the image values. For an extended source, the total flux density S can be approximated by $\mathit{S}\mathrm{~}\frac{\mathrm{1}}{{\mathit{A}}_{\mathrm{beam}}}\sum _{\mathit{j,k}}^{{\mathit{A}}_{\mathrm{patch}}}{\mathit{I}}_{\mathit{j,k}}$(B.3)where A_patch represents the subset of image pixels summed over, and A_beam is the “area” in pixels of the beam, defined by

${\mathit{A}}_{\mathrm{beam}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{b}}_{\mathrm{0}\mathit{,}\mathrm{0}}}\begin{array}{c}\sum \end{array}{\mathit{b}}_{\mathit{j,k}}\mathit{.}$The problem we have is to calculate the uncertainty in the S values given by Eq. (B.3), starting with a knowledge of b and the value of σ_I. This becomes easy once we realize that Eq. (B.3) can be thought of as a further convolution

${S}\mathrm{=}\frac{{I}\mathit{\star }{p}}{{\mathit{A}}_{\mathrm{beam}}}\mathit{,}$where the convolver p equals unity for pixels within the patch or subset contributing to the sum, and zero otherwise. Bold font has been used to emphasize that all these quantities are images; the desired scalar value of spectrum flux density is in this notation given by S_0,0. The uncertainty in S can thus be related directly to the original σ by

${\mathit{\sigma }}_{\mathit{S}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\sigma }}^{\mathrm{2}}}{{\mathit{A}}_{\mathrm{beam}}^{\mathrm{2}}}\mathrm{⟨}\mathrm{\left(}{B}{{B}}^{\mathrm{\ast }}\mathrm{\right)}\mathrm{\left(}{P}{{P}}^{\mathrm{\ast }}\mathrm{\right)}\mathrm{⟩}\mathit{,}$where P is the discrete Fourier transform of p. More directly useful to us is to make use of Eq. (B.2) to express σ_S in relation to the rms σ_I of the cube planes:

${\mathit{\sigma }}_{\mathit{S}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathit{I}}^{\mathrm{2}}}{{\mathit{A}}_{\mathrm{beam}}^{\mathrm{2}}}\frac{\mathrm{⟨}\mathrm{\left(}{B}{{B}}^{\mathrm{\ast }}\mathrm{\right)}\mathrm{\left(}{P}{{P}}^{\mathrm{\ast }}\mathrm{\right)}\mathrm{⟩}}{\mathrm{\Sigma }{\mathit{b}}^{\mathrm{2}}}\mathit{.}$This can be calculated for any given p, but it is of interest to look at its behaviour when A_patch ≫ A_beam. In this case, we expect P to be much more compact than B, thus we can approximate $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{S}}^{\mathrm{2}}& \mathrm{~}& \frac{{\mathit{\sigma }}_{\mathit{I}}^{\mathrm{2}}}{{\mathit{A}}_{\mathrm{beam}}^{\mathrm{2}}}\frac{\mathrm{(}{\mathit{B}}_{\mathrm{0}\mathit{,}\mathrm{0}}{{\mathit{B}}_{\mathrm{0}\mathit{,}\mathrm{0}}^{\mathrm{\ast }}}^{\mathrm{)}}\mathrm{⟨}{P}{{P}}^{\mathrm{\ast }}\mathrm{⟩}}{\mathrm{\Sigma }{\mathit{b}}^{\mathrm{2}}}\\ & \mathrm{~}& \frac{{\mathit{\sigma }}_{\mathit{I}}^{\mathrm{2}}}{{\mathit{A}}_{\mathrm{beam}}^{\mathrm{2}}}\frac{{\mathrm{(}\mathrm{\Sigma }{\mathit{b}}^{\mathrm{)}}}^{\mathrm{2}}{\mathit{A}}_{\mathrm{patch}}}{\mathrm{\Sigma }{\mathit{b}}^{\mathrm{2}}}\mathrm{·}\end{array}$For a Gaussian b,

$\mathrm{\Sigma }{\mathit{b}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{b}}_{\mathrm{0}\mathit{,}\mathrm{0}}}{\mathrm{2}}\mathrm{\Sigma }\mathit{b}\mathit{.}$So we arrive at a final approximation

${\mathit{\sigma }}_{\mathit{S}}\mathrm{~}\frac{{\mathit{\sigma }}_{\mathit{I}}}{{\mathit{b}}_{\mathrm{0}\mathit{,}\mathrm{0}}}\sqrt{\frac{\mathrm{2}{\mathit{A}}_{\mathrm{patch}}}{{\mathit{A}}_{\mathrm{beam}}}}\mathrm{·}$Thus, with b_0,0 set to 1 according to usual interferometer practice, the ratio between the output and input uncertainties σ_S and σ_I is given, not by the square root of the “number of beams” A_patch/A_beam, but by this quantity times root 2.

Appendix B.6: Deciding the order of the baseline model

In fitting to the EBHIS profiles we had to deal with spectra with non-flat baselines. For this purpose a sequence of Chebyshev polynomials T_j was added to the model, as described in Sect. 3.3:

${\mathit{s}}_{\mathrm{base}}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{a}}_{\mathit{j}}{\mathit{T}}_{\mathit{j}}\mathrm{\left(}\mathit{u}\mathrm{\right)}$where

$\mathit{u}\mathrm{=}\mathrm{2}\left(\frac{\mathit{v}\mathrm{-}{\mathit{v}}_{\mathrm{lo}}}{{\mathit{v}}_{\mathrm{hi}}\mathrm{-}{\mathit{v}}_{\mathrm{lo}}}\right)\mathrm{-}\mathrm{1.}$One then has to decide at what value of n to truncate the sequence. The Bayesian methodology offers a way to decide this.

The posterior function p(q | y) giving the probability that the data y is fitted by model parameter values q is given in Eqs. (5) and (7). So far we have not concerned ourselves with the evidence E, the denominator of Eq. (5), which is the integral over the parameter space of the numerator. If we change the number of parameters, however, the value of E will vary, and can be used to decide the optimum number.

In the present case the evidence integral has 6 + n dimensions, and considering for example Fig. B.5, we can see that the integrand, at least in the directions of the 6 profile-only parameters, is unlikely to be well behaved. However, the model is linear in the Chebyshev functions, which means that in these n directions, given Gaussian-distributed noise as we have, the likelihood function has also a Gaussian shape in these dimensions. In addition, the correlation between the n Chebyshev parameters and the 6 profile parameters may be expected to be small, particularly if we include rather more baseline than there are channels under the line profile. Therefore we can approximate the posterior probability function for the whole model by a product between the profile and baseline contributions: $\begin{array}{ccc}\mathit{p}\mathrm{\left(}{q}\mathrm{\right|}{y}\mathit{,n}\mathrm{)}& \mathrm{=}& \mathit{p}\mathrm{(}{{q}}_{\mathrm{profile}}\mathrm{+}{{q}}^{\mathrm{\prime }}\mathrm{|}{y}\mathit{,n}\mathrm{)}\\ & \mathrm{~}& \mathit{p}\mathrm{\left(}{{q}}_{\mathrm{profile}}\mathrm{\right|}{y}\mathrm{)}\mathrm{\times }\mathrm{[}\mathit{p}\mathrm{\left(}{{q}}^{\mathrm{\prime }}\mathrm{\right|}\mathit{n}\mathrm{\left)}\hspace{0.17em}\mathit{p}\mathrm{\right(}{y}\mathrm{\left|}{{q}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\mathrm{]}\mathit{.}\end{array}$(B.4)Here q′ represents the n baseline parameters, and the dependence of the probabilities on n has been made explicit.

Deciding the priors p(q′ | n) presents its usual difficulties. Since the Chebyshev polynomials have a significant degree of orthogonality, we can approximate p(q′ | n) by a product of n separate priors. We choose the same very broad prior for each Chebyshev term. “Broad” of course implies slowly-varying compared to the likelihood function; thus each baseline prior may be replaced by its value at the posterior maximum, which we designate P_base. Since all n of them are separate and equal, $\mathit{p}\mathrm{\left(}{{q}}^{\mathrm{\prime }}\mathrm{\right|}\mathit{n}\mathrm{)}\mathrm{=}{\mathit{P}}_{\mathrm{base}}^{\mathit{n}}$.

The Gaussian likelihood for the n Chebyshev parameters has, following Eq. (7), the form $\mathit{p}\mathrm{\left(}{y}\mathrm{\right|}{{q}}^{\mathrm{\prime }}\mathit{,n}\mathrm{)}\mathrm{=}\mathit{p}\mathrm{(}{y}\mathrm{\left|}{{q}}_{\mathrm{0}}^{\mathrm{\prime }}\mathit{,n}\mathrm{\right)}\exp \left(\mathrm{-}\mathrm{[}{{q}}^{\mathrm{\prime }}\mathrm{-}{{q}}_{\mathrm{0}}^{\mathrm{\prime }}{\mathrm{]}}^{\mathrm{T}}\frac{{{C}}^{\mathrm{\prime }\mathrm{-}\mathrm{1}}}{\mathrm{2}}\mathrm{[}{{q}}^{\mathrm{\prime }}\mathrm{-}{{q}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{]}\right)\mathit{.}$(B.5)Here ${{q}}_{\mathrm{0}}^{\mathrm{\prime }}$ represents the values of the baseline parameters at the posterior maximum, and C′ is the covariance matrix of the posterior with respect to those parameters. The T superscript indicates transpose.

From Eqs. (6), (B.4) and (B.5) the evidence evaluates to

${\mathit{E}}_{\mathit{n}}\mathrm{=}{\mathit{\kappa }}_{\mathit{n}}{\mathit{P}}_{\mathrm{base}}^{\mathit{n}}\sqrt{\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}}\det \mathrm{\right(}{{C}}^{\mathrm{\prime }}\mathrm{)}}$where

${\mathit{\kappa }}_{\mathit{n}}\mathrm{=}\mathit{p}\mathrm{\left(}{y}\mathrm{\right|}{{q}}_{\mathrm{0}}^{\mathrm{\prime }}\mathit{,n}\mathrm{)}\mathrm{\int }\mathit{d}{{q}}_{\mathrm{profile}}\hspace{0.17em}\mathit{p}\mathrm{(}{{q}}_{\mathrm{profile}}\mathrm{)}\hspace{0.17em}\mathit{p}\mathrm{(}{y}\mathrm{\left|}{{q}}_{\mathrm{profile}}\mathrm{\right)}\mathit{.}$Because of the weak coupling between the Chebyshev contribution and the profile model proper, we assume that κ_n can be decomposed into

${\mathit{\kappa }}_{\mathit{n}}\mathrm{~}\mathit{Kp}\mathrm{\left(}{{q}}_{\mathrm{0}}\mathrm{\right|}{y}\mathit{,n}\mathrm{)}$where K is a constant and p(q₀ | y,n) is the maximum value of the posterior, which can be measured at the end of the fit process. We cannot easily calculate K, but since it will be the same for any n, for purposes of comparing different n values we can simply evaluate ${\mathit{E}}_{\mathit{n}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{E}}_{\mathit{n}}}{\mathit{K}}\mathrm{=}\mathit{p}\mathrm{\left(}{{q}}_{\mathrm{0}}\mathrm{\right|}{y}\mathit{,n}\mathrm{)}{\mathit{P}}_{\mathrm{base}}^{\mathit{n}}\sqrt{\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{n}}\det \mathrm{\right(}{{C}}^{\mathrm{\prime }}\mathrm{)}}\mathit{.}$(B.6)Intuitively, as n increases, we would expect the fit to get better and better, reflected in an increased value of the posterior maximum, and thus of κ_n in Eq. (B.6); but since P_base for a broad baseline prior is significantly less than 1, ${\mathit{P}}_{\mathrm{base}}^{\mathit{n}}$ will of course become geometrically smaller, and will eventually offset the rise in κ_n. This is what is known as the “Ockham’s Razor” feature of the Bayesian methodology, in that it weights against a model with too many parameters (D’Agostini 2003; Gregory & Loredo 1992).

Shown in Fig. B.7 is the value of ${\mathit{E}}_{\mathit{n}}^{\mathrm{\prime }}$ for different values of n. Two Gaussian priors for the Chebyshevs are tested, one having a width equal to 10 times the noise standard deviation (square points), the other being 10 times wider again (crosses). It can be seen that the width of the prior has small effect, provided only it is wide. In both cases the evidence shows a peak at n = 6. This was therefore the value chosen for the baseline fits.

Appendix B.7: Baseline priors for the EBHIS profiles

Many of the EBHIS cubes show significant baseline ripple. This must be fitted to obtain the most accurate estimates of the total flux (and to a lesser extent, the other parameters) of spectral lines. Of course the baseline variations continue smoothly through the spectral line itself, but if one includes the spectral line in an unconstrained baseline fit, one risks that the line itself will distort the baseline estimation. Traditionally this is accommodated by manually excluding the line channels from a separate baseline fit. Formally speaking however it is better to fit the same model to every channel, but to constrain the baseline parameters by priors. This has been our approach.

Gaussian priors were chosen for each of the Chebyshev amplitudes in the baseline model. The centre and spread for each were obtained in the following way. Firstly, the baseline model was fitted over the same channel range to the spectrum at every spatial pixel of the data cube, excluding only those pixels either where a significant source contribution was obvious, or at the edges of the cube where the noise was higher than average. Then, for each Chebyshev order, the centre of the prior was taken from the mean of the whole-cube fitted amplitudes, multiplied by the number of spatial pixels summed over the source at issue; and the spread (sigma value) of the prior was assigned from the standard deviation of these amplitudes, suitably corrected for the existing correlation between spatial pixels which is caused by the breadth of the telescope beam.

Appendix C: Graphs of the fits to all spectra

All Tables

All Figures