Mode width fitting with a simple bayesian approach.
Application to CoRoT targets HD 181420 and HD 49933††thanks: The CoRoT space mission, launched on
2006 December 27, was developed and is operated by the CNES, with
participation of the Science Programs of ESA, ESA’s RSSD, Austria,
Belgium, Brazil, Germany and Spain.
Key Words.:statistics – detection – data analysis – asteroseismology
Aims:We investigate the asteroseismology of two solar-like targets as observed with the CoRoT satellite, with particular attention paid to the mode fitting. HD 181420 and HD 49933 are typical CoRoT solar-like targets (156 and 60-day runs). The low signal-to-noise ratio (SNR) of about prevents us from unambiguously identifying the individual oscillation modes. In particular, convergence problems appear at the edges of the oscillation spectrum.
Methods:We apply a Bayesian approach to the analysis of these data. We compare the global fitting of the power spectra of this time series, obtained by the classical maximum likelihood (MLE) and the maximum a posteriori (MAP) estimators.
Results:We examine the impact of the choice of the priors upon the fitted parameters. We also propose to reduce the number of free parameters in the fitting, by replacing the individual estimate of mode height associated with each overtone by a continuous function of frequency (Gaussian profile).
Conclusions:The MAP appears as a powerful tool to constrain the global fits, but it must be used carefully and only with reliable priors. The mode width of the stars increases with the frequency over all the oscillation spectrum.
The CNES CoRoT satellite has been orbiting the Earth at an altitude of 800 km, in a polar orbit, since late December 2006. Its scientific goals are divided into two themes: asteroseismology and the detection of transiting exo-planets. In this paper, we focus on the analysis of seismological data, in particular on the analysis of the power spectra of the solar-like targets. The main parameters extracted from observations, in order to constrain stellar models, are the global mode frequencies, heights, widths, and parameters such as the rotational frequency splitting and the star inclination with respect to the line of sight. These parameters are extracted by fitting the power spectrum with a mathematical model, which consists of a sum of Lorentzian functions, justified by assuming that the oscillations are stochastically excited (Duvall & Harvey 1986).
Such a technique has been commonly used with helioseismic measurements such as GONG, BiSON, GOLF, LOI, MDI (for a review of these instruments see Howe 2009). However, thanks to the excellent signal to noise ratio (SNR) of solar data (40 with LOI, 300 with GOLF), the mode fitting could be done individually for each pair of modes (odd or even). Moreover, in the Sun’s case, the inclination is well known as is the rotational splitting. The difficulty of fitting a stellar power spectrum is firstly that the SNR is less than 10 (e.g. Appourchaux et al. 2008) and secondly that both inclination and splitting are generally unknown, even if they have been roughly estimated with the help of complementary radial velocity measurements, for which the estimates depend strongly on the stellar model (size, mass).
Therefore, in order to extract the maximum information from the solar-like stellar spectra obtained with CoRoT, it appeared unavoidable to globally fit the power spectra, otherwise, the global parameters such as the splitting, the inclination and the stellar background could not be properly fitted (Appourchaux et al., 2008). Consequently, the number of parameters to be simultaneously estimated increases substantially. As an example, in the case of the solar-like target HD 49933, where 14 overtones have been taken into account, a set of 75 free parameters was needed (see Appourchaux et al. 2008). The fitted values are obtained with the Maximum Likelihood Estimator (MLE) (see Sect. 1). However, as already attempted with several CoRoT solar-like targets (HD 49933: Appourchaux et al. 2008, HD 181420: Barban et al. 2009, HD 181906: Garcia & et al. 2009), the MLE struggles to fit the power spectra when the SNR is lower than a given threshold (example in Sect. 1); the risk is to get a biased estimate of global parameters, because of the bad fitting of the “external modes” (at the edges of the power spectrum).
In this context, Appourchaux (2008) suggested adopting a Bayesian approach in order to make the most of all of the available a priori information we have on the star. Following those suggestions, Benomar et al. (2009) applied a Bayesian approach to fit the first run of HD 49933. The estimate of the parameters was determined by scanning the probability density function (PDF) associated with each of the 75 parameters with the numerical “Monte-Carlo Markov Chain” (MCMC) method. The result of this analysis gave a significantly more reliable result, particularly in preventing some demonstrably erroneous estimates (e.g. “Dirac-like” modes, see Sect. 1). The only real drawback of such a method is the long computational time (several days for a global fitting of HD 49933, Benomar, private communication).
In this paper, we propose an a minima Bayesian approach, whose objective is to prevent erroneous solutions and to quickly extract the main trend of the stellar parameters. In that context, we use the Maximum A Posteriori (MAP) approach, which roughly consists of adding prior information to the MLE, then maximizing it in the same way. In other words, it is a regularized MLE. The main difference with respect to the MCMC method is that we do not explore the posterior PDF but simply determine the local maxima of the posterior PDF. The second purpose of the paper, which can be seen as a particular case of the MAP approach, is to reduce the total number of parameters to be fitted. The work is organized as follows: Sect. 2 contains a summary of the MLE approach and gives an example of the typical problems that may appear; Sect 3 presents the Bayesian approach and discusses the risks we run by using it improperly; Sect. 4 is dedicated to the reduction of the number of free parameters and its application to the GOLF solar data; finally, Sect. 5 shows the application to 2 CoRoT solar-like targets: HD 181420 and HD 49933. We particularly pay attention to the resulting mode width dependence on the frequency.
2 Inability of MLE in the case of noisy spectra
2.1 Global power spectrum fitting
With the assumption of stochastic excitation of the modes, the power spectrum of the eigenmodes is distributed around a mean profile following a with 2 degrees of freedom (d.o.f). It can be proved by considering the photometric time series as distributed normally around its mean value. The power spectrum is the square of the modulus of the Fourier transform of the time series, that is the sum of the square of the real part and the square of the imaginary part , which gives a with 2 d.o.f. In other words, each point of the power spectrum is considered as a random variable which follows a distribution with 2 d.o.f. Duvall & Harvey (1986) showed that assuming that the excitation of the modes by turbulent convection is stochastic, the power spectrum must present a mean profile defined by a sum of Lorentzian functions as a function of the frequency:
where is the frequency, , and the height, width and eigenfrequency (1 per pair ), the splitting frequency and the inclination angle. The azimutal order is taken into account in the inclination angle dependence of the mode height (Gizon & Solanki 2003).
The likelihood function is the probability of the data set D with respect to the model characterised by the parameter set . The best fit corresponds to the model for which the data set probability is maximum. By assuming each point of the power as independent, which is wrong if smoothed, is:
where indicates the observed power spectrum value corresponding to the frequency bin . In practice, it is equivalent and easier to minimize the quantity . We used the finite-differencing (or Quasi-Newton) minimization routine of MATLAB to perform our likelihood calculations; we also provided the analytic expression of the derivative to optimize the routine efficiency.
Finally, the model function (Eq. 1) would be true if the power were only due to the stellar oscillations; actually, an additional term has to be added, in order to take into account the stellar noise associated with its surface granulation. This additional term may be written as a sum of a uniform noise level and a sum of background components as usually used in the helioseismology literature:
The power term , is generally fixed at 2 or 4. The SNR of a single mode is defined as the ratio of its height to the background .
2.2 When the MLE struggles to converge to a reasonable solution
Let us consider the fitting of the power spectrum of HD 181420, the first solar-like target to have been observed with CoRoT over a long run (156 days). The analysis of its power spectrum with the MLE has been carried out by Barban et al. (2009 hereafter B09), who included the results of several independent groups. In their work, they fitted the power spectrum by considering 14 overtones and degrees . Fitting 14 overtones until (Fig. 1) means fitting 75 parameters: 14 heights, 14 widths, 42 frequencies, a splitting, a stellar inclination, 2 stellar noise parameters (1 noise profile with a slope fixed at 2).
In Fig. 1 we present the global fitting of the power spectrum used by B09, obtained with the MLE for their scenario 1 (see also Sect. 5.2). The input parameters are close to the estimated values obtained in B09. It arises that when the signal-to-noise ratio of the power spectrum is too low, the MLE converges to a manifestly incorrect solution. In B09, the mode widths and heights of the external modes (lower and higher frequency) were rejected, in order to keep only the reliable fit values obtained with 10 overtones.
When the mode height becomes too low, the estimator tends to exactly fit some individual spikes of the power spectrum, that is, extremely high and narrow (very apparent in the frequency range [1100-1300] µHz). Hereafter we speak of “Dirac-like” convergence. Such a behavior is not surprising since the MLE consists of maximizing the similarity between the data and the model. The Bayesian approach allows us to avoid such solutions.
3 The maximum a posteriori approach
3.1 Bayesian spirit in the sky
Bayesian methods have been widely used in the deconvolution of astronomical images since the time of the restoration project for HST (Bertero et al. 1995), and also in cosmology (e.g. Trotta 2008) in order to fully exploit available reliable prior knowledge. Bayesian reasoning often arises unconsciously: e.g. in the solar case, if the MLE estimator gives a value of the rotational splitting near 0.82 Hz instead of Hz and an inclination of 0 instead of , we reject that solution because we know it is wrong, not on the basis of a statistical criterion. The same occurred in B09 when eliminating the estimates of the external mode heights and widths. They seemed wrong although they were statistically correct, since they occurred at a local minimum of the parameter space.
The Bayesian approach consists of using all of the prior information “I” that we have. Bayes theorem can be expressed as:
The Bayesian approach, compared with MLE, modifies the likelihood into a posterior distribution which takes into account the prior information . The is a normalization term. In this paper we consider only the simplest way of applying the Bayesian approach, i.e. the maximum a posteriori estimator (MAP):
Note the normalization term is not included and will not be taken into account in the following. Often, the prior is written as a Gaussian function centered around the expected value of the free parameter:
Therefore, the function to be minimized can be expressed as:
The MAP is the crudest Bayesian approach, but it allows us to improve the fits noticeably.
3.2 Prudence is required to avoid the traps laid along the way
In this section, we discuss the importance of carefully selecting the priors. Indeed, as can be anticipated when looking at Eq. 5, putting a stronger prior (i. e. narrow variance) may introduce a bias: the risk is to obtain exactly the value that we were expecting.
As an example, Fig. 2 (top) presents a simulation of the power spectrum of a single overtone featuring degrees , with SNR, which corresponds to the maximum SNR of HD 181420 (B09). We have tested the MAP fitting procedure by using a set of priors of the same value but with linearly increasing variance. The only free parameter that we have kept is the splitting frequency. Fig. 2 (bottom) shows the value of the estimated splitting frequency as a function of the prior variance. The “true” splitting corresponds to the value used to build the noisy power spectrum, while the “prior” splitting is the value that we think to be true and the “MLE” splitting is the value obtained with the MLE. When tends to 0, the output estimate tends to the prior value, whereas as tends to infinity the estimated output tends to the MLE estimate. So, the a priori probability does not help us find theÊ“true” value when the prior is wrong. It just tends towards the prior value as a function of the prior’s strength.
In other terms, the MAP does not push the solution to converge towards the local minimum that we may expect to be present at the “true” value of the parameter. Indeed, because of the stochastic nature of the power spectrum and because we have only one data set (the time series is one), the local minimum associated with the “true” value does not exist. The “true” value is only the mean value obtained on a large set of experiments. Thus, i) the use of MAP has to be performed only with a reliable prior (i.e. solar splitting frequency Hz) ii) the MAP has to be used only to prevent the fit from converging towards an incorrect solution (i.e. negative splitting).
4 Application to the global fittings
4.1 Putting a prior on the frequency splitting
The most common problem met with fitting the CoRoT solar-like oscillation spectra is the “Dirac-like” convergence (Sect. 2.2). According to the Bayesian approach, we must use all the information that may be deduced through reliable complementary observations or physical arguments. A proxy of the splitting frequency often can be deduced from the power spectrum itself. Indeed, on the power spectra of the 3 CoRoT targets HD 49933, HD 181420 and HD 181906, an excess power at very low frequency (Hz) is observed (Fig. 3 in A08 and B09, Fig. 4 in Garcia et al. 2009). This feature is associated with the rotational frequency, which appears directly in the power spectrum thanks to the motion of stellar spots. It corresponds to the star surface rotation. For physical reasons it is clear that the internal stellar rotation does not differ strongly with respect to the surface rotation, otherwise the star would not stay in state of equilibrium (all are main sequence stars). Hence, the splitting prior is determined by this low frequency excess power, with a confidence interval that has to be chosen ad-hoc, as a function of the width of the low-frequency peak (e. g. % confidence for Benomar et al. 2009).
4.2 Putting a constraint on the height
The second prior addresses the original purpose of the paper. We started by considering that the mode amplitude varies as a continuous function of the frequency. This is justified by theoretical studies (Houdek et al. 1999, Samadi & Goupil 2001 and Samadi et al. 2007) and has been observed in solar seismological data (e.g. Chaplin et al. 1998). The mode amplitude of a single p-mode is defined by its integral . As shown in Toutain & Appourchaux (1994), the mode height and width determinations are strictly correlated. So, constraining one of them is sufficient to constrain the amplitude.
In solar seismological data (e.g. GOLF), the height presents a “bell” shaped dependence as a function of the frequency (Fig. 3), while the width presents a more complex behavior: 2 slopes and a plateau (Fig. 4). On CoRoT data (e.g. HD 49933 in A08) the smoothed power spectrum shows that the height mainly follows a bell shaped profile, as for the Sun. It reflects the fact that the oscillation spectrum presents a beginning, an ending and a maximum in between. On the other hand, contrarily to the GOLF solar data in which the mode width is measurable with the naked eye, in the CoRoT data it is impossible to give an estimate of the width trend, because of the SNR.
It appeared logical to take such information into account. The difficulty is to find a simple mathematical and computational translation of it. The first idea is to fix the height close to a bell shaped profile. The drawback of such an approach is to introduce several hyper parameters to describe the bell shaped profile in the fitting process, which goes against the aim of reducing both the total number of parameters and the computational time. Therefore, we propose to replace the one height per overtone estimate by a continuous analytical function of the frequency, whose parameters are determined in the fitting process.
Figure 3 presents the solar mode heights as measured in 3 years of GOLF data, corresponding to a period of solar minimum activity. We have estimated the mode parameters by fitting the power spectrum with the MLE in groups of 4 overtones: 27 modes were fitted in the frequency range Hz. The purpose of the plot is to show that the height is well approximated by a simple Gaussian profile in the central part of the oscillation spectrum, that is, in the frequency range Hz, which corresponds to 10 overtones. A better fit would be obtained by complicating the fitting function, such as a combination of a Gaussian and an Lorentzian profile (pseudo-Voigt), but it is useless for the following analysis. As a comparison, in the HD 181420 CoRoT case, the frequency range in which the modes are high enough to be fitted is Hz, i.e. 14 overtones. So, by assuming the solar-like star to behave like the Sun, the “Gaussian height approximation” (hereafter GHA) is reliable for most of the power spectrum.
To test the validity of this approximation, a global fit taking into account 18 solar overtones was performed on 90 days of GOLF data in the frequency range Hz, which is larger that the confidence interval of the Gaussian fit. The comparison between the quadruplet mode fitting on the 1000 days of data and the global fit with the 90 days of data is done in Fig. 4. The width estimate is different in the details, but the global trend is conserved. The mean slope of the bandwidth on a log-log scale is equal to 3.84 with the not-approximated approach, whereas it is equal to 3.99 with the GHA, which represent an incorrect isestimate of almost 4%. Note that such an operation is a MAP approach with ; the probability that height follows a Gaussian profile is equal to 1.
The possibility of putting priors on the frequency localization was not considered, even though it could be done reliably by using the position of the power maxima measured on a smoothed power spectrum. The reason is that it has been observed that the mode frequencies are parameters for which no problems of convergence are met. Moreover, such a constraint on the frequency position implies a prior on the small separations and , for which we do not have valid constraint: until now, in all CoRoT solar-like observations, it was impossible to determine which of the peaks corresponds to degree or .
5 Application to CoRoT targets HD 181420 and HD 49933
We applied such a method to the global fitting of the two CoRoT solar-like targets that were globally fitted with MLE: HD 181420 (cf B09) and HD 49933 (cf A08). In this section, we present the specific choices of the priors, then the astrophysical results.
5.1 Priors on the splitting frequency, the Harvey profile and the height
According to the excess power observed in both stars at low frequency, the splitting prior was set at the maximum value of the peaks. For HD 181420, a quite large excess power appears in the frequency range Hz (cf B09), with a maximum at 4.2 Hz; the prior splitting frequency was fixed at Hz. For HD 49933, the excess power appears as a peak at almost Hz (cf A08); the prior splitting frequency was fixed at Hz.
Secondly, the low frequency stellar noise is well fitted by a sum of 3 background components, i.e. 3 Lorentzian functions of the frequency, centered on 0 and with 3 different amplitudes and widths (Harvey 1985). In our case, we do not consider the too low frequency range (Hz), hence one background component is sufficient. The power term (see Eq. 3) is set to 2. Despite the simple analytic expression, the background is the most sensitive contribution of the global fitting to be determined. Indeed, the consequences of an incorrect fitting are damaging for the oscillation parameter estimates. As an example, overestimating the background leads to an underestimate of the mode amplitudes and may strongly change the output values for the height even more for the width. Therefore, to diminish the risk of incorrectly estimating the background, we fitted it alone over a wide frequency range Hz, where the oscillation mode energy is negligible and the interval large enough to constrain the background parameters. Then, when fitting the oscillation spectrum, these parameters are set as free with the previous estimates as priors. The tolerance is set to 5% of the prior values, to permit small fluctuations due to the influence of the oscillation power when fitting the background profile alone.
For the GHA, the prior value of the width of the Gaussian function is set to the approximate width at half maximum estimated on a smoothed power spectrum over 25 Hz. For HD 181420 and HD 49933 it was set as Hz. Moreover, a constraint was placed on the centering of the Gaussian function on the frequency axis. In the same way, it arises from the position of a smoothed power spectrum (over 350 Hz). For HD 181420 it was set at Hz, while for HD 49933 it was set at Hz.
5.2 Results: the mode width increases with frequency
For each star, two mode identifications are tested because of our inability to identify with the naked eye the odd from the even modes. For HD 181420, the scenario deals with the identification of the peak at almost 1500 Hz: it corresponds to the pair in scenario 1 and to in scenario 2. For HD 49933, the scenario deals with the identification of the first considered peak, at almost 1250 Hz: it corresponds to the pair in scenario 1 and to in scenario 2. Figs. 5, 6, 7 and 8 present the results of the global fits obtained with the GHA for HD 181420 and HD 49933.
Beyond the fact that the Dirac-like convergence logically disappeared, the main result is the width trend as function of the frequency: in the 2 stars, whatever the scenario is, the mode width increases continuously. For HD 181420 the mean slope of the width on a log-log scale is equal to 2.8 in scenario 1 and 3.1 in scenario 2, while for HD 49933 the mean slope is equal to 3.2 and 2.5. These values are lower than the GOLF solar value, where the mean slope was estimated to be 4 in the central part of the oscillation spectrum. However, in Fig. 9 the widths of the two stars and the Sun are represented as functions of the radial order , instead of the frequency, to take into account the scaling factor due to the different size and density of the 3 stars. With such a representation, the mode width frequency dependences look more similar. Note that the slight excess widths observed at the maximum power of the 2 star spectra (radial order ) are the probable signature of the fact that the height profile differs slightly from a perfect Gaussian trend. It is probably sharper in this frequency range.
In addition, note that the mode frequency estimates changed noticeably between the MLE and the MAP fits, as illustrated for scenario 1 of HD 181420 (Fig. 1). In the MAP approach, modes with degrees and are clearly disconnected (mean small separation of about 5-10 Hz) while with the MLE the two modes interlace. Also, a careful analysis of the small separation , particularly in scenario 1 of HD 181420 and scenario 2 of HD 49933, shows that modes and are often inverted. A careful prior should be inserted in the fitting process in order to fix the small separation in both cases ( at left of and vice versa); then the likelihood value would determine which solution is the most likely. Such an analysis has to be carried out very carefully and the MCMC approach would be preffered to more accurately explore the parameter space. This is beyond the scope of the present work.
According to the suggestions of Appourchaux (2008) we have applied the maximum a posteriori approach to the global fitting of solar-like power spectra. It consists of a regularized maximum likelihood estimator, and hence is easy to implement. In particular, we have hypothesized that the mode height follows a Gaussian profile as function of the frequency. This was justified by the inability of the MLE to correctly fit the whole oscillation spectrum (Fig. 1). In other words, it is reasonable to differently process data with SNR than solar data with SNRÊ.
This method allowed us to globally fit the power spectra of 2 solar-like targets HD 181420 and HD 49933 over 14 overtones. The main result of such a method is that the mode width increases with frequency, as for the Sun. Moreover, it appears that the mean slope of the width as a function of the frequency on a log-log scale is similar to the solar value.
About the method, we conclude that with a power spectrum where the SNR lies between 3 and 10, a Bayesian approach is unavoidable, unless we are only able to properly fit the central part of the spectrum. The MAP estimator is very easy to apply, but as highlighted in Sect. 3.2, it has to be handled carefully. Such a method is a powerful tool to quickly and efficiently fit large samples such as those from Kepler (NASA) or for PLATO (ESA Cosmic Vision program).
We thank the CNES for the short term contract which allowed P. Gaulme to carry out the present work. The authors also thank F.X. Schmider and A. Ferrari (Fizeau, OCA, Nice) for discussions about the global fitting of oscillation spectra in the solar case and J. Leibacher for the patience in carefully reading the present manuscript.
|HD 181420||Scenario 1||Scenario 2|
|C ( s)||208.6||+9.6/-9.1||209.0||+18.8/-17.2|
|HD 49933||Scenario 1||Scenario 2|
- Appourchaux (2008) Appourchaux, T. 2008, Astronomische Nachrichten, 329, 485
- Appourchaux et al. (2008) Appourchaux, T., Michel, E., Auvergne, M., et al. 2008, A&A, 488, 705
- Barban et al. (2009) Barban et al. 2009, submitted to A&A
- Benomar et al. (2009) Benomar, O., Appourchaux, T., & Baudin, F. 2009, submitted to A&A
- Bertero et al. (1995) Bertero, M., Boccacci, P., & Maggio, F. 1995, International Journal of Imaging Systems Technology, 6, 376
- Chaplin et al. (1998) Chaplin, W. J., Elsworth, Y., Isaak, G. R., et al. 1998, MNRAS, 298, L7
- Duvall & Harvey (1986) Duvall, Jr., T. L. & Harvey, J. W. 1986, in Seismology of the Sun and the Distant Stars, 105–116
- Garcia & et al. (2009) Garcia, R. & et al. 2009, submitted to A&A
- Gizon & Solanki (2003) Gizon, L. & Solanki, S. K. 2003, ApJ, 589, 1009
- Harvey (1985) Harvey, J. 1985, in ESA Special Publication, Vol. 235, Future Missions in Solar, Heliospheric & Space Plasma Physics, ed. E. Rolfe & B. Battrick, 199–+
- Houdek et al. (1999) Houdek, G., Balmforth, N. J., Christensen-Dalsgaard, J., & Gough, D. O. 1999, A&A, 351, 582
- Howe (2009) Howe, R. 2009, ArXiv e-prints
- Samadi et al. (2007) Samadi, R., Georgobiani, D., Trampedach, R., et al. 2007, A&A, 463, 297
- Samadi & Goupil (2001) Samadi, R. & Goupil, M.-J. 2001, A&A, 370, 136
- Toutain & Appourchaux (1994) Toutain, T. & Appourchaux, T. 1994, A&A, 289, 649
- Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71