Discussion and Conclusions
Results with 5-sigma outlier rejection
Results with 4-sigma outlier rejection
Details
The iteration
The data set I am working on here is the set of "Sorted" files from Tom's analysis that were posted from July to November 2002. I have analysed them using the correction method discussed for DS24 (Orion data set.) This involves the subtractions of spatial polynomials. I used a 4th order fit to the whole data set and 1st order corrections for the individual images. In the previous work with the 7 x 7 grid of images, I used 2nd order spatial fits to individual images but with the very much poorer spatial coverage of this data set, this was not attempted. There really isn't enough data to fit even 1st order gradients in declination for most parts of the data set. I did it anyway ...
The method I used consists of iteratively purging the data set of high-scatter images while fitting the spatial polynomials, using progressively tighter tolerances as the iteration converges. This process is completely automatic, controlled by adjustable parameters which set the point at which outliers are rejected and are basically multiples of the standard deviation. Yes - this does mean that I am assuming that the underlying measurement errors for the best images are close to Gaussian. That's a circular definition of "best images"!
Verdict first; evidence later. The pretty pictures follow this section and I've relegated the details to the end where you may safely ignore them.
Images vary a lot in quality. The algorithms presented here allow the use of automated methods to reject the worst images. Apart from setting the analysis parameters, no manual input is required. I didn't look at a single image.
There seems to be general agreement that the problem is the sky and that Tom needs to move away from Chicago! The spatial correction used here beats the residual standard deviation down to 0.0081 magnitudes in V-band and 0.0060 magnitudes in I-band. This is by no means as good as it sounds! Individual images with standard deviations considerably higher than this are currently left in the data. Remember that the error distribution is far from Gaussian and 4-sigma errors are not rare. Peak errors of plus or minus 4-sigma then give you an error band around 0.06 magnitudes! Ouch!! And it happens: what's more, as the sky changes it can and does mimic plausible looking star variability. User beware. I can tighten up the image rejection criterion some more but more data will be thrown away. We had better face up to it; nobody is going to detect planet transits in this data.
Can one do better? Maybe a little. My work on the 7 x 7 grid showed a small but significant improvement on going to quadratic correction of each image. This is not possible with this data set, which barely defines linear corrections. A data set with consistent overlaps between images in both RA and declination would be better. One might wring a little more by carrying the overall fit to higher order. But basically it is not going to get noticeably better than what we are doing now, at least, not without moving to a better site.
How good is that? Well, my code eventually identified and eliminated short period variables with 0.3 magnitudes range with essentially 100% certainty; these can obviously be picked up again by the appropriate processing using the final fits. At the 0.1 magnitude level, things are still a bit uncertain. That 0.06 magnitude error band is real! I think the code is detecting most of the 0.1 magnitude short period variables but it picks up a lot of bogus ones too so tedious manual processing is needed - or a smarter algorithm to pick them out better. At 0.06 magnitudes, the current code picks at least 10 bogus variables for every real one. I have still not found a long period variable.
This shows the degree of improvement in V-band obtained by sytematically throwing out 5-sigma outliers and applying the spatial corrections. The first iteration applies a null correction and has very broad limits for rejecting residuals so this is essentially raw data.
There are initially 4,436 stars seen 40 or more times in a total of 342,341 accepted observations on 4,992 images. By the end of the iteration 3,197 stars are left with 182,921 accepted observations on only 2,604 images.
Spatial corrections are somewhat larger in I-band so the degree of improvement is somewhat greater.
The biggest improvement of all is in the Welch-Stetson variability parameter. I still have a small "cooking factor" in there that reduces the computed values for the brightest sources. I'm not sure what is reducing values for the faintest sources. If my sums are right, the standard deviations supplied with the data are actually overestimates by some 10 or 20%, which might have something to do with it.
Correction is obviously still not perfect! I have examined all 82 of the sources giving Welch-Stetson statistic greater than 2.0 and I reckon that only 6 of them are really variable ... maybe fewer. Remember, however, that the code has thrown out all the large amplitude variables. The ones found have amplitudes 0.3 magnitudes or less.
This time, at the end of the iteration 2,882 stars are left rather than 3,197. 161,868 observations are accepted instead of 182,921 and only 2,270 images remain. In return for throwing away all this data, the cigar shaped spread of standard deviations is appreciably narrower.
The scatter in the standard deviation is still too large. For Gaussian statistics, my trusty F-distribution tables tell me (with a bit of interpolation) that 99.9% of standard deviation estimates for 40-sample sets should fall in a range of 2.1:1. Except at the low brightness end where the noise is dominated by the measurement noise, it's a lot more than that - more like 6:1.
Same thing in I-band. Nothing dramatic but definitely less spread than before. Worse scatter than V-band; maybe 8:1
A considerable reduction in scatter here. Where there were 82 stars with Welch-Stetson parameter above 2.0, there are now 39. Three of the six claimed variables have also disappeared but this is expected because the algorithm is deliberately trying to eliminate variables! Somewhat disconcertingly, 10 of the 39 variables are new. None of them looks like a real variable to me and three of them are close together and show almost identical I-band structure that is obviously the result of poor correction of one night's observations. Clearly, I need to tighten up the criteria and throw out more images!
Here is the final Welch-Stetson plot on a bigger scale. Totally unretouched:
that's all the points. Notice that the points scatter around a mean of 0.32
rather than zero. This is a measure of how well the consistent
image to image sky variations that affect both colours are being taken out.
Not perfect but a good deal better than the raw data where the vastly greater
scatter has a mean of 3.9.
In order to compute and apply spatial corrections, it is necessary to know where the centre of each image is. I derived this from the data set using the method that Michael Richmond posted some time ago, that is, finding the maximum and minimum RA and declination appearing on each image and averaging to find the centre. This appears to be entirely satisfactory. The program that does this provides a list of image ID's, RA and Dec. I used a spreadsheet to format this into an image list, adding image flags, initially all set to accept the image and initial null image corrections. It turned out that no manual processing was required so the image list could have been produced entirely automatically ... once I got the bugs out of the program.
I used a preselected finding list consisting of Tycho2 stars brighter than magnitude 10.5 in V-band and brighter than 10.0 in I-band. This list contained 12,927 stars in a band that included the total range of declination of all the images. Less than half of these are actually covered by this data set. The 5 months of data were preselected into a working file, retaining only measurements of stars in the finding list and keeping only sets with more than a minimum number of observations.
All the fitting used weighting in accordance with the estimated standard deviations included in the data set. This has the advantage of giving the "best" least squares estimates for Gaussian error distributions. The error distributions are clearly not Gaussian so I don't place much importance on this. The important advantage is that it makes the whole process insensitive to the data selection used. I chose to use only bright Tycho2 stars but the precise magnitude cutoffs do not make much difference. The large number of faint stars would contribute something towards the fit even though they are weighted down - I left them out to save on procesing time. I initially used unweighted fitting, which made the choice of cutoff magnitudes very critical, but, once I bit the bullet I was surprised how few changes were needed to the code to convert to weighted fits.
Previous work had run into a noise floor at around 0.008 or 0.01 magnitudes for the brightest stars. This is several times the estimated standard deviations so that bright stars appeared to be grossly overweighted in the fits. One consequence of this was that the Welch-Stetson variability parameter tended to come out very large for the brightest sources - almost all of them appeared to be "significantly" variable. I therefore added a "cooking factor" by increasing all variances by an amount chosen to match the observed noise floor. This code is still in there. As the fits improve, I can adjust the factor downwards. It is still a little high for the results presented here, as can be seen from the Welch-Stetson plots which show reduced scatter for the brightest sources. Current values for the added standard deviations are 0.003 magnitudes for V-band and 0.006 for I-band. This I-band value is much too high but the fitting method is robust with respect to this!
The data set is obviously heterogeneous, some images having much larger scatter than others. The method adopted was to fit the spatial polynomials iteratively. At each stage, magnitudes spatially corrected using the previous fit were computed, ignoring data from images already flagged as having excessive scatter. All data for one star of the finding list was collected. The mean was computed and the magnitudes converted to residuals from the mean. Outliers were then rejected at a multiple of the standard deviation of the previous fit. This was done iteratively until all residuals were within tolerance. If too few points remained after this, the star was rejected. This obviously loses real large amplitude variables but that's what is needed at this point. They can be put back in at the end. No: my current code doesn't do this properly. Standard deviations and Welch-Stetson statistics were computed at this point, logged to an output file and compared to thresholds to determine if the star was to be used in the spatial fitting.
Star data sets that survived this preliminary processing were used to accumulate the Normal Equations for the least squares fitting of the spatial polynomials. Regular least squares methods were used instead of the more desirable Singular Value Decomposition methods because of storage constraints. When all the stars have been read in, the Normal Equations are solved and various statistics are derived. The residuals from the overall 4th order fits are used to estimate the new reference standard deviations for the next iteration (if required) and the residuals from the individual image fits are compared to expected values. Images showing excessive residuals are flagged so as not to be used again. Once flagged, an image is not retested to see if things have improved with the improvement in the fit whereas stars are retested each time round the iteration. That's just the way the program is organised - it is not easy to retest an image.
Iteration is continued until improvements in both the V-band and I-band overall fits are insignificant. The actual criterion used is based on the change in variance of the residuals. For a random data set, this change, appropriately normalised, is distributed as chi-squared on a number of degrees of freedom equal to the number of parameters fitted (15 for the 4th order fit.) On the other hand, when the iteration is completely converged, the change is zero. In practice, the iteration is deemed to be converged when the change is reduced to some amount related to the expected chi-squared distribution - the factor used is yet another adjustable parameter. For this data set, convergence is very slow - around 20 iterations - at least partly because the gradients, especially in declination, are poorly constrained by the data. A data set with good overlaps should converge more rapidly. The 7 x 7 set converged, even with quadratic image correction, in 5 or 6 iterations.
Initial weighted mean standard deviations were 0.0197 magnitudes in V-band and 0.0211 magnitudes in I-band. These are not the values that would be obtained for raw data since all wildly discrepant measurements were rejected even in this inital pass.
Final standard deviations using 4-sigma rejection were 0.0081 magnitudes in V-band and 0.0060 magnitudes in I-band.