Kepler Presearch Data Conditioning II - A Bayesian Approach to Systematic Error Correction

It is standard practice when removing systematic errors in stellar data to use robust least-squares (LS) on a set of basis vectors, as is used in methods such as SYSREM (Tamuz 2005) and TFA (Kovacs 2005). A robust LS approach, as outlined below in § 2.2, can find a chance linear combination of the systematic error model components that reduces the bulk rms at the expense of distorting the intrinsic stellar variations and introducing additional noise on short timescales. The fundamental problem with this approach is the fact that the implicit model fitted to the data for each star is incomplete. LS cotrending projects the data vector onto the selected basis vectors and removes the components that are parallel to any linear combination of the basis vectors. This process is guaranteed to reduce the bulk rms residuals, but may do so at the cost of injecting additional noise or distortion into the flux time series. Indeed, this occurs frequently for stars with high intrinsic variability, such as RR Lyrae stars, eclipsing binaries, and classical pulsators. For example, if one of the model terms is strongly related to focus variations, and the long-term trend is for the width of the stellar point spread function (PSF) to broaden over the observation interval, then the flux for all stars should decrease over time. A LS fit, however, may invert the focus-related model term for a star whose flux increases over the observation interval, thereby removing the signature of intrinsic stellar variability from this light curve because there is a coincidental correlation between the observed change in flux and the observed change in focus. Given that the star would be expected to dim slightly over time, if anything, due to the focus change, PDC should be correcting the star so that it brightens slightly more than the original flux time series would indicate.

The situation is analogous to opening a jigsaw puzzle box and finding only 30% of the pieces present. LS gamely tries to put the jigsaw pieces together in order to match the picture on the box cover by stretching, rotating, and translating the pieces that were present in the box. The result is a set of pieces that roughly overlap the picture on the box cover, but one where the details do not necessary match up well, even though individual pieces may obviously fit. In order to improve the performance of robust LS, we need to provide the fitter with constraints on the magnitudes and signs of the fit coefficients. These constraints can be obtained by using the ensemble behavior of the stars to develop an empirical model of the underlying physics. For example, the photometric change that can be induced by a pointing change of 0.1" must be bounded, and this bound can be estimated by looking at how the collection of stars behaves for a pointing change of this magnitude.

As an example of this analysis, and to demonstrate systematic trends in the Kepler data, take channel 2.1, which is the most thermally sensitive CCD channel in Kepler's focal plane. Nearly all stars on this channel exhibit obvious focus-and pointing-related instrumental signatures in their pixel time series and flux time series. Figure 1 shows several characteristic light curves for typical targets¹⁰ on channel 2.1 during the Kepler Quarter 7 data season normalized by the median flux value. Note the long-term increase for all flux curves over the 90-day interval. This is due to seasonal changes in the shape of the telescope, and therefore its focus, as the Sun rotates about the barrel of the telescope while the spacecraft orbits the Sun and maintains its attitude fixed on the field of view (FOV). All light curves exhibit these long-term trends, but to differing degrees. Also present are some short-term oscillations evident but mainly obscured in variable targets, which are due to focus changes driven by a heater cycling on and off to condition the temperature of the box containing reaction wheels 3 and 4 on the spacecraft bus. This component was receiving more and more shade throughout this time interval, and the thermostat actuated more frequently over time. A target that is varying on levels and periods similar to these systematic effects can obscure the systematics, making identification difficult. Looking at a single quiet target shown in Figure 2 (one of the targets shown in solid blue inFig. 1) allows us to more clearly see the systematic trends which are exhibited in all the targets in Figure 1, but obscured by variability. Each Earth point, one at the beginning of the quarter (cadence index 0) and after each monthly downlink (cadences 1500 and 2800) results in a heating of different sides of the telescope as the spacecraft reorients the antennae to downlink data. The Earth points themselves are gaps in the data. They result in periods of local heating and cooling, distorting the telescope. A characteristic recovery time is also evident. The other trends as described above are also clearly evident. A final short data gap is also evident at cadence 3950, but this was not due to a reorientation of the spacecraft, so no thermal recovery is present.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a counter example, Figure 3 shows the same highly variable target in solid blackin Figure 1. Note how this highly variable star almost completely obscures the long-term trend. The targets shown in Figures 2 and3 will hereby be referred to as the "quiet target" and the "variable target" and used as canonical example targets in § 3.

Zoom In Zoom Out Reset image size

How can we separate intrinsic stellar variability from instrumental signatures? We do not expect intrinsic stellar variability to be correlated from target to target, except for rare coincidences, and even then one would not expect a high degree of correlation for all time scales. However, we do expect instrumental signatures to be highly correlated from target to target and can exploit this observation to provide constraints on the cotrending that PDC performs.

Figure 4 shows a histogram of the absolute value of the correlation coefficient for 1864 targets on channel 2.1. The targets' light curves are highly correlated, as evidenced by the near complete pile-up near an absolute correlation coefficient of 1. Examination of individual light curves indicates that these light curves are contaminated to a large degree by instrumental signatures, as evidenced in Figures 1 and 2. However, not all the targets are dominated by systematic errors. The trick is to come up with a method that can distinguish between intrinsic stellar variability and chance correlations with linear combinations of the diagnostic time series used to cotrend out systematic errors.

Zoom In Zoom Out Reset image size

A Bayesian approach called the maximum a posteriori (MAP) method allows us to provide PDC with constraints on the fitted coefficients to help prevent over-fitting and distortion of intrinsic stellar variability. The application of a Bayesian MAP approach to Kepler data was first proposed in (Jenkins 2012). In this exposition we follow the notation of (Kay 1993).

The PDC-MAP technique examines the behavior of the robust LS fit coefficients across an ensemble of targets on each CCD readout channel in order to develop a description for the "typical" value for each model term. This description is a probability density function (PDF) that can be used to constrain the coefficients fitted in a second pass. To develop this approach, we build on a maximum likelihood approach.

The maximum likelihood approach models each light curve, y, as a linear combination of instrumental systematic vectors, referred to as cotrending basis vectors or CBVs, arranged as the columns of a design matrix, H, plus zero-mean, Gaussian observation noise, w:

The maximum likelihood estimator (MLE) seeks to find the solution, , that maximizes the likelihood function, p(y; θ ), given by

where C_w is the covariance of w and N is the number of data points. Taking the gradient of the log of equation (2), setting it equal to zero, and solving for θ yields the familiar LS solution,

This solution assumes the model H is a complete model to the data. We will show that the Bayesian model accounts for an incomplete model, which is the common case when removing systematics from stellar signals.

Adopting the Bayesian approach allows us to incorporate side information, such as knowledge of prior constraints on the model, in a natural way. Bayesianists view the underlying model as being drawn from a distribution and the data as being one realization of this process. In this case, we wish to find the MAP estimator of the model coefficients, given the observations (data):

where we have applied Bayes' rule (D'Agostini 2003) to simplify the expression. In this equation, p(θ ) is the prior PDF of the model coefficients. The mathematical form for p(y|θ ) is the same as for the non-Bayesian likelihood function p(y; θ ) in equation (2).

For illustration purposes, if we adopt a Gaussian form for the coefficient distribution, θ , then p(θ ) takes a closed form solution,

where C_θ and μ_θ are the covariance and mean of θ , respectively, and we assume that the coefficients are uncorrelated (which will hold true for orthogonal basis functions). We can then maximize equation (4), using equation (5), by maximizing its log likelihood,

Taking the gradient of equation (6) with respect to θ , setting it to zero, and solving for θ yields

If the observation noise, w, is zero-mean, white Gaussian noise with variance σ², then equation (7) can be rewritten as

The key to this Bayesian technique is to determine when to preference, or weight, the prior PDF over the conditional PDF. If the variance in the data is large compared to the "spread" allowed by the prior PDF for the model, then the MAP estimator gives more weight to the prior so that as σ² → ∞. This case would correspond, for example, to targets with large stellar variability, such as with the target given in Figure 3. In this case, the MAP weighting constrains the fitter from distorting the light curve and introducing noise on a short time scale. Conversely, if the variance in the data is small compared to the degree to which the prior PDF confines the model, the MAP estimator "trusts" the data over the prior knowledge, and as σ² → 0. This case would correspond to targets with small stellar variability, such as with the target in Figure 2 where there is little risk of over-fitting and distortion of the light curves, and it is a "safe" bet to use the conditional, LS fit.

The cotrending basis vectors are obtained using singular value decomposition. In order to have equal representation for all light curves independent of their absolute magnitude, we first normalize the targets by their median flux values [Δflux/median(flux)]. We then select the 50% most highly correlated targets based on the median absolute Pearson correlation. This cut generates a set that exhibits the strongest trends in the data. It mostly removes targets with large variability, but not completely. A variable star exhibiting a strong trend can still remain in the reduced list. We therefore first make a cut on the estimated variability of each target.

An estimate of the intrinsic stellar variability of each target must be found. Herein lies the fundamental chicken-and-egg problem of the cotrending method. We need to know the stellar variability of each target in order to know how much to rely on the prior. But if we already knew the stellar variability, then we would have no need for the prior—the cotrending solution would simply be the intrinsic stellar variability subtracted from the light curve, minus a Gaussian noise estimate. This issue is not specific to this particular cotrending method, either. Whenever a system is characterized with an incomplete model, there exists the problem of identifying the components in the system not represented in the model. We fortunately do not need to absolutely know the variability, we only need an estimated metric in order to weigh the prior. This estimate can be obtained by comparing a third-order polynomial to the light curve. The polynomial will remove any long-term trends leaving behind a roughly detrended curve. The standard deviation of this polynomial-removed light curve results in a rough calculation of the variability of the target. Removing a low-order polynomial is essentially a high-pass filter, we are therefore assuming any long-term trends are systematic and short-term trends are stellar. There are numerous counter-examples of short-term trends that are actually systematic—reaction wheel zero crossings is a good example. However, short-term systematic trends tend to be small in magnitude, whereas long-term systematics tend to result in large diversions in the flux amplitude. Likewise, there are examples of intrinsic long-term trends, but they are generally smaller than the systematic trends. Since we are only concerned with the relative amplitude of stellar versus systematic variability, we are using the low-pass filter to distinguish two characteristic realms of influence: long-term trends dominated by systematics and short-term trends dominated by intrinsic stellar variation. An example is shown in Figure 6. Here, a highly variable target is compounded with a long-term DVA and thermal trend. For periods less than 400 cadences, the variance in the flux is dominated by stellar features. The long-term variance, and the general trend to higher flux values, is due to systematics. The variance of the residual after removing the polynomial fit, labeled as "coarsely detrended light curve" in the figure, gives a rough estimate of the stellar variability of this target. Note that there are still systematic features in the detrended light curve. They are, however, small in magnitude compared with the stellar variability.¹¹

Zoom In Zoom Out Reset image size

The variability, V, is measured using

where is the standard deviation of the third-order polynomial detrended light curve, Δy is the uncertainty of the flux data as determined by the PA pipeline component (Clarke 2010), and is the median variability over all light curves in the sample. The normalization by the uncertainty is to ensure the noise in the data is not included in the stellar variability. The normalization by the median variability is so that a variability of 1 is considered typical, thereby simplifying the analysis parameterization. Figure 7 shows a histogram of the measured variability for all targets on channel 2.1. The median of all of these values is evidently 1, and the distribution is typical for all channels where most are close to typical variability but with a long tail to high variability (note the log scale for the x-axis). There are two cutoff thresholds plotted, as well. The upper (dashed red ) is the threshold to determine if a target is "highly variable." The lower (solid green ) is to determine if a target is "very quiet." The very quiet targets have such a low amount of variability that using the prior PDF when generating the fit has been found to be problematic. Any targets above the high variability threshold are removed from the reduced list. The remaining targets are sorted with respect to median absolute correlation, and the 50% most highly correlated are used for SVD.

Zoom In Zoom Out Reset image size

Due to all targets being normalized by their median, most targets pass through zero amplitude at the midpoint, which can be seen as a "node" in the light curves at cadence 2200 in Figure 1. If SVD was performed on this set "as is", then all the strong cotrending basis vectors would have zero amplitude at the midpoint. The basis vectors would therefore be unable to remove systematics in the minority of targets that do not pass through zero at the midpoint. The light curves are therefore "dithered" slightly by a zero-mean Gaussian dithering magnitude in order to slightly "spread" the light curves about the zero-flux value. Since the dithering is zero-mean, this has no effect on the resultant basis vectors other than to remove the artificial zero crossing node at the midpoint. Note that the dithering is only used to generate the basis vectors. The cotrending is performed on the non-dithered, but still median-normalized, light curves.

Figure 8 shows the singular values from the singular value decomposition. This figure is characteristic of all channels, showing two or three strong singular values, then a slowly tapering region for about another dozen values, until finally asymptotically approaching zero (as is expected with SVD). The first several left-singular vectors are selected (typically the first eight) to become the cotrending basis vectors. These first singular vectors exhibit the principle trends in the data due to DVA, pointing errors, impulses due to Argabrightenings (Witteborn 2011), focus errors, and reaction wheel zero crossings, among other trends. The number of basis vectors used is generally eight, however a SNR test is performed where the S/N is determined by

with A_signal and A_noise being the rms of the light curve and noise floor, respectively. The noise floor is approximated by the first differences between adjacent flux values. Any of the eight basis vectors with a S/N below a threshold of 5 dB are removed, but only a small number of basis vectors over the entire FOV are removed by the S/N test. Most have high S/N. There are other sophisticated methods to find the dimensionality of an eigensystem, such as Bayesian model selection (Minka 2000). These are not used because they tend to pick too high a dimensionality in this particular situation. We wish to find only the singular vectors with systematics; the lesser singular vectors do contain light-curve signal information, but not necessarily systematics, and we have found including them in the MAP fit adds no value yet slows down the algorithm.

Zoom In Zoom Out Reset image size

For a minority of basis vectors, a few target light curves can dominate the signal. The normalization process attempts to "equalize" the strength of all targets, but a small number of light curves can be over-represented in the singular vectors from SVD. To eliminate this, we calculate an entropy metric for each basis vector using the following entropy calculation

where p(x) is a PDF created from the right-singular vectors from SVD (referred to as the " V-matrix"),

The V-matrix contains the contribution of the signal in the basis vector from each target light curve. We must first normalize the entropy calculation to a Gaussian distribution, which has the highest entropy of any continuous distribution with the same 2nd moment. The entropy of a Gaussian is

σ being the second central moment of V_ki for fixed i. The resultant relative entropy is therefore

If one (or a few) targets are domineering, then they will have much larger values in the V-matrix then all the other targets. A negative value of the entropy calculation will identify this condition. Bad entropy is somewhat arbitrary, but we have found that a value below -0.7 is poor. For any basis vectors with identified poor entropy, the V-matrix column for that basis vector is examined for stand-out targets. The offending targets are removed, and SVD is re-computed on the remaining targets. The process is iterated until the entropy of all basis vectors is below -0.7. Typically, no more than a couple iterations is necessary, and fewer than 20 targets are removed (out of 2500 total targets).

Figure 9 shows the first eight cotrending basis vectors generated for channel 2.1, and Figure 10 shows just the first basis vector. Trends can be found in all the vectors, but it is useful to concentrate on the first and strongest. Here the most characteristic trends and systematics in the data can be found. The general trend to higher flux is due to the seasonal change and solar orientation. The short recovery periods at cadence indices 0, 1500 and 2800 are due to monthly downlinks. The short spikes at 700 and 1450 are due to artifacts from correcting cosmic rays near reaction wheel zero crossing periods.¹² The periodic oscillation is due to heater cycling. Notice how the basis vector in Figure 10 closely follows the flux light curve in Figure 2. This signifies that virtually all the features in this flux light curve are due to systematic effects and not intrinsic stellar variability. In theory, any features in the light curves in Figure 1 that are not represented in the basis vectors in Figure 9 are intrinsic to the target. However, a simple LS projection of the light curves on the basis vectors will not produce desirable results for all targets, as will be shown below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once the cotrending basis vectors are found, a robust LS fit is performed on each target. This creates the empirical data used to generate the prior PDF.

The prior PDF is based on the distribution of robust fit coefficients of the basis vectors for all light curves using the method described in (Bowman 1987). This method computes a probability density estimate of the sample data based on a normal kernel function using a window that is a function of the number of points in the data sample. The form of the prior PDF will depend on the parameterization of the robust fit coefficients. We must thus decide how to parameterize the coefficients to best extract the correlations. Some systematic effects are caused by focal plane irregularities and instrumental vibrations which are stronger near the edges of the CCD frame (nearer to the spacecraft housing). There are also other issues that are dependent on the physical position of each pixel on the CCD. Therefore, the targets' locations in the sky as characterized by right ascension and declination are reasonable parameters to characterize target location with respect to the sources of systematic effects. The targets' influence by systematic effects is also directly related to the stellar magnitude, since different magnitude targets result in different saturation levels of the CCD pixels. For example, the readout electronics for the CCDs are sensitive to temperature drift, but the sensitivity is non-linear with respect to CCD flux levels. So brighter targets are affected by instrument temperature differently than dim targets. We therefore parameterize the prior PDF with three independent variables: (1) stellar magnitude (K_p), (2) right ascension (R.A.), and (3) declination (decl.).

Figures 11 13 show the robust fit coefficients for a basis vector plotted against K_p, R.A. and declination. The blue star data is for all targets, whereas, the redcircle data is just for those targets remaining for SVD after the cuts discussed in § 3.1. The solid blueand dashed redcurves in Figures 12 and13 are the travelling window means of the bluestar and redcircle data, respectively. The cuts clearly produce a bimodal distribution in K_p for this basis vector. A simple Gaussian fit would not reproduce this and demonstrates that the systematic trends are correlated with K_p, but variable targets are masking the true correlation when a simple robust fit is performed. The correlations in R.A. and decl. are also evident but to a lesser extent. Notice also that the mean (solid bluecurves) are biased compared to the dashed redcurve. This is again because the variable targets are masking the true trends in the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some basis vectors exhibit stronger trends in K_p, R.A., or decl., but not necessarily all three simultaneously, as is expected if the different systematics represented by the basis vectors have different instrumental sources. Plotting different basis vectors and/or channels reveals different trends and correlations.

We want to mainly rely on targets in the neighborhood around the target we are fitting, referred to as the target under study (TUS), in R.A.-, decl.-and K_p-space when generating the prior PDF. If we simply found an evenly weighted PDF, then a large cluster of targets with a certain coefficient value, even if nonlocal to the TUS, would always dominate the peak of the prior PDF. We therefore use a weighted probability density estimate based on the standardized Euclidean distance between targets x and y,

where Λ is a diagonal matrix whose diagonal elements give the relative weighting for each dimension. A straight normalization in each dimension by its standard deviation would result in equal weighting of all three dimensions, but we wish to overemphasize the prior PDF in dimensions that exhibit greater correlations, and in our case the robust fit coefficients exhibit a stronger correlation in K_p than in R.A. or declination. The Λ matrix diagonals are therefore

where mad(θ_i) is the median absolute deviation of the coefficient distribution along dimension i and S_i is the scaling factor for dimension i, S_i = {2, if i⇒K_p or 1,otherwise.} The above weighting results in the K_p dimension weighted twice as much as R.A. and decl. when generating the prior PDF. That is, targets further away in the K_p dimension are weighted proportionately less than in R.A. and declination. This effectively results in taking a tighter cut in K_p-space to emphasize the greater correlation in that dimension. Since the PDF is weighted by this distance metric, the PDF will emphasize the correlation in K_p and yet still be sensitive to the trends in R.A. and declination. The median absolute deviation is used instead of the standard deviation in order to ignore outliers. The effect of the scaling factor S is not great compared with the overall distance weighting, however, it does improve results slightly. As development continues, this functionality may reveal itself to be more useful when the number of independent variables increases.

The weighting by equation (17) and how it affects the prior PDF is illustrated in Figures 14 16, the latter two being the prior PDFs for the same two targets in Figures 2 and3. The bluehistogram in all three figures is exactly the same, since they are generated from the same distribution of coefficients. However, the prior PDF (solid redcurve) is dramatically different. In Figure 14 a bimodal PDF is evident due to targets nearby to the TUS containing two clusters, around -1.3 and -0.85, and suggests that the coefficient value for the TUS should be one of these two values. Which value that is actually chosen will be dependent on the form of the conditional PDF and the weighting of the prior PDF as discussed in § 3.3. In Figure 16 the targets near the TUS have coefficients clustered around -0.34, which is far from the peak in the unweighted PDF. Using the unweighted PDF would have completely missed the actual systematic trend in the data near the TUS. The log of the prior PDF is plotted in these figures for direct comparison with equation (10), which results in a compression of the PDF near the top.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, the prior PDF is developed by generating a three-dimensional weighted distribution of robust LS fit coefficients in R.A.-, decl.-and K_p-space. This methodology makes no assumptions on the form of the PDF, Gaussian or otherwise, and allows PDC-MAP to identify and characterize the form of the systematic trends across the full distribution of targets.

For each light curve, the weighting parameter, W_pr in equation (10) is an empirical weighting parameter that is principally based on the variability of each target. The greater the variability, the greater we need to constrain the LS fit. However, there is another complication. A fit to the prior PDF is not always a good fit to the trend in the target. The reason for this disagreement is currently being investigated. One factor influencing the "goodness" of the prior fit is the sparseness of the targets in certain regions in R.A.-, decl.-and K_p-space. A sparse distribution will result in poor prior statistics. There are also other unknown causes resulting in poor priors for some targets, and so an additional parameter in the prior weighting is an evaluation of the "goodness" of the prior fit. The goodness is evaluated using a method similar to the variability calculation above. The prior fit is compared to a third-order polynomial fit to the light curve with a soft-wall cutoff using the following equation:

where G_raw is the "raw" goodness given by

and F pr and F poly are the prior PDF fit and the third order polynomial fit to the data, respectively. Normalization by the median absolute deviation (mad) of the polynomial fit-removed light curve allows for a comparison of the difference between the polynomial fit and the prior fit with respect to the variance of the target. The soft cutoff is to ensure that small changes in the light curve will not have dramatic changes in the weighting. The scaling parameter α_G is determined when the deviation of the prior fit to the polynomial fit becomes too poor to be useful in constraining the posterior fit and is tuned empirically by maximizing the PDC goodness metric which is described in a companion paper (Stumpe 2012). An example of a poor prior fit is given in Figure 17. Note how both the long-term trend and the Earth-point recoveries are much larger in the prior fit than in the actual data. Examples such as this are in the minority, but frequent enough to require the additional test for prior goodness. It could be proposed that this target is trending downward, cancelling out the upward trend of the prior fit. This is unfortunately not the case. Examination of Kepler season quarters 6 and 8 reveal that this target is not experiencing a general trend in quarter 7. The prior fit is indeed poor and should not be used to any large degree. The unfortunate side effect of the prior goodness test is that PDC-MAP is less sensitive to very long-term trends in the data. A true long-term trend in the data that cancels out the systematic trend can confuse the prior goodness metric which would interpret the fit as a bad prior. The only way to surely know the actual long-term trend is to examine multiquarter data. Future versions of the PDC module may indeed provide this functionality.

Zoom In Zoom Out Reset image size

The resultant full form to the prior weighting is

The parameters β_V and β_G being scaling factors for the variability and prior goodness, respectively, and are tuned to maximize the PDC goodness metric, like α_G above. Future work includes fully characterizing the conditions where the prior fit is poor, and thereby removing the unfortunate need for an empirical prior "goodness" metric.

In cases where the prior goodness is near zero, the fit reverts to a reduced robust fit where the number of basis vectors is limited to just the first several (default is four). A MAP fit has a pleasant feature where a large number of basis vectors can be used. The prior PDF restricts the fit from drifting drastically in the function space, searching the large set of basis vectors for a combination that reduces the bulk rms at the expense of distorting stellar features. If the prior cannot be used, then there is no such restriction, and the posterior PDF becomes a LS fit, so a more limited number of basis vectors must be used in order to constrain the fit. The first several basis vectors have very strong trends in most of the data, and also have low noise components, so they are generally safe to use, even with an unrestricted LS fit. It is also generally true that a target with a bad prior is so because the target is quiet, and any small deviation in the prior from a true trend is very noticeable, and the prior is neither necessary or desirable to use.

If the target is below the variability threshold shown in Figure 7, then the target is very quiet, and in many cases the use of the prior fit only worsens the fit over a LS fit and so the prior weighting is zeroed. This is due to the prior fit never being an exact match to the target trend, and even small deviations can "pull" the posterior fit away from a good fit. In such cases there is little risk of a quiet target biasing a LS fit away from a proper cotrending fit. The majority of targets do not fall into either of the above two cases, and the prior is used to the degree dictated by the prior weight, and a Bayesian MAP fit is performed.

Once the prior weighting is determined, the posterior PDF can be assembled using equation (10). Due to the empirical (and therefore non-analytical) form of the prior PDF, the posterior must be maximized numerically. In general, a multidimensional maximization is difficult and time consuming due to the risk of only finding local maxima. Fortunately, due to the use of SVD, the basis vectors are all orthogonal, so the various coefficients can be maximized sequentially. The process is therefore straightforward. The strongest singular vector is maximized first and then all subsequent singular vectors are maximized in turn.

Following along with the same two examples of a quiet and a variable target, Figures 18 and 19 show the final posterior PDF along with the prior and conditional PDFs. The blackdots and magentastars are the maxima of the prior and conditional PDFs, and the bluecircle is the maximum of the posterior PDF. Due to the varying scales of the three curves, the prior and conditional curves have been renormalized to the same scale as the posterior for illustration purposes. The conditional fit has the smooth quadratic form characteristic of a LS fit. The prior appears Gaussian on this scale but in general is not. The title of each plot gives the prior weight. For the quiet target, the variability is low so the prior weight is only 4.65, resulting in a minor correction to the conditional curve, and so the conditional and posterior curves are virtually identical. Also, in cases where the width of the prior PDF is large, the maximum is low and it makes little contribution to the posterior. A wide prior is equivalent to saying that the prior information results in little added information to a good fit, the extreme case being a flat prior PDF which provides no additional information. This is in contrast to the variable target where the conditional maximum is near the prior PDF, meaning the fit almost completely relies on the prior data, due to the high weight on the prior PDF of 976, resulting in only a slight influence of the conditional fit. In the case of the variable target in particular, if the prior PDF was not determined using a weighted distribution (in K_p-, R.A.- and decl.-space) then the maximum of the prior would have been at about 0.6, as shown by the unweighted histogram in Figure 16. This would have resulted in a poor fit to the systematics. The actual prior fit takes into account the location of the TUS and the systematic trends in targets nearby to the TUS.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Once the maximum of the posterior PDF is found for each basis vector, the MAP fit is a linear combination of the basis vectors. The resultant fits are in Figure 20 and 21. For the quiet target, all three fits roughly overlap the actual trend in the data. The prior fit is not an exact match, and the slight disagreement is to be expected since the prior is purely formulated using targets other than the TUS. For a quiet target such as this one, highly weighting the prior PDF would result in a degradation of the fit, and so instead the PDF relies mainly on the conditional PDF (i.e., the reddashed and greensolid curves overlap). The resultant light curve after the trend removal is in the bottom figure for both the MAP fit and the conditional fit, the later being a LS fit. For the quiet target, note how the resultant curve is near featureless above the noise floor. Some slight artifacts are not fully removed and methods to correct these are discussed in § 4. In the case of the variable target, the conditional PDF results in a fit that attempts to remove all features in the light curve, whereas the prior PDF correctly identifies just the systematic trends in the data. In this case, it is beneficial to rely principally on the prior PDF. The prior cannot be well discerned in the figure because it lies under the MAP fit. The conditional fit also introduces a considerable amount of noise into the corrected light curve, due to it being constructed more from the lesser basis vectors (shown in Fig. 9) which contain larger noise components.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Propagation of uncertainties is not necessarily straightforward because a covariance matrix is difficult to formulate for an empirical prior PDF. As a first approximation, the propagation can be assumed to be through a LS solution—which is close to the solution for most targets. If C_raw and C_cot denote the covariance matrices for the temporal samples of the raw and cotrended flux time series for a given target, then the uncertainties may be propagated (disregarding the uncertainty in the mean level which can be considered to be negligible) by

where the transformation T_cot is defined by

H being the same design matrix as in equation (10). This is overly conservative since the posterior PDF is more constrained than a simple LS fit. A more accurate propagation of uncertainties would take into account the attenuation of the uncertainties due to the prior PDF.