at the  towards other earths ii meeting held in porto portugal in 2014 september , one theme emerged that bares directly on the subject of this paper .
intrinsic stellar activity ( sa ) has become the main limiting factor for planet searches in both transit and radial velocity ( rv ) data .
sa includes stellar oscillations , granulation , spots and faculae / plagues and long - term magnetic activity cycles ( see , e.g. , saar & donahue 1997 ; santos et  al .
the presence of active regions on the stellar surface can show periodicities and amplitudes similar to those induced by real planetary signals .
. currently the best rv precision is at the m  s level but new spectrographs are under development like espresso and expres that aim to improve rv precision by a factor of approximately 10 over the current best spectrographs , harps and harps - n . clearly , the success of these developments hinges on our ability to distinguish true planetary signals from sa - induced signals . at the same time
xavier dumusque challenged the community to a large - scale blind test using simulated rv data at the 0.7  m  s level of precision , to understand the limitations of present solutions to deal with stellar signals and to select the best approach .
the rv challenge data were generated using a modification of the soap code , which is described in boisse , bonfils & santos ( 2012 ) and dumusque , boisse & santos ( 2014 ) .
this paper describes a new approach to rv analysis using apodized keplerian ( ak ) models to distinguish between planetary and sa - induced rv signals .
the methodology is developed in the next section using the initial challenge test data set for which the signals present were known in advance .
this is followed by our analysis and results for the first five rv challenge data sets .
the first step in the new approach employs ak models in which the semi - amplitude of the kepler rv model is multiplied by a symmetrical gaussian of unknown width  and with an unknown centre of the apodizing window ta .
the same model is used to fit both the planetary and sa - induced rv signals . since
a true planetary signal spans the duration of the data , the apodization time ,  , will be large while sa - induced signals generally vary on shorter time - scales .
one exception to this is the rapidly rotating ( p = 0.59  d ) active m4 dwarf gj 1243 which exhibits a starspot signal that remains stable over the entire baseline of the kepler mission ( davenport , hebb & hawley 2015 ) . as an additional check on the validity of planetary candidate signals , we employed a control diagnostic periodogram which is explained in section 2.4 .
in addition to the rv measurements , each of the 15 challenge data sets includes simultaneous observation of three sa diagnostics .
two of these come from additional information on the spectral line shape that are extracted from the cross - correlation function ( ccf ) , the average shape of all spectral lines of the star .
these two shape parameters are the ccf width ( full width at half - maximum , fwhm ) and bisector span ( bis ) .
the third diagnostic log  ( rhk ) is based on the ca  ii h&k line flux that is sensitive to active regions on the stellar surface .
a preliminary analysis of the first five data sets plus test data indicated a strong correlation between rv and the log  ( rhk ) diagnostic and slightly reduced correlation with the fwhm diagnostic . for the data sets analysed in this study , the correlation with the bis diagnostic
the general model for m ak signals also included a linear regression term between rv and the sa diagnostic log  ( rhk ) as well as an extra gaussian noise term with unknown standard deviation s. the correlation term is expected to be particularly useful in removing sa signals associated with the star 's long - term magnetic cycle ( dumusque et  al .
it is also expected to partially remove sa associated with the rotation of active regions ( spots and plagues ) . for an m signal ak model ,
the predicted rv is given by 
 ( 1)\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } \begin{eqnarray * } y(t_i ) & = & v + \sum _ { j=1}^m \big [ k_j \exp \big[-\frac{(t_i - t_{ja})^2}{2 \tau _ j^2}\big ] \nonumber \\ & & \times ( \cos \lbrace \theta ( t_i+\chi p_j)+\omega _ j\rbrace + e_j \cos \omega _ j ) \big ] \nonumber \\ & & + \ ,
the model also includes an extra gaussian white noise term with unknown standard deviation s which can account for stellar jitter .
the equation involves 7  m + 3 unknown parameters 
 v = a constant velocity; = the log  ( rhk ) correlation parameter;kj = velocity semi - amplitude for the jth signal;j = apodization time constant for the jth signal;tja = centre of apodizing window for the jth signal;pj = the period for the jth signal;ej = the eccentricity for the jth signal;j = the longitude of periastron for the jth signal;j = the fraction of the period , prior to the start of data taking , that periastron occurred at . for a planetary signal
, jpj = the number of days prior to ti = 0 that the star was at periastron , for an orbital period of pj days . the time ti = 0 corresponds to the unweighted average of the observation times.s = standard deviation of extra gaussian noise term .
v = a constant velocity ;  = the log  ( rhk ) correlation parameter ; kj = velocity semi - amplitude for the jth signal ; j = apodization time constant for the jth signal ; tja = centre of apodizing window for the jth signal ; pj = the period for the jth signal ; ej = the eccentricity for the jth signal ; j = the longitude of periastron for the jth signal ; j = the fraction of the period , prior to the start of data taking , that periastron occurred at . for a planetary signal ,
jpj = the number of days prior to ti = 0 that the star was at periastron , for an orbital period of pj days .
the item (ti + jpj ) = the true anomaly , the angle of the star in its orbit relative to periastron at time ti .
as described in more detail in gregory ( 2011 ) , we employed a re - parametrization of  and  to improve the markov chain monte carlo ( mcmc ) convergence speed motivated by the work of ford ( 2006 ) .
the two new parameters are  = 2 +  and  = 2  . parameter  is well determined for all eccentricities . although  is not well determined for low eccentricities , it is at least orthogonal to the  parameter .
we use a uniform prior for  in the interval 0 to 4 and uniform prior for  in the interval 
2 to + 2. this ensures that a prior that is wraparound continuous in (  ,  ) maps into a wraparound continuous distribution in (  ,  ) . to account for the jacobian of this re - parametrization , it is necessary to multiply the bayesian integrals by a factor of ( 4 ) , where np = the number of periods in the model
the ak models were explored using an automated fusion mcmc algorithm ( fmcmc ) , a general purpose tool for non - linear model fitting and regression analysis ( gregory 2011 , 2013 ) .
fmcmc is a special version of the metropolis mcmc algorithm that incorporates parallel tempering , genetic crossover operations and an automatic simulated annealing .
each of these features facilitates the detection of a global minimum in chi - squared in a highly multi - modal environment . by combining all three
the fmcmc is controlled by a unique adaptive control system that automates the tuning of the mcmc proposal distributions for efficient exploration of the model parameter space even when the parameters are highly correlated .
the ak models combined with the fmcmc algorithm constitute a multi - signal ak periodogram .
a more detailed description of fmcmc is available in chapter 1 of the  supplement to bayesian logical data analysis for the physical sciences , a free supplement available in the resources section of the cambridge university press website for my textbook  bayesian logical data analysis for the physical sciences : a comparative approach with mathematica support. a mathematica implementation of fmcmc is also available from the resource section .
chapter 1 also includes a comparison of three marginal likelihood estimators used for bayesian model comparison and concludes in favour of the nested restricted monte carlo ( nrmc ) estimator which is used in this work .
the supplement includes a detailed discussion of the priors adopted by this author for exoplanet rv analysis .
two different frequency priors were discussed : ( a ) \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ p(f|i ) = 1/\sqrt{f } \times 0.5/(\sqrt{f_h } -\sqrt{f_l})$\end{document } and ( b ) p(f|i ) = 1/f  1/ln  ( fh / fl ) , where fh and fl are the prior upper and lower frequency bounds , respectively . in this work , ( a ) was used primarily for parameter estimation purposes which helps in the detection of shorter period signals while the scale - invariant form ( b ) was adopted for model comparison purposes which utilized nrmc method for marginal likelihood estimation .
a new eccentricity prior was utilized in this work intermediate between kipping 's low - period ( p < 382.3  d ) case and his average for all 396 planets ( kipping 2013 ) .
the equation for the new eccentricity prior is 
 ( 2)\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } \begin{equation * } p(e|i ) = 7.354*(1 - e^{0.3})^{1.5}. \end{equation*}\end{document}the two additional parameters for a one - signal ak model are the apodization time ,  , and the apodization function centre time ta .
we used a scale - invariant prior for  and a uniform prior for ta .
the upper limit on  was set at four times the duration of the data set which corresponds to a maximum variation in the window height of approximately 1  per  cent across the data set for a window with a centre at the middle of the data or by 3  per  cent for one centred at either end of the data .
the lower limit was set to the data duration divided by 40 which corresponds to 37  d. some of the calculations were redone using an upper limit of eight times the duration and a lower limit of 10  d with no noticeable improvement in our ability to distinguish sa and planetary signals . much more analysis would be needed to select an optimum lower limit for . the upper and lower limits on the uniform ta prior were t to + t relative to the mean data sample time , where t is the data duration .
as we will shall see , the primary role of the ak models is to distinguish planetary signal candidates from sa signals .
suppose the results indicate that k of the signals are sa and m  k are planetary candidate signals .
final model parameter estimates and model comparisons are based on subsequent runs using a model of k aks and m  k kelperians .
the test data set contained a single planetary signal in a circular orbit with p = 16.0  d , k = 1.5  m  s and phase = 0.0 , together with a variety of sa signals . fig . 
1 shows a comparison between the raw rv test data and two diagnostics , the fwhm and log  ( rhk ) , and their generalized lomb 
, the three signals appear highly correlated but close inspection of the periodograms indicates differences between the rv data and the diagnostics in the vicinity of the known 16  d planetary signal .
the gls periodogram p( ) measures the relative -reduction , p( ) , as a function of frequency  and is normalized to unity by \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \chi _ 0 ^ 2$\end{document } ( the  for the weighted mean of the data ) .
the quantity p( ) is normalized to unity with p( ) = 0 indicating no improvement of the fit and p( ) = 1 a 
for all the gls periodogram plots , the abscissa is expressed as a period where \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ p = \frac{2 \pi } { \omega } $ \end{document}. frequentist p - values of 0.1 , 0.01 , 0.001 are shown in the upper - right periodogram .
the left column of panels are plots of the raw rv test data and two diagnostics , the fwhm and log  ( rhk ) .
the test data set is known to have a single planetary signal with a period of 16  d. for each plot , the mean value has been subtracted .
for all the gls periodogram plots , the abscissa is expressed as a period where \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ p = \frac{2 \pi } { \omega } $ \end{document}. in fig .  2
, the upper - left panel shows the rv data after removing the best linear regression fit with log  ( rhk ) as the independent variable .
this is referred to as rv ( rhk corrected ) in the panel caption . removing the correlation with log  ( rhk ) reduced the standard deviation of the rv data from 8.55 to 2.72  m  s. we also tried removing a linear regression model with log  ( rhk ) and the fwhm as two independent regression variables .
this resulted in only a very slight reduction in the standard deviation of the residuals to 2.71  m  s. for this reason , in the ak model fitting described below , we only include a single log  ( rhk ) regression term as indicated in equation ( 1 ) .
the left column of panels are plots of the test data rv and fwhm after removing the best linear regression fits to log  ( rhk ) ( rhk corrected ) .
it proved useful to employ the fwhm diagnostic data as a control after removing the best linear regression fit to log  ( rhk ) .
this reduced the standard deviation of the raw fwhm diagnostic from 26.2 to 5.6  m  s. these control data appear in the lower - left panel of fig .  2 and
it provides another check on validity of planetary candidate signals identified with the ak model .
a valid planetary signal is not normally expected to have a significant counterpart in the gls periodogram of the control data .
the gls periodograms of the rv ( rhk corrected ) data and control are shown in the two rhs panels .
the periodogram of rv ( rhk corrected ) exhibits a dominant peak with a period of 16  d and a forest of smaller peaks .
the periodogram of the control exhibits multiple lower level peaks which were interpreted as residual sa signals or instrumental systematics .
it also proved useful to construct a differential form of the gls periodogram of the rv residuals for selected period regions as is shown in fig . 
the black trace is from the upper - right periodogram in fig .  2 .
the dark grey trace is the negative of the periodogram of the control and the light grey trace shows the difference , the black trace plus the dark grey trace . the differential gls is shown for six period intervals .
the strongest feature occurs at 16  d and on either side are the two 1  yr aliases .
1 d aliases can be seen in the vicinity of p = 0.94  d and 1.06  d. for all three , the light grey trace and the black trace coincide closely near the positive peaks , consistent with a planetary signal or its alias .
signals in common to both black and dark grey traces ( e.g. as occurs near p = 6.3  d ) indicate sa . at each stage in the analysis , the gls periodogram ( fig .  2 ) and the differential gls periodogram ( fig .  3 ) of the current rv residuals provide guidance on the choice of the next period to include in the ak model .
a differential gls periodogram of selected period regions for the test data based on the rv and fwhm control ( both rhk corrected ) .
the lower - left panel shows the span of the apodization window within the overall data window for each signal ( grey trace for the maximum a posteriori ( map ) values of  and ta , black for a representative set of samples which in this case is completely hidden below the grey ) .
the lower and upper boundaries of the apodization window are defined by ta   and ta +  , respectively .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the one - signal ak periodogram of the test data .
the upper right shows log10[prior  likelihood ] versus period showing the one period detected .
the middle left shows the values of the one unknown period parameter versus iteration number .
the lower left shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is hidden below the grey ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
5 shows the one - signal ak model residuals [ i.e. m = 1 in equation ( 1 ) ] .
the p = 16  d signal and its 1  yr and daily aliases have been successfully removed .
the feature at p = 6.3  d is now the highest peak in the residuals .
the panels show the one - signal ak model residuals for the test data set ( top left ) , its gls periodogram ( top right ) and the differential periodogram for selected period ranges .
the analysis proceeds by repeatedly adding an additional ak signal to the model with an initial period estimate based on the previous residuals and examining the fmcmc results plus the gls and differential periodograms of the residuals .
note that the linear log  ( rhk ) correlation is not removed prior to each fit but is accomplished by the    rhk(ti ) term in equation ( 1 ) as part of each fit .
this process continues until the p - value > 0.01 for the highest peak in the gls periodogram of the residuals , and the fmcmc results for the next more complicated model indicate that there is no well - defined solution for an additional signal . in the interest of space , we jump to our final five ak signal model results .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the five - signal ak periodogram of the test data .
the upper right shows log10[prior  likelihood ] versus period showing the five periods detected .
the middle left shows the values of the five unknown period parameter versus iteration number .
the lower - left panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
the lower and upper boundaries of the jth apodization windows are defined by tja  j and tja + j , respectively .
the lower - right panel is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
based on the span of the apodization window , all of the signals , with the exception of p = 16  d , are consistent with sa . as a check on the planetary interpretation of the 16  d signal , fig . 
3 indicates that p( ) for the control trace is within 0.03 of zero for a period of 16  d. at this point in the analysis , the peak in the gls periodogram of the five - signal ak fit residuals at p = 0.95  d was down to p( ) = 0.045 .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the five - signal ak periodogram of the test data .
the upper right shows log10[prior  likelihood ] versus period showing the five periods detected .
the middle left shows the values of the five unknown period parameter versus iteration number .
the lower left shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is hidden below the grey ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta , where the symbols refer to the different signal periods specified in the top panels . to obtain the best p = 16  d parameter estimates , the m ak model terms in equation ( 1 ) were replaced by four aks plus one ordinary keplerian to model the 16  d planetary signal , and the data refit .
the starting fmcmc parameter values were taken from the best - fitting five ak signal model .
7 shows a sample of resulting 16  d parameter correlation plots , which in a bayesian analysis reflect our current state of knowledge for the 16  d planetary signal in the test data set .
parameter correlation plots for the p = 16  d planetary signal in the test data set based on a model consisting of four aks plus one ordinary keplerian .
the p and k values have been normalized by dividing by their true values to convert to ratios . since the true eccentricity = 0
, the relevant phase for comparison is our re - parametrized phase term  = 2 + . in our analysis , the phase is referenced to the mean observation time .
this was converted to a phase referenced to the time of the first sample and added to an additional /2 to compare with the xavier dumusque 's definition of phase which refers to a sine function while ours refers to a cosine function .
it is not too meaningful to draw any conclusions about parameter biases on the basis of this one data set . later in section 4 , we compare the k and p parameter estimates derived from the m ak model fits to the k ak plus ( m  k ) kepler model fits for the six data sets analysed ( including the test data ) .
it is clear that ak model estimates of the k parameter are significantly biased to larger values compared to their known true values . as a result
, our final parameter estimates for candidate planets were based on the k ak plus ( m  k ) kepler model fits , where k is the number of signals attributed to sa .
ideally , we would like to see the bayesian marginals span the true values of p , k and phase slightly better than indicated by fig .  7 .
this might be an indication that the true scatter is underestimated possibly as a result of correlated noise in the residuals .
one way of exploring the residuals is to compute the autocorrelation function , (j ) , 
 ( 3)\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } \begin{equation * } \rho ( j ) = \frac{\sum _ { \rm overlap } [ ( x_{i}-\overline{x})\ ( x_{i+j}-\overline{x})]}{\sqrt{\sum _ { \rm overlap } ( x_{i}-\overline{x})^2 } \times \sqrt{\sum _ { \rm overlap } ( x_{i+j}-\overline{x})^2 } } , \end{equation*}\end{document}where xi is the ith residual , j is the lag and \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \overline{x}$\end{document } is the mean of the samples in the overlap region . the top row of fig . 
8 shows the autocorrelation function of the ( rhk corrected ) rv test data and the raw log  ( rhk ) diagnostic .
note that the highest peak in the autocorrelation function for the log  ( rhk ) diagnostic corresponds roughly with the star 's true rotation period of 25.05  d ( vertical dashed line ) .
the left - hand panel in the second row shows the autocorrelation function for the rv residuals for the four ak plus one kepler model including the log  ( rhk ) regression . clearly , after removing
a model consisting of five signals and the log  ( rhk ) regression , the autocorrelation of the residuals is looking much closer to white noise but there is still evidence for a small positive autocorrelation for a lag of 1 d. the remaining panel shows a plot of the number of such sample pairs as a function of the lag . because the data are not uniformly sampled , for each lag all sample pairs that differed in time by this lag 0.1  d were utilized .
the top row shows the autocorrelation function of the ( rhk corrected ) rv test data and the raw log  ( rhk ) diagnostic .
the left - hand panel in the second row shows the autocorrelation function of the rv residuals after removing a model consisting of four ak signals plus one kepler signal and the log  ( rhk ) regression .
the bottom panel shows a plot of the number of such sample pairs as a function of the lag . in this section ,
we summarize the steps in the methodology and illustrate its success in modelling the both long - term and short - term sa - induced rv variations .
remove the best linear regression fit from the raw rv data to obtain the rv ( rhk corrected ) data , using the sa diagnostic log  ( rhk ) as the independent regression variable.compute the gls periodogram of rv ( rhk corrected ) and p - value levels corresponding to 0.1 , 0.01 , 0.001 .
the highest peak in this periodogram is used as the initial estimate for the period of the ak model signal.remove the best linear regression fit from the raw fwhm data to obtain the control fwhm ( rhk corrected ) data , using the log  ( rhk ) diagnostic as the independent regression variable.compute the gls periodogram of control fwhm ( rhk corrected).compute the differential gls periodogram for selected period regions , showing on the same plot the periodogram of rv ( rhk corrected ) in black , the periodogram of the control fwhm ( rhk corrected ) plotted negatively and shown in dark grey , and the their sum shown in grey .
planetary candidate signals for which the grey trace and black trace are significantly different ( i.e. the dark grey trace has a p( ) more negative than 0.05 ) indicate that the signal is possibly sa in origin.fit the model described by equation ( 1 ) to the raw rv data , in our case using the fmcmc algorithm .
use the highest peak in the gls periodogram of the rhk - corrected rv data as the initial estimate for the period of the ak model signal , regardless of whether the differential gls periodogram indicates it is sa in origin . we are using the aks to model both types of signals.compute the gls and differential periodograms of the fit residuals.repeat the last two steps by including an additional ak signal with an initial period estimated from the largest peak in the gls periodogram of the previous fit residuals.repeat the previous step until the p - value of the highest peak in the residuals is greater than 0.01 and the fmcmc fit indicates that there is no well - defined solution for an additional signal.examine the span of the apodization window for the map parameter values and representative samples for each signal extracted .
the lower and upper boundaries of the jth apodization windows were chosen to be tja  j and tja + j , respectively.if the apodization window spans the data duration , and p( ) for the dark grey trace in the differential periodogram of the residuals is less negative than 0.05 , then treat the signal as a candidate planetary signal .
if the apodization window spans the data duration and p( ) for the dark grey trace is more negative than 0.05 , then treat the signal as a possible candidate planetary signal . if the apodization window does not span the data duration , classify it as a sa signal.maintain a list of the order in which the signals were extracted and the signal designation as p for planetary candidate , p ? for possible planetary candidate and sa for stellar activity signal.suppose
there are m signals in total , k of which are sa signals and ( m  k ) are planetary candidates or possible planetary . fit a model with k ak signals and ( m  k ) ordinary keplerian signals together with the   
rhk(ti ) term of equation ( 1 ) to derive final parameter estimates for the planetary and possible planetary candidates.compute bayes factors as a further test of the validity of any planetary candidate .
remove the best linear regression fit from the raw rv data to obtain the rv ( rhk corrected ) data , using the sa diagnostic log  ( rhk ) as the independent regression variable .
compute the gls periodogram of rv ( rhk corrected ) and p - value levels corresponding to 0.1 , 0.01 , 0.001 .
the highest peak in this periodogram is used as the initial estimate for the period of the ak model signal . remove the best linear regression fit from the raw fwhm data to obtain the control fwhm ( rhk corrected ) data , using the log  ( rhk ) diagnostic as the independent regression variable . compute the gls periodogram of control fwhm ( rhk corrected ) .
compute the differential gls periodogram for selected period regions , showing on the same plot the periodogram of rv ( rhk corrected ) in black , the periodogram of the control fwhm ( rhk corrected ) plotted negatively and shown in dark grey , and the their sum shown in grey .
planetary candidate signals for which the grey trace and black trace are significantly different ( i.e. the dark grey trace has a p( ) more negative than 0.05 ) indicate that the signal is possibly sa in origin .
fit the model described by equation ( 1 ) to the raw rv data , in our case using the fmcmc algorithm . use the highest peak in the gls periodogram of the rhk - corrected rv data as
the initial estimate for the period of the ak model signal , regardless of whether the differential gls periodogram indicates it is sa in origin .
compute the gls and differential periodograms of the fit residuals . repeat the last two steps by including an additional ak signal with an initial period estimated from the largest peak in the gls periodogram of the previous fit residuals . repeat the previous step until the p - value of the highest peak in the residuals is greater than 0.01 and the fmcmc fit indicates that there is no well - defined solution for an additional signal .
examine the span of the apodization window for the map parameter values and representative samples for each signal extracted .
the lower and upper boundaries of the jth apodization windows were chosen to be tja  j and tja + j , respectively .
if the apodization window spans the data duration , and p( ) for the dark grey trace in the differential periodogram of the residuals is less negative than 0.05 , then treat the signal as a candidate planetary signal .
if the apodization window spans the data duration and p( ) for the dark grey trace is more negative than 0.05 , then treat the signal as a possible candidate planetary signal . if the apodization window does not span the data duration , classify it as a sa signal . maintain a list of the order in which the signals were extracted and the signal designation as p for planetary candidate , p ? for possible planetary candidate and sa for stellar activity signal .
suppose there are m signals in total , k of which are sa signals and ( m  k ) are planetary candidates or possible planetary .
fit a model with k ak signals and ( m  k ) ordinary keplerian signals together with the   
rhk(ti ) term of equation ( 1 ) to derive final parameter estimates for the planetary and possible planetary candidates .
this is discussed further starting in section 3.2 . in the upper - left panel of fig . 
9 , the ( open circles ) represent the raw rv test data minus the map fit for the 16  d planet which we designate as rv1plancorr .
the dots show how much of the rv1plancorr can be accounted for by a correlation with the sa diagnostic log  ( rhk ) , based on a weighted least - squares linear regression fit using log  ( rhk ) as the independent regression variable .
the starred points show the difference between the open circles and the dots . removing the correlation with log  ( rhk ) reduced the standard deviation of rv1plancorr from 8.38 to 2.29  m  s.
clearly , the log  ( rhk ) regression fit tracks the long - term sa very well .
close inspection of the seasonal panels indicates that the regression fit only provides a first - order tracking of the short - term variations in rv arising from individual active regions ( spots and plagues ) .
the open circles in the upper - left panel show the raw rv test data minus the map fit for the 16  d planet .
the dots show how much of this remainder can be accounted for by a correlation with the sa diagnostic log  ( rhk ) .
the other panels show expanded versions , one for each observing season . in the upper - left panel of fig . 
the dots show the map fit of the full model consisting of four aks plus one ordinary keplerian ( to model the 16  d planetary signal ) , together with the linear regression term using log  ( rhk ) as the independent regression variable .
the starred points show the difference between the open circles and the dots . employing the full model reduced the standard deviation of from 8.55 to 1.35  m  s. the other panels show expanded versions , one for each observing season .
clearly , the full model provides a much better fit to the rv data including the short - term variations arising from individual active regions .
the open circles in the upper - left panel show the raw rv test data .
the dotted points show how much of each rv measurement can be accounted for by the full model ( equation 1 ) fit .
the procedures outlined in section 2 were applied to the first five rv challenge data sets and the results are given below .
table  1 shows the order in which the signals emerged in successive ak models together with the signal designation as planetary ( p ) or stellar activity ( sa ) in the final column .
order of signals extracted from rv 1 together with approximate parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) .
comparison of the raw rv and fwhm and log  ( rhk ) diagnostics indicated that there was a clear correlation between the three .
the raw rv 1 data had a standard deviation of 5.6  m  s. after removing the best linear regression fit with log  ( rhk ) as the independent variable , the standard deviation was reduced to 3.0  m  s. the panels on the left of fig . 
11 show the rv data and fwhm ( control ) with the best linear regression fits to log  ( rhk ) removed ( rhk corrected ) .
note that for each ak model fit , the correlation parameter  between rv and the sa diagnostic log  ( rhk ) was included as a fit parameter according to equation ( 1 ) .
12 shows the differential gls periodogram for the rv 1 data ( rhk corrected ) for selected period ranges .
the strongest signal at p = 33  d has no significant counterpart in the control .
the rv 1 data and fwhm ( control ) after removing the best linear regression fits to log  ( rhk ) ( rhk corrected ) together with their gls periodograms on the right .
the differential gls periodogram for the rv 1 data ( rhk corrected ) for selected period ranges .
the lower - left panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
the lower and upper boundaries of the jth apodization window are defined by tja  j and tja + j , respectively .
it is clear from this that there are three planetary candidates with periods of 9.89 , 23.4 , 33.3  d whose apodization windows span the duration of the data .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the six signal ak periodogram of the rv 1 data .
the upper right shows log10[prior  likelihood ] versus period showing the six periods detected .
the middle left shows the values of the six unknown period parameters versus iteration number .
the lower left shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta , where the symbols refer to the different signal periods specified in the top panels .
the weighted rms residual = 1.44  m  s for the six ak signal model .
the standard deviation of the extra gaussian white noise term s = 1.26  m  s. the mean measurement uncertainty was 0.674  m  s. fig . 
14 shows the rv 1 six - signal ak model residuals [ i.e. m = 6 in equation ( 1 ) ] , its gls periodogram ( top right ) and the differential periodogram for selected period ranges .
the dashed line representing the p - value = 0.01 is not seen as it is just above the top of the graph .
these results are consistent with the conclusion that the p = 9.89 , 33.35 , 23.4  d signals are bona fide planetary candidates .
the panels show the rv 1 six - signal ak model residuals ( top left ) , its gls periodogram ( top right ) and the differential periodogram for selected period ranges .
table  2 shows the true parameters of the five planetary signals that were employed in the rv 1 simulation .
clearly , the three planetary signals with k > 1 were recovered in our analysis but there is no clear evidence for the two weaker planetary signals at p = 112  and  273  d in the gls residuals of fig .  14 . the true p , k , e parameters of the five planetary signals that were employed in the rv 1 simulation . the top row of fig . 
15 shows the autocorrelation function of the raw rv data and the log  ( rhk ) diagnostic . in this particular case , the highest peak in the autocorrelation function for the log  ( rhk )
diagnostic corresponds closely with the star 's rotation period of 25.05  d. sa signals frequently occur at the rotation period and harmonics .
the second row shows the autocorrelation function for the rv data with just the ln(rhk ) regression removed ( rhk corrected ) and the six ak model residuals . in this case , after the six ak signals and log  ( rhk ) regression are removed , the autocorrelation of the residuals is looking close to white noise .
the top row shows the autocorrelation function of the raw rv 1 data and the log  ( rhk ) diagnostic .
the second row shows the autocorrelation function for the rv data with log  ( rhk ) regression removed ( rhk corrected ) and the six ak model residuals . table  3 gives the final parameter estimates for the three planetary candidates derived from a model consisting of three aks for the sa signals and three keplerians for the planetary candidates .
this analysis yielded the same set of periods found in the six ak model fit .
the quoted eccentricity values are for the mode while for the other parameters we quote the median .
the errors on semi - major axis , a , and msin  i do not include the uncertainty in the mass of the host star .
the last row gives the date of the periastron passage that occurs just prior to our reference date jdb = 55855.666  d which is the unweighted mean observation time for rv 1 .
rv 1 parameter estimates and 68  per  cent credible boundaries for three planetary candidates discussed in the text .
the quoted eccentricity values are for the mode while for the other parameters we quote the median .
the errors on semi - major axis , a , and msin  i do not include the uncertainty in the mass of the host star . see the note on periastron passage in the text . to test the significance of the three planetary candidates
, we computed bayes factors comparing the following four models : ( a ) zero keplerian signals with only the log  ( rhk ) regression removed , ( b ) one 33.35  d keplerian signal plus the log  ( rhk ) regression , ( c ) two ( 33.35 and 9.89  d ) keplerian signals plus log  ( rhk ) regression , and ( d ) three ( 33.35 , 9.89 , 23.4  d ) keplerian signals plus log  ( rhk ) regression .
relative to model ( d ) , the bayes factors are 
 \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } 
 \begin{equation * } 
 3.8 \times 10^{-66}:7.2 \times 10^{-45}:1.6 \times 10^{-26}:1.0 . 

\end{equation * } 
 \end{document}the bayes factors were computed from model marginal likelihoods based on the nrmc method ( gregory & fischer 2010 ; gregory 2011 , 2013 , and in more detail in section 1.6 of the  supplement to bayesian logical data analysis for the physical sciences , available in the resources section of the cambridge university press website for my textbook  bayesian logical data analysis for the physical sciences : a comparative approach with mathematica support ) .
comparison of the raw rv and fwhm and log  ( rhk ) diagnostics indicated that there was a clear correlation between the three .
the raw rv 2 data had a standard deviation of 8.58  m  s. after removing the best linear regression fit with log  ( rhk ) as the independent variable , the standard deviation was reduced to 3.95  m  s. the two rows of fig . 
16 show the rhk - corrected rv and fwhm control together with their gls periodogram on the right .
17 shows the differential periodogram for the rhk - corrected rv data for selected period ranges .
the strongest signal at p = 10.64  d has no significant counterpart in the control .
the rv 2 data and fwhm ( control ) after removing the best linear regression fits to log  ( rhk ) together with their gls periodograms on the right .
the differential gls periodogram for the rhk - corrected rv 2 data for selected period ranges .
the lower - left panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
it is clear from this that there are three planetary candidates with periods of 3.77 , 10.64 and 75.56  d whose apodization windows span the duration of the data .
18 although the apodization time constant of 800  d for the 10.64  d signal makes it a borderline p candidate . also , the presence of sa activity at a period of 11.32  d , in the close vicinity of the 10.64  d signal , perhaps called into question a planetary interpretation of the 10.64  d signal . for the competition ,
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the eight - signal ak periodogram of the rv 2 data .
the upper right shows log10[prior  likelihood ] versus period showing the eight periods detected .
the middle left shows the values of the eight unknown period parameters versus iteration number .
the lower left shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
the weighted rms residual = 1.42  m  s for the eight ak signal model .
the standard deviation of the extra gaussian white noise term s = 1.25  m  s. the mean measurement uncertainty was 0.674  m  s. fig . 
19 shows the rv 2 eight - signal ak model residuals ( i.e. m = 8 in equation 1 ) , its gls periodogram ( top right ) and the differential periodogram for selected period ranges .
the highest spectral peak in the residuals has a p - value just < 0.1 .
the dashed line representing the p - value = 0.01 is not seen as it is above the top of the graph .
these results are consistent with the conclusion that the p = 3.77 , 10.64 , 75.56  d signals are bona fide planetary candidates . the panels show the rv 2 eight - signal ak model residuals ( top left ) , its gls periodogram ( top right ) and the differential periodogram for selected period ranges .
table 4 lists ( in order of extraction ) the signals , their nominal parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) based on the mcmc parameter estimates . for the competition ,
the 10.64  d signal was reported as a probable planetary signal for the reasons mentioned above .
a bayes factor analysis concluded that all three signals were highly significant compared to the next simpler model .
order of signals extracted from rv 2 together with approximate parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) .
table 5 shows the true p , k , e parameters of the five planetary signals that were employed in the rv 2 simulation .
clearly , the three planetary signals with k > 1 were recovered in our analysis but there is no clear evidence for the two weaker planetary signals at p = 5.79  and  20.16  d in the gls residuals of fig .  19 . the true p , k , e parameters of the five planetary signals that were employed in the rv 2 simulation .
the same set of periods was obtained using a model consisting of five ak signals and three keplerian signals .
final parameter estimates for the three planetary candidates were derived from this latter model and are given in table  6 .
the quoted eccentricity values are for the mode while for the other parameters we quote the median .
the errors on semi - major axis , a , and msin  i do not include the uncertainty in the mass of the host star .
the last row gives the date of the periastron passage that occurs just prior to our reference date jdb = 55855.6693  d which is the unweighted mean observation time for rv 2 .
rv 2 parameter estimates and 68  per  cent credible boundaries for three planetary candidates discussed in the text .
the errors on semi - major axis , a , and msin  i do not include the uncertainty in the mass of the host star . to test the significance of the first two planetary candidates , we computed bayes factors comparing the following models : ( a ) zero keplerian signals with only the log  ( rhk ) regression removed , ( b ) one 3.77  d keplerian signal plus the log  ( rhk ) regression , ( c ) two ( 3.77 and 10.64  d ) keplerian signals plus log  ( rhk ) regression . the bayes factors relative to model ( c ) are 
 \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } 
 \begin{equation * } 
 5.6 \times 10^{-60}:9.3 \times 10^{-37}:1.0 . 

\end{equation * } 
 \end{document}the three - planet model was compared to the two - planet model by comparing the two models : ( a ) three aks ( p = 8.33 , 12.5 , 24.8  d ) plus two keplerians ( 3.77 , 10.64  d ) plus the log  ( rhk ) regression to ( b ) three aks ( p = 8.33 , 12.5 , 24.8  d ) plus three keplerians ( 3.77 , 10.64 , 75.7  d ) plus the log  ( rhk ) regression . the bayes factor relative to model ( b ) is 3.4  10:1.0 .
comparison of the raw rv and fwhm and log  ( rhk ) diagnostics indicated that there was a clear correlation between the three .
the raw rv 3 data had a standard deviation of 10.85  m  s. after removing the best linear regression fit with log  ( rhk ) as the independent variable , the standard deviation was reduced to 5.49  m  s. the two rows of fig . 
20 show the rhk - corrected rv data and fwhm control together with their gls periodogram on the right .
21 shows the differential periodogram for the rhk - corrected rv data for selected period ranges .
the strongest signal at p = 48.8  d has no significant counterpart in the control .
strong yearly and daily aliases of the 48.8  d signal are also evident . the rv 3 data and fwhm ( control ) after removing the log  ( rhk ) diagnostics ( rhk corrected ) together with their gls periodograms on the right .
the differential gls periodogram for the rhk - corrected rv 3 data for selected period ranges .
table 7 lists ( in order of extraction ) the signals , their nominal parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) based on the mcmc parameter estimates .
the weighted rms residual = 1.79  m  s for the six ak signal model .
the standard deviation of the extra gaussian white noise term s = 1.66  m  s. the mean measurement uncertainty was 0.674  m  s. three signals
( p = 17 , 48.8 and 1100  d ) were initially classified as planetary although there was a noticeable signature ( p( ) > 0.05 ) of the 17  d signal in the fwhm control trace which resulted in a probable p classification .
order of signals extracted from rv 3 together with approximate parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) .
22 , after signals at 12.5 , 17 , 48.8 and 333  d were first removed .
the same set of periods was obtained using a model consisting of three ak signals plus three keplerian signals and the log  ( rhk ) regression term .
the one difference from the six ak signal fit is that now the signal near p  1100  d has a low - eccentricity tail extending to approximately double that period . to test the significance of the first two planetary candidates
, we computed bayes factors comparing the following models : ( a ) zero keplerian signals with only the log  ( rhk ) regression removed , ( b ) one 48  d keplerian signal plus the log  ( rhk ) regression , ( c ) two ( 48 and 17  d ) keplerian signals plus log  ( rhk ) regression .
the bayes factors relative to model ( c ) are 
 \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } 
 \begin{equation * } 
 2.7 \times 10^{-101}:1.2 \times 10^{-54}:1.0 . 

\end{equation * } 
 \end{document}the three - planet model was compared to the two - planet model by comparing the two models : ( a ) two aks ( p = 12.5 , 333  d ) plus two keplerians ( 48 , 17  d ) plus the log  ( rhk ) regression to ( b ) two aks ( p = 12.5 , 333  d ) plus three keplerians ( 48 , 17 , 1100  d ) plus the log  ( rhk ) regression . the bayes factor relative to model ( b ) is 1.3  10 : 1.0 .
final parameter estimates for the planetary candidates were derived from this three ak signals plus three keplerian signals plus log  ( rhk ) regression model .
the last row gives the date of the periastron passage that occurs just prior to our reference date jdb = 55855.6693  d which is the unweighted mean observation time for rv 3 .
rv 3 parameter estimates and 68  per  cent credible boundaries for three planetary candidates discussed in the text .
the quoted eccentricity values are for the mode while for the other parameters we quote the median .
the errors on semi - major axis , a , and msin  i do not include the uncertainty in the mass of the host star .
table  9 shows the true p , m , k , e parameters of the seven planetary signals that were employed in the rv 3 simulation .
the true p , m , k , e parameters of the seven planetary signals that were employed in the rv 3 simulation . clearly , two planetary signals with k > 1 were recovered in our analysis and the p/2 harmonic of the p = 2315  d signal .
there is evidence for a signal at p   590  d in the rhk - corrected data ( fig . 
it does not survive the removal of the larger 1200  d signal , possibly because it is essential a harmonic of the 1200  d signal .
24 shows the residuals , their gls and differential gls periodogram , for a model consisting of three ak signals plus three keplerian signals plus log  ( rhk ) regression .
the highest residual is at p = 1.036d and has a p - value 0.001 .
weak features are in evidence close to p = 1.1188 and 26.3  d but nothing close to 201  d. the top - left panel shows the rv 3 residuals for a model consisting of three ak signals , three keplerian signals plus the log  ( rhk ) regression .
the other panels show the gls periodogram ( top right ) and the differential periodogram for selected period ranges . a portion of the differential gls periodogram illustrating the 1200  d signal for the rv 3 data after the 12.5 , 17 , 48.8 and 333  d signals plus the log  ( rhk ) regression were removed .
figure for rv 3 data shows the eccentricity versus period parameters for the six signals obtained using a model consisting of three ak signals , three keplerian signals plus log  ( rhk ) regression .
comparison of the raw rv and fwhm and log  ( rhk ) diagnostics indicated that there was a clear correlation between the three .
the raw rv 4 data had a standard deviation of 8.30  m  s. after removing the best linear regression fit with log  ( rhk ) as the independent variable , the standard deviation was reduced to 3.32  m  s. the two rows of fig .  25 show that the rhk - corrected rv data and fwhm control together with their gls periodogram on the right .
26 shows the differential periodogram for the rhk - corrected rv data for selected period ranges .
the rv 4 data and fwhm ( control ) after removing the best linear regression fits to log  ( rhk ) together with their gls periodograms on the right .
the differential gls periodogram for the rhk - corrected rv 4 data for selected period ranges .
the lower - left panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) . the map apodization window for the p
= 11.75 signal clearly spans the duration of the data and on this basis alone it would be classified as a planetary candidate .
28 ) has a p( ) = 0.06 at the 11.75  d period , which is the main reason it was classified as ( p ? ) .
also , the proximity of the sa signals at 11.33 suggested that the 11.75  d signal might be another sa component .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the eight - signal ak periodogram of the rv 4 data .
the upper right shows log10[prior  likelihood ] versus period showing the eight periods detected .
the middle left shows the values of the eight unknown period parameters versus iteration number .
the lower left shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
the rv 4 residuals , their gls periodogram and differential periodogram for the eight ak signal model . based on the apodization window
, the p = 0.943  d signal was also initially classified as a possibly planetary ( p ? ) .
bayesian model comparison ( see below ) indicates that the p = 0.943  d signal is not significant .
the weighted rms residual = 1.52  m  s for the eight ak signal model .
the standard deviation of the extra gaussian white noise term s = 1.36  m  s. the mean measurement uncertainty was 0.674  m  s. table 10 lists ( in order of extraction ) the signals , their nominal parameter values and designation of signal type as planetary ( p ) , ( p ? ) or stellar activity ( sa ) based on the mcmc parameter estimates .
( p = 11.75 and 0.943  d ) which as mentioned above were initially classified as possibly planetary candidates . order of signals extracted from rv 4 together with approximate parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) . to test the significance of the possible p = 0.943  d planetary candidate
, we computed the bayes factor comparing the following models : ( a ) six ak signals plus the log  ( rhk ) regression , and ( b ) six ak signals plus one keplerian ( 0.943  d ) plus the log  ( rhk ) regression .
the bayes factor of ( a ) relative to model ( b ) is 64 : 1.0 which clearly indicates that the 0.943  d signal is not significant .
that left one possible planetary candidate at p = 11.75  d. when the results of the challenge were announced , it turned out that rv 4 data set contained no planetary signal .
comparison of the raw rv and fwhm and log  ( rhk ) diagnostics indicated that there was a strong correlation between the three .
the raw rv 5 data had a standard deviation of 8.92  m  s. after removing the best linear regression fit with log  ( rhk ) as the independent variable , the standard deviation was reduced to 2.59  m  s. the two rows of fig . 
29 show the rhk - corrected rv data and fwhm control together with their gls periodogram on the right .
30 shows the differential periodogram for the rv data ( rhk corrected ) for selected period ranges .
the rv 5 data and fwhm ( control ) after removing the best linear regression fits to log  ( rhk ) together with their gls periodograms on the right .
the differential gls periodogram for the rhk - corrected rv 5 data for selected period ranges .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the six - signal ak periodogram of the rv 5 data .
the upper right shows log10[prior  likelihood ] versus period indicating multiple six period solutions .
the middle - left panel which plots period versus mcmc iteration indicates that the 0.9613 , 15.3 , 21.2 , 40 , 179.7 , 335  d solution is dominant .
the lower left shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
the upper - left panel is a plot of the log10[prior  likelihood ] versus iteration for the six - signal ak periodogram of the rv 5 data .
the upper right shows log10[prior  likelihood ] versus period indicating multiple six period solutions .
the lower left shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is mainly hidden below the grey ) .
the lower right is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
table 11 lists ( in order of extraction ) the signals , their nominal parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) based on the mcmc parameter estimates .
the weighted rms residual = 1.18  m  s for the six ak signal model .
the standard deviation of the extra gaussian white noise term s = 0.98  m  s. only one signal ( p = 0.9610  d ) was initially classified as a possible planet although an examination of the apodization window span ( lower - left panel in fig . 
order of signals extracted from rv 5 together with approximate parameter values and designation of signal type as planetary ( p ) or stellar activity ( sa ) . to test the significance of a possible p = 0.9610  d planetary candidate
, we computed the bayes factor comparing the following models : ( a ) five ak signals plus the log  ( rhk ) regression , and ( b ) five ak signals plus one keplerian ( 0.961  d ) plus the log  ( rhk ) regression .
the periods of the five apodized planets are the first five listed in table 11 .
the bayes factor of ( a ) relative to model ( b ) is 7.9  10 : 1.0 which clearly indicates that the 0.961  d signal is not significant . on the basis of this analysis
, we can not make a good case for any planetary signal in this data set .
what can we learn looking back at our analysis from knowledge of the true planetary signals used in the simulation ?
table 12 shows the true p , m , k , e parameters of the six planetary signals that were employed in the rv 5 simulation .
the 0.961  d signal is consistent with a sidereal alias of 26.4237  d which is close to the value of the true 26.2  d period .
if our signal at 179  d is a consequence of the true 173  d signal , then it is apparent that for a true amplitude of 0.59  m  s the apodizing model is failing to discern this as a planetary signal .
the weighted rms residual = 1.18  m  s of the six ak signal model while the mean measurement uncertainty is 0.674  m  s. the true p , m , k , e parameters of the six planetary signals that were employed in the rv 5 simulation .
we now turn to consider how successful the ak method has been in the five challenge simulations that were analysed in this way . out of a total of 10 planetary signals with k > 1  m 
s , the method correctly identified six as class ( p ) , two as ( p ? ) and claimed another as ( p ) which turned out to be the first harmonic of a true planet with period longer than the data set .
it missed a true planet with p = 596  d and k = 1.91  m  s. the ak method did not successfully detect any of the seven planets with k < 1  m  s. although the mean measurement error for all five of the data sets was 0.674  m  s , the fit residuals ranged from 1.2 to 1.7  m  s. on the other hand , the ak method did not lead to any false detections .
two signals were initially listed as ( p ? ) in data set 4 but one was strongly ruled out by the bayes factor . for the other , the differential gls periodogram played an important role in restricting the classification to ( p ? ) .
another signal was initially listed as ( p ? ) in data set 5 but again was strongly ruled out by the bayes factor . for the challenge data sets analysed to date , the standard deviation of the data is clearly dominated by sa .
table 13 summarizes the performance of the ak approach regarding its ability to penetrate the fog of sa .
the last column gives the ratio of the initial data standard deviation to the standard deviation of the final residuals which has an average value of 5.9 ; however , the residuals are on average a factor of 2.3 the mean measurement uncertainty .
the average reduction that is due to the linear regression term alone is a factor of 2.5 .
as was mentioned in section 3 , the ak method was employed to distinguish between sa and planetary candidates .
32 shows the individual k parameter marginal distributions for the eight planets detected plus one test data set planet based on ak model fits plus a log  ( rhk ) regression term .
set , the results are based on model fits consisting of m aks plus a log  ( rhk ) regression , where m is the total number of signals fit .
33 shows the individual k parameter marginal distributions for the eight planets detected plus the one test data set planet , based on model fits consisting of k aks plus m  k keplerians plus a log  ( rhk ) regression .
m is the total number of signals and ( m  k ) is the number classified as a planetary candidate in the analysis . in this case
set , the results are based on model fits consisting of k aks plus ( m  k ) keplerians plus a log  ( rhk ) regression , where m is the total number of signals .
there is a known bias towards higher eccentricity values for keplerian models of rv data , e.g. shen & turner ( 2008 ) and zakamska , pan & ford ( 2011 ) .
is it really true that the ak models are biased towards higher k values , or do they have a bias towards higher e , which would also inflate the k values ?
34 show a comparison of the marginal eccentricity distributions for the three planets detected in the rv 1 data set .
the solid curves are the distributions obtained with the ak model and the dashed are for the pure kepler model . if anything , the eccentricity bias of the ak model fits is towards lower values when compared to the pure kepler fits .
the lower row of panels show the corresponding comparison for k values for the rv 1 planets .
the top row of panels show a comparison of the marginal eccentricity distributions for the three planets detected in the rv 1 data set .
the solid curves are the distributions obtained with the ak model and the dashed are for the pure kepler model .
figs 35 and 36 show similar plots for the period parameter . in fig . 
36 , only two of the nine marginals overlap the true period within the 68  per  cent credible region .
this could be a consequence of the remaining undetected planets with smaller k values and/or inadequate modelling of the sa signals , leading to correlated residuals .
set , the results are based on model fits consisting of m aks plus a log  ( rhk ) regression , where m is the total number of signals fit .
set , the results are based on model fits consisting of k aks plus ( m  k ) keplerians plus a log  ( rhk ) regression , where m is the total number of signals .
the rv challenge winning team of mikko tuomi and anglada escude included a first - order moving average ( ma ) in their keplerian analysis along with a correlation with all sa observables .
it therefore seemed logical to explore whether augmenting our ak models with a ma term would lead to a noticeable improvement in results .
to investigate we did several runs on the rv 3 data set that included a ma term ( tuomi et  al .
our general pth order ma model rv is given by z(ti ) as follow : 
 ( 4)\documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } \begin{eqnarray * } z(t_i ) & = & y(t_i ) + \sum _ { \eta = 1}^p \big [ \gamma _ { \eta } \times{\rm exp } \big ( -\frac{(t_i - t_{i-\eta } ) } { \lambda _ { \eta } } \big ) \nonumber \\ & & \times ( v(t_{i-\eta } ) - y(t_{i-\eta } ) ) \big ] , \end{eqnarray*}\end{document}where y(ti ) is given by equation ( 1 ) , v(ti ) is the measured radial velocity , p is the order of the moving average , and  and  are the unknown ma parameters for  = 1 to p. of the three ma models , the second - order ma was judged to yield the best results . the upper panel of fig . 
37 shows the eccentricity versus period parameters for an rv 3 six - signal ak model plus log  ( rhk ) regression which includes a second - order ma term .
it exhibits a broad low - eccentricity peak centred on p = 2300  d. recall that the true period is 2315  d. the other periods are similar to our earlier results with no ma term .
the lower panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is partially hidden below the grey ) .
although the apodization window for the p = 2300  d signal spans the data , the coverage for the p = 17 and 48  d is judged to be not as good as that obtained for the no ma case . the middle is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
the upper panel shows the eccentricity versus period parameters for an rv 3 six - signal ak model plus the log  ( rhk ) regression and includes a second - order ma term .
the middle panel shows a plot of the apodization time constant ,  , versus apodization window centre time , ta , where the symbols refer to different signal periods specified in the top panel .
the lower panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is partially hidden below the grey ) .
38 shows the eccentricity versus period parameters for an rv 3 six - signal model consisting of three ak signals , three keplerian signals plus log  ( rhk ) regression and includes a second - order ma term .
the lower panel shows the span of the apodization window within the overall data window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is partially hidden below the grey ) .
the middle panel is a plot of the apodization time constant ,  , versus apodization window centre time , ta .
the values of the ma parameters and the extra noise parameter are \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \gamma _ 1 = 0.59 \pm 0.05 , \lambda _ 1 = 3.0_{-1.4}^{+0.9}$\end{document }  d , \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{upgreek } 
 \usepackage{mathrsfs } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document } 
 } { } $ \gamma _ 2 = 0.35 \pm 0.05 , \lambda _ 2 = 68_{-63}^{39}$\end{document }  d and s = 1.39  0.05  m  s. finally , we show a comparison of k and p parameter marginal distributions for the no ma case ( black ) and the second - order ma case ( grey ) .
39 show k parameter marginal distributions for the 17 ( left ) and 48  d ( right ) planets detected in rv 3 .
set , the results are based on model fits consisting of k aks plus ( m  k ) keplerians plus a log  ( rhk ) regression , where m is the total number of signals and ( m  k ) is the number classified as a planetary candidate .
the upper panel shows the eccentricity versus period parameters for an rv 3 six - signal model consisting of three ak signals , three keplerian signals plus the log  ( rhk ) regression and includes a second - order ma term .
the middle panel shows a plot of the apodization time constant ,  , versus apodization window centre time , ta , where the symbols refer to different signal periods specified in the top panel .
the lower panel shows the apodization window for each signal ( grey trace for map values of  and ta , black for a representative set of samples which is partially hidden below the grey ) .
the upper panels show k parameter marginal distributions for the 17 and 48  d planets detected in rv3 .
the black trace is for the no ma case and the grey corresponds to the second - order ma case .
set , the results are based on model fits consisting of k aks plus ( m  k ) keplerians plus a log  ( rhk ) regression , where m is the total number of signals .
the second - order ma analysis was repeated using a wider prior search range for  extending from 10  d to eight times the data duration .
this resulted in a change in period of one of the sa signals from 8.3 to 29  d ( with  = 10  d ) and a change in the parameters of the 17  d signal such that apodization window no longer spanned the data duration . based on this very limited comparison , it is not clear that including a ma term has improved our ability to distinguish sa and planetary signals in this one particularly challenging data set .
for the rv 3 data set , the ma analysis suggests that there is a significant correlation in the residuals of the six ak signal fit .
the results reported in this paper indicate that an ak model provides a useful way to distinguish planetary signals from sa - induced signals in precision rv data . the general model for m ak signals
includes a linear regression term between rv and the sa diagnostic log  ( rhk ) , as well as an extra gaussian noise term with unknown standard deviation . in the current implementation ,
the ak method achieved a reduction in sa noise by a factor of approximately 6 .
the analysis also employed a differential version of the gls periodogram that uses a control diagnostic to provide an additional means of distinguishing sa signals and to help guide the choice of new periods .
the ak model estimates for k were found to have a bias to larger values .
for this reason , final parameter estimates for the planetary candidates were derived from fits that include ak signals to model the sa components and simple keplerians to model the planetary candidates .
preliminary results are also reported for ak models augmented by a ma component that allows for correlations in the residuals .
based on a very limited comparison , it is not clear that including a ma term has improved our ability to distinguish sa and planetary signals .
the final residuals are still a factor of 2.3 higher on average than the measurement uncertainties .
a simple symmetrical gaussian ak is only expected to provide an approximate model for some of the sa signals .
 a new apodized keplerian ( ak ) model is proposed for the analysis of precision radial velocity ( rv ) data to model both planetary and stellar activity ( sa ) induced rv signals . 
 a symmetrical gaussian apodization function with unknown width and centre can distinguish planetary signals from sa signals on the basis of the span of the apodization window . 
 the general model for m ak signals includes a linear regression term between rv and the sa diagnostic log  ( rhk ) , as well as an extra gaussian noise term with unknown standard deviation . 
 the model parameters are explored using a bayesian fusion markov chain monte carlo code . a differential version of the generalized lomb  scargle periodogram that employs a control diagnostic provides an additional way of distinguishing sa signals and helps guide the choice of new periods . 
 results are reported for a recent international rv blind challenge which included multiple state - of - the - art simulated data sets supported by a variety of sa diagnostics . in the current implementation , 
 the ak method achieved a reduction in sa noise by a factor of approximately 6 . 
 final parameter estimates for the planetary candidates are derived from fits that include ak signals to model the sa components and simple keplerians to model the planetary candidates . 
 preliminary results are also reported for ak models augmented by a moving average component that allows for correlations in the residuals . 
INTRODUCTION
METHODOLOGY
RESULTS
DISCUSSION
CONCLUSIONS