The analysis of indexed astronomical time series – VI. Covariances of amplitude ratios and phase differences estimated from multicolour photometry of a pulsating star

One of the tools used to identify the pulsation modes of stars is a comparison of the amplitudes and phases as observed photometrically at different wavelengths. Proper application of the method requires that the errors on the measured quantities, and the correlations between them, be known (or at least estimated). It is assumed that contemporaneous measurements of the light intensity of a pulsating star are obtained in several wavebands. It is also assumed that the measurements are regularly spaced in time, although there may be missing observations. The amplitude and phase of the pulsation are estimated separately for each of the wavebands, and amplitude ratios and phase differences are calculated. A general scheme for estimating the covariance matrix of the amplitude ratios and phase differences is described. The first step is to fit a time series to the residuals after pre-whitening the observations by the best-fitting sinusoid. The residuals are then cross-correlated to study the interdependence between the errors in the different wavebands. Once the multivariate time-series structure can be modelled, the covariance matrix can be found by bootstrapping. An illustrative application is described in detail.


I N T R O D U C T I O N
One of the methods commonly used to identify the oscillation modes in pulsating stars is the use of photometric observations obtained in different wavebands. Pulsation theory predicts certain relationships between the amplitudes (or phases) of the pulsation light curves, as seen in the different wavebands. The theory often contains undetermined parameters, and the predictions are commonly presented as regions in two-dimensional plots of, e.g., the ratio of the pulsation amplitudes in two different wavelengths, versus the phase difference between the light curves seen through the two different filters. The observed amplitude ratio and phase difference then provide a single point in the diagram, which hopefully lies in a region in the plane that corresponds to some particular pulsation mode.
In principle it is obviously advantageous to obtain light curves in as many wavebands as possible. There are, however, potential problems with applying the method described above if observations are made through three or more filters: for three filters, there are nine possible combinations of amplitude ratios and phase differences; for four and five filters the numbers grow to 36 and 100. Not only is it not feasible to study so many plots, but it is also unclear how conflicting results ought to be reconciled. Clearly what is needed is a way of combining the many parameters in some way so as to obtain a single solution that is optimal in some sense. Koen et al. (1999;hereafter KRWM99) proposed minimizing the sums of squares of the differences between the observed and theoretical quantities. [The same fundamental approach was taken by Balona & Evers (1999), although the details of their method differ considerably from those presented below.] The minimization is with respect to the unknown pulsation mode of the star, as well as any other unknown parameters (e.g. those measuring the departure of the oscillations from adiabaticity). Clearly such a minimization needs to incorporate knowledge of the relative accuracies of the observed quantities, as well as correlations amongst these. Such knowledge is summarized in the covariance matrix S of the observables. This paper is devoted to a method for the estimation of S.
If the standard errors of the estimated quantities (and their correlations) were known, the problem would be relatively simple (KRWM99). One of the complicating factors is the estimation of the standard errors of estimated amplitudes and phases, in situations where correlated noise is superposed on the signal of interest (Schwarzenberg-Czerny 1991; Montgomery & O'Donoghue 1999;Koen 1999, hereafter Paper V). In the present context there is the added complication of correlations between noise processes measured through different filters.
The complexity of the problem does not encourage an analytical approach; therefore the method proposed is bootstrapping, also known as computer resampling (e.g. Efron & Tibshirani 1993). The various steps in implementing the bootstrap are as follows.
(1) The observations are pre-whitened by the sinusoid(s) arising from the pulsations. This leaves a noise process for each filter.
(2) A time-series model is fitted to the noise process associated with each of the filters. The model residuals, or`innovations', should be free of autocorrelation.
(3) The innovation sequences associated with the different wavebands are cross-correlated, and the results are used to construct a covariance matrix for all the innovations.
(4) The innovations are transformed into a new set of residuals, which are free of both autocorrelation and cross-correlations. Since these new residuals do not have any sequential coherence at all, any permutation of them is statistically completely equivalent.
(5) The final residuals are used to re-construct sets of observations that are statistically equivalent to the originals. The parameters of interest (amplitude ratios, phase differences) are estimated for the synthetic data sets.
( 6) Step 5 is repeated a large number of times, and the results are used to establish the distributional properties of the parameters of interest.
The general strategy is described in some detail in Section 2, followed by a complete analysis of a`real life' example. Section 4 contains a few remarks.

M E T H O D O L O G Y
The situation considered is as follows: contemporaneous observations in several wavebands are obtained of a star that is varying sinusoidally with a known angular frequency v. (Extension to the case of several different periodicities is straightforward, but the notation more complicated.) The observations are described by where k indexes the filter and j is the time of observation. The notations C(k) and f(k) are used for the best-fitting (in the leastsquares sense) amplitude and phase of the variations in waveband k. The residuals left after subtracting the best-fitting sinusoid from the observations y j (k) are denoted by e j (k). The observables of direct interest are the amplitude ratios u1 C2aC1Y u2 C3aC2Y ¼Y uK 2 1 CKaCK 2 1Y 2 and the phase differences It is convenient to stack the u( j ) and v( j ) in a single vector z, with its M 2K 2 2 components given by Of interest here is the covariance matrix S of the vector z. Two important remarks need to be made: first, if the K noise sequences {e 1 kY e 2 kY ¼Y e N k} are white noise processes which are also uncorrelated with one another, calculation of S is fairly straightforward. The usual regression covariance matrix of the C(k) and f(k) is reliable in this case, and may be used to calculate a good approximation for S (KRWM99). In what follows it is assumed that the K individual noise processes may be autocorrelated, as well as cross-correlated with one another. This invalidates use of the ordinary regression covariance results. The second point concerns the assumption that the observations are made at regularly spaced discrete time intervals. This is obviously quite restrictive: it typically applies only to multichannel high-speed photometry. Another observational configuration that falls within this ambit is that of continuously cycling through a filter sequence at a set speed, as when doing CCD photometry. In such cases the fundamental time unit is that required to complete one filter sequence.
We propose to calculate S either by simulation, or by a technique called`bootstrapping' (see below, and e.g. Efron & Tibshirani 1993). Both methods require statistical models for the individual noise processes {e j (k)}, as well as for the interrelationships amongst them. The next several paragraphs are devoted to a description of how to accomplish this.
A simplifying, but non-essential, assumption is that the sequences of residuals {e j (k)} are statistically stationary, i.e. have properties (such as the mean value and variance) that do not change with time. (The principles involved do not change if the series are non-stationary, but additional treatment of the sequences may be required. The example observations analysed in Section 3 are in fact non-stationary in variance, but the data can be divided into segments that are stationary.) This innocuous assumption has the important consequence that the sequences of residuals may be modelled as autoregressive moving average (ARMA) processes ± see Koen & Lombard (1993, Paper I in this series). The most general form of an ARMA process is where the sequences {h 1 kY h 2 kY ¼Y h N k} are white noise. The notation`ARMA(p, q)' is used, indicating that the process has autoregressive and moving average orders of respectively p and q.
In practice, the great majority of processes are of forms with 0 # p # 2Y 0 # q # 1X It would be out of place to describe the identification of the orders p and q in detail here; the interested reader is referred to the classical text by Box & Jenkins (1976), or more modern treatments such as that by Brockwell & Davis (1991). The actual fitting of models is also dealt with in these texts, as well as in almost any other modern book on time-series analysis. A large number of commercial software packages capable of full ARMA analyses are available, and the accompanying manuals often contain brief descriptions of the procedures necessary to identify, estimate and test the optimal model. A possible complication is that there may be`missing values', i.e. observations are not available at all N time points. (For example, data streams might be interrupted in order to measure sky background brightness, or because of poor weather or instrument failures.) Commercial software is rarely equipped to deal with such data: often missing values are simply filled in with the overall mean values. This may not be serious if time series are almost complete, but may pose a substantial problem if there are many missing observations. Fortunately, there are statistical methods available for ARMA analysis of large data sets from which substantial numbers of values are missing: Paper I contains a description of a technique due to Jones (1980). Jones's (1980) algorithm is efficient in the sense that it based on maximizing the likelihood of the data. It avoids having to deal with the full covariance matrices of the data by drawing on a Kalman filter algorithm. Again, it would take us too far afield to deal with this in detail; the interested reader is referred to the papers mentioned.
Once the optimal ARMA model has been selected, it is necessary to test whether the model residuals are indeed uncorrelated. This is easily done by calculating the appropriate autocorrelation function, and checking that it is consistent with the residuals (estimated h) being white noise. If the residuals are correlated, a different model ± often suggested by the pattern of the residual autocorrelations ± is fitted. The process is continued until satisfactory models have been fitted to the K individual series of observations. Note that missing observations require that some care be taken with the calculation of the autocorrelation function: equation (24) in Paper V is an appropriate formula.
The next step in the analysis is to enquire whether the individual series are related. It might be thought that this could have been done very simply at the outset by calculation of cross-correlation functions between the different sets of residuals defined in (1). This happens to be fraught with danger: cross-correlation of two completely unrelated series, which are themselves autocorrelated, can often lead to spurious evidence of strong correlation between them (see e.g. Jenkins 1979;Koen 1994). What is required is first to remove the autocorrelation from the individual series. Once appropriate ARMA models have been found for the individual series, this is trivial: the ARMA model is subtracted from the original sequence of observations (a procedure referred to as`prewhitening'), leaving the correlation-free set of {h j (k)}. Crosscorrelating the {h (k)} for different values of k will then allow conclusions to be drawn about the relatedness of the error processes associated with the different filters.
Only the case of limited-span cross-correlation is considered in this paper: by this is meant that the different h(k) are only correlated if they are close to being contemporaneous. Symbolically, for finite L. Typically L will have values of 1 or 2. It is useful to define the vectors We also define Note that A`is independent of the index j, by the assumption of statistical stationarity. Also, A 2`; A t`X Using the notation above, the covariance matrix of all the h j (k) can be written quite concisely; for the case L 1 it is The form for larger L has a similar band diagonal structure. The matrix S is a covariance matrix, and hence it is symmetric and positive definite; this implies that it can be factorized as where L is a lower triangular matrix (Cholesky decomposition ± see e.g. Searle 1971). Define then, from (10), where I is the identity matrix. The last equation implies that the components of the vector Z are uncorrelated, and constitute the basic building blocks from which the observed noise time series can be reconstructed.
To summarize: the cross-correlation function is used to establish the value L of the maximum lag, and to estimate the entries in the matrices A`; the matrix S is constructed from the A( equation 9); the matrix L follows from (10). The vector Z follows from (11); its components are fully uncorrelated.
It remains to discuss the generation of the synthetic data sets to be used in the bootstrapping procedure. The first step is to generate sets of vectors that are statistically equivalent to Z. This is easy, as the components of Z are uncorrelated; the equivalents Z H are formed by randomly drawing (with replacement) the necessary number of components from Z.  (7), we then have a set of N K-vectors w j H , from which the full K bootstrap noise processes {e j H (k)} may be created ± see (5). The simulated observations follow immediately from equation (1).
As mentioned in the Introduction, the procedure described in the previous paragraph is repeated a very large number of times, and the amplitudes C(k) and phases f (k) estimated for each bootstrap sample. The covariances of the amplitude ratios and phase differences are then easily calculated from the collection of simulation results.

AN EXAMPLE
KRWM99 published the results of UBVRI photometry of the d Scuti star HD 95321. The star is apparently monoperiodic, with P 0X213 658 dX The observations were obtained by CCD photometry, continuously cycling through the set of filters. One cycle of UBVRI observations lasted 248.44 s, and this is used as the unit of time in discretizing the time variable. The observing run spanned six nights, one of which was not usable; the individual observing sessions ranged from 4.5 to 5.3 h. The number of useful observations ranged from 225 in the I band to 304 in the U band. An example of the noise processes in the different wavebands is shown in Fig. 1: these are the residuals after subtracting the best-fitting sinusoid with P 0X213 656 d from each of the sets of differential photometry obtained on a particular night.
There are of course large time lapses between the nightly observing runs, and these gaps are not integer multiples of 248.44 s. There is thus some ambiguity in the`cycle count' during the gaps. This is of little consequence in dealing with the identification and simulation of the noise processes, as these are essentially independent from night to night. Fitting a short-period sinusoid over the entire set of times is a different matter, but need not be of great concern: we simply define a new time series on the estimated integer times, with a period equal to that of HD 95321. Although this is not an exact representation of the observed situation, we are only interested in the variability of the amplitude and phase estimates (i.e. relative, rather than absolute, values): the newly defined time series will clearly be an adequate proxy for the real thing as far as this aim is concerned.
We proceed with an analysis of the residuals left after prewhitening the observations by the best-fitting sinusoid with the appropriate period. The first task is to fit models of the form (5) to these data. If the residuals from each night are modelled individually, a total of 25 models are required (five nights of observations, for each of the five filters used). There are several reasons why, if it is at all viable, it is desirable to fit a single model to all the residuals for a given filter. These include simplicity, reduction in the modelling effort, and a gain in accuracy through fitting to a larger data set. Note that this assumption is also implicit in the description, and notation, used in Section 2 above.
The statistical structure of the noise in Fig. 1 is typical of all the residuals. Clues such as the obviously long correlation memory suggests an ARMA(1,0) model as a first attempt. Examination of the set of innovations {h j (k)} [see (5)] showed this model to be inadequate: innovations were strongly autocorrelated at a lag of one time unit. Such a limited-length autocorrelation structure is typical for an ARMA(0,1) time series, which, combined with the original model, implies that an ARMA(1, 1) model may fit the noise well. This turned out to be so: the autocorrelation functions of the innovations of ARMA(1,1) model fits were found to be consistent with white noise, implying that the models gave good representations of the entire correlation structure in the K individual noise processes {e j (k)}. It is reassuring that an ARMA(1, 1) is the optimal model (according to the Bayes information criterion ± see Paper I) for the data for all of the filters: it suggests that this model may be describing common underlying noise processes. (See also Paper V, in which similar conclusions were reached for a large number of photometric noise processes.) The models are summarized in Table 1.
An examination of plots of the innovations against time shows that the scale of the scatter is different during different nights. Since the innovations are not necessarily Gaussian (see below), care has to be taken in applying standard tests for differences in the variances over different nights. For example, the standard Bartlett's test (Bartlett 1937) is known to be sensitive to deviations from normality, and should not be used. Fortunately, the innovations are uncorrelated, and hence use can be made of distribution-free permutation tests. We use the simple statistic (Hartley 1950) Hk max j s 2 j kamin j s 2 j kY 13 where j indexes the different nights, and s 2 j k is the variance of the innovations during night j in waveband k. The statistic H(k) is first calculated for the original set of innovations. The innovations  Table 1. The results of fitting ARMA(1,1) models to the data for each of the filters. All nightly data for a given filter have been combined. The standard deviation s h of the innovations is given in millimagnitudes. The first three columns are the initial estimates; the other columns contain estimated values after correcting the data for heteroscedasticity. The last column shows the minimum standard deviation of the innovations, followed by scalefactors for each of the five nights. from all the nights are then put together in a single collection, which is shuffled into random order, and partitioned into groups of the same sizes as the data sets for the individual nights. Finally, the Hartley (1950) statistic is calculated for the permutation. The full set of innovations is again shuffled, partitioned, and H(k) calculated; this procedure is repeated a large number of times (typically a few thousand times), and the observed H statistic is ranked in the distribution of permutation Hs in order to obtain a significance level. Based on the results of 5000 permutations per set of innovations, it was concluded that all H(k) statistics except that for the U-band data were highly significant. This so-called heteroscedasticity of the data may be modelled in a very simple way by assuming that s h is constant over each night, i.e. behaves as a step function S Sj of time. The modification required to the algorithm of Jones (1980) is surprisingly modest: the variance s 2 in his equation (3.2) needs to be replaced by s 2 St 1X (A detailed discussion of Kalman filtering of heteroscedastic data is given in Koen 1996, Paper IV.) Comparison of the right-hand side of Table 1 with columns (2) and (3) makes it clear that a 1 and b 1 are affected very little when provision is made in the estimation process for the heteroscedasticity of the data. None the less, it is noteworthy that the estimated innovations may be significantly affected.

Filter
Several of the innovations are larger than 4 standard deviations in absolute value; this raises the issue of testing for outliers. A sourcebook of many outlier tests is that by Barnett & Lewis (1994); we use the simple two-sided test statistic where the (standard) notation h (j) indicates the jth largest innovation, and s(k) is the standard deviation of the innovations associated with waveband k. Unfortunately, the significance level tables in Barnett & Lewis (1994) only extend as far as sample sizes of 120, so that it is necessary to determine values for our data by simulation.
It was found that the B, V and R data sets each had one outlying innovation (T statistics of 4.87, 4.78 and 4.18; significance levels of 0.0004, 0.0004 and 0.005); after removal of these, all T-values have significance levels larger than 0.08.
An assumption of the outlier test is that the innovations have a Gaussian distribution. Since this is also a point of general interest, formal significance tests for deviation from normality were carried out. Three so-called`empirical distribution function' statistics were calculated for each of the five sets of innovations: the Anderson±Darling A 2 , Crame Âr±Von Mises W 2 and Watson U 2 statistics [see D 'Agostino & Stephens (1986) for definitions, significance level tables, etc.]. The V-band innovations deviate from normality more significantly than the 5 per cent level; the I-band innovations are marginally non-normal (U 2 is significant at the 5 per cent level; the other two statistics are non-significant); and the innovations associated with the other three wavebands have distributions consistent with Gaussians. The results imply that some caution needs to be exercised when basing deductions on the assumption of normality, at least for the V-band data.
A point regarding the tests for outliers and for normality worth remarking on is that the mean values of the sets of innovations were assumed known. This has a marked influence on the test results: the tests are a lot less stringent than they would be if sample means were used. However, as the innovations are effectively regression residuals, the assumption that the mean values are equal to zero is certainly plausible.
The results of cross-correlating the innovations associated with the various filters are given in Table 2. The standard errors of the cross-correlations are estimated as M 21/2 , where M is the number of terms summed in calculating the cross-correlation. All correlations larger than two standard errors, and with lags j`j # 5Y are listed. As expected, the strongest correlations occur at lags of zero and one time units. There are two reasons for doubting correlations at lags of 2 units or more. First, it seems implausible on physical grounds that there could be strong correlations at a larger lag, without significant correlation also occurring at smaller lags. Secondly, with only two exceptions, the correlations at lags larger than one unit are only marginally significant. In what follows, only correlations with j`j # 1 are considered meaningful, i.e. it is assumed that L 1 in (6).
The two non-zero matrices A`(with magnitudes in mmag) are There are a number of important practical issues in the calculation of the matrix L in (10). First, the data are incomplete, since some of the components of w i are sometimes unavailable. In Table 2. Significant cross-correlations between h j (k 1 ) and h j+`( k 2 ). Since covh j k 1 Y h j` k 2 ; covh j k 2 Y h j2` k 1 Y not all of the results are shown. such cases, the relevant rows and columns of A 0 and A 1 are deleted. There are also instances in which entire vectors w i are missing (as, for example, in the gaps between observing runs). In these cases successive w j are not correlated, and hence the matrix A 1 must be replaced by the null matrix. The second issue of some importance is the size of the matrix S; in the present instance the total number of elements in W is 1347, giving more than 1.8 million entries in S. Of course, only half this number needs to be stored, as S is symmetrical, but the computer memory requirement is still enormous. One possible way of dealing with the problem is to make use of the sparsity of S and L, and to keep careful track of the matrix positions of the non-zero elements. There is also a much simpler approach which can be used in the present case: the longest uninterrupted stretch of data is about 40 time units long, and each such block of residuals is uncorrelated, or very weakly correlated, with the rest. This means that S can be treated as consisting of a number of relatively small uncorrelated blocks, each of which can be independently Cholesky-decomposed as required in (10). For each independent block, the relevant part of the vector Z follows from (11). A slightly different strategy is followed in implementing (12), as the ultimate aim is to simulate the noise processes e j (k) described by (5). This requires having values of e H j2i (k), i.e. noise values from time periods t j 2 i preceding the epoch t jX We deal with this by calculating matrices S H as in (9), but (i) padding each of the independent blocks of which S is composed by an extra 10 rows of the form in (9), and (ii) ignoring the possibility of missing observations. Thus, if the block of observations is of length n, we simulate the covariance matrix S H of n 10 observations with the required covariances. The corresponding W H follows from (12), and the required {e H j (k)} from (5). Finally, the e H j (k) corresponding to unobserved e j (k) are ignored. Results based on 5000 bootstrap replications of the original data are now discussed. Table 3 compares the estimates of the standard errors from the ordinary linear least-squares (LLS) theory and those obtained from bootstrapping, for the amplitudes and phases of the sinusoid observed in each of the five wavebands. It appears that the LLS amplitude standard errors underestimate the true errors by a factor of about 2, while the LLS phase standard errors underestimate by a factor of about 3. Defining the vectors of estimated amplitudes and phases, the correlation matrices are also of interest. Some of these correlations are clearly nonnegligible, and will play a role when comparing theoretical and observed amplitude ratios and phase differences. Correlations between the amplitude and phase estimates are quite small: the largest absolute value is less than 0.1, implying that the quantities in equations (2) and (3) are essentially decoupled from one another.
The final results are given in Table 4, and in the correlation matrices (20) and (21). Perhaps the most notable result in Table 4 is that, with the now more dependable errors, none of the phase differences is significantly non-zero. In view of the results in Table  3, it is not surprising that the LLS errors quoted in Table 4 underestimate the true errors by factors of 2 to 3.
Defining vectors U and V with components u(i) and v( j) as in equations (2) and (3), the amplitude ratio and phase difference correlation matrices are As expected, the correlations between the u(i) and v( j) are small in absolute value, the largest being 0.06.

A F E W R E M A R K S
The casual reader may be left with a sense of disappointment that the author has not provided a simple, clearcut, universally applicable solution to the problem discussed above. It may therefore be useful to state what the author's intentions were: first, to attempt to demonstrate that, at least in principle, the required error covariance matrices can be calculated; secondly, to map out a general line of attack on the problem, which should work in most cases; and thirdly, to develop specific results which should provide some sense of the difference in the final answers between a full solution and the simple-minded approach used by KRWM99. The scheme used above was to some extent dictated by the availability of software, rather than sound theoretical considerations: it is certainly desirable to fit a complete multivariate timeseries model to all the data simultaneously, rather than use a piecemeal strategy such as done here. It may be surmised that a VARMA(1, 1) model [i.e. a multivariate ARMA(1, 1) ± see e.g. Lu Ètkepohl 1993] would furnish a good description of the data. However, development of the software necessary to fit such models to data with missing values would be non-trivial.
A considerable part of the labour involved in calculating (20), (21) and the bootstrap standard errors in Table 4 stems from the requirement that the statistical relationships between the innovation processes associated with the different filters be faithfully reproduced in the bootstrapping machinery. If the processes in the different wavebands were unrelated, A 0 would be a diagonal matrix, A j would be null for j ± 0Y and S would be diagonal; the triangular matrix L would reduce to diagonal form, with elements equal to the square root of the elements of S. Clearly, the computer memory requirements would be very modest indeed, and the programming considerably simplified. It therefore seemed worthwhile to re-calculate standard errors and correlation matrices of the amplitude ratios and phase differences, ignoring the interrelated nature of the five series of observations. Interestingly, the standard errors are increased by 5±18 per cent, a typical figure being 15 per cent. Comparison of (22) and (23) with (20) and (21) shows that some of the correlations are quite strongly affected, so that neglecting the inter-dependence of the observed series may be risky: We close with a comment on the impact of the above results on the identification of the pulsation mode of HD 95321. KRWM99 proposed a weighted least-squares method for the identification of the degree`of the pulsation: the quantity to be minimized is SS z 2 z p S 21 z 2 z p t Y 24 where z and S are defined in Section 2 above, and z p is the theoretically predicted value of the vector z. Table 5 compares the results for covariance matrices S based on ordinary least-squares (KRWM99) and the bootstrapping technique of this paper. The ordering of the likelihood of the modes in terms of SS is the same in the two cases, but the contrasts are much reduced in the case of S based on bootstrapping.
The last observation has an important consequence for meaningful mode discrimination. The Gaussian log likelihood of z is L 24 ln 2p 2 1 2 ln jSj 2 1 2 z 2 z p S 21 z 2 z p t 24 ln 2p 2 1 2 ln jSj 2 1 2 SSX 25 This means that the likelihood ratio statistics for comparing the results for different`reduce to l`1Y`2 2ln L`1 2 ln`2 SS`2 2 SS`1X 26 In order to gain some insight into whether the results for the different modes in Table 5 differ meaningfully, we use the fact that the l(`1,`2) with`1 corresponding to the optimal solution is asymptotically chi-squared-distributed with three degrees of freedom. (The three arise from the fact that not only`but also the two non-adiabaticity parameter R and C are estimated for each  Table 4. The amplitude ratios u and phase differences v for the different wavebands in which HD 95321 was observed, as given in KRWM99. The indices of the amplitudes and phases in (2) and (3)  solution.) The interested reader is referred to any general introductory text on statistics (e.g. Mood, Graybill & Boes 1974) for more information.
The likelihood ratio statistics are also shown in Table 5. In the case of the ordinary least-squares entries, solutions for`$ 3 can be rejected at a confidence level far better than 1 per cent; for the bootstrapping entries the worst and best solutions differ at only the 20 per cent level. The implication is that, with realistic error estimates, clearcut identification of the pulsation mode of HD 95321 from the KRWM99 multicolour photometry may be impossible. The KRWM99 data are rather more extensive, in terms of both wavelength coverage and time coverage, than those that have often used for mode identification. A repetition of the exercise of this paper for other stars is called for, in order to establish whether mode identification from multicolour photometry alone is currently viable.