Non-parametric and least squares Langley plot methods

Langley plots are used to calibrate sun radiometers primarily for the measurement of the aerosol component of the atmosphere that attenuates (scatters and absorbs) incoming direct solar radiation. In principle, the calibration of a sun radiometer is a straightforward application of the Bouguer–Lambert–Beer law V = V0e −τ , where a plot of ln(V ) voltage vs. m air mass yields a straight line with intercept ln(V0). This ln(V0) subsequently can be used to solve for τ for any measurement of V and calculation of m. This calibration works well on some high mountain sites, but the application of the Langley plot calibration technique is more complicated at other, more interesting, locales. This paper is concerned with ferreting out calibrations at difficult sites and examining and comparing a number of conventional and non-conventional methods for obtaining successful Langley plots. The 11 techniques discussed indicate that both least squares and various non-parametric techniques produce satisfactory calibrations with no significant differences among them when the time series of ln(V0)’s are smoothed and interpolated with median and mean moving window filters.


Introduction
Langley plots are used to determine the instrumental constant V 0 , i.e., to calibrate, sun radiometers from a series of measurements V i at various air masses m i .According to the Bouguer-Lambert-Beer (BLB) law, the optical depth τ is determined from pairs of points (V i , m i ) that are fit to the linear equation ln (V ) = ln (V 0 ) − τ • m. (1) If τ is constant, the equation defines a straight line; the graph is called a Langley plot.When data are not perfect and contain outliers (τ is not always the same for all measurements when time t and air mass m(t) change), the Langley plot is obtained after removing the outliers.Thus, one can still obtain V 0 from Eq. (1).The derived instrumental constant V 0 , if valid, is used to retrieve the optical depth τ (m) for any measured V and calculated m: τ (m) = −m −1 ln(V /V 0 ).The main purpose of sun radiometry is to retrieve the optical depth of atmospheric constituents, mainly aerosols, but also O 3 , CO 2 , SO 2 , H 2 O, etc. from the direct beam's atmospheric transmittance V /V 0 .
After multiplying the transmittance by the extraterrestrial solar flux, one obtains the flux measured by the instrument.The flux at the site of the measurements can be used, e.g., to validate radiative transfer models (Mlawer et al., 2000).The multi-filter rotating shadowband radiometer (MFRSR) (Harrison et al., 1994) and the rotating shadowband spectroradiometer (RSS) (Harrison et al., 1999) measure the diffuse components of the flux, as well as the direct components.These instruments were calibrated via the Langley method and with standard lamps (Kiedron et al., 1999).Schmid and Wehrli (1995) concluded that when retrieving optical depth the calibration with Langley plots (Langleys) is superior to the calibration based on standard light sources; however, the quality of calibration by Langleys depends on the site at which the instrument is located.
The Langley plot method of calibration consists of locating a subset of raw data points (ln(V i ), m i ) to which a straight line can be fit.The intercept of the straight line estimates the calibration constant ln(V 0 ).This is a consequence of the BLB law.Often, already processed data of retrieved optical depth Published by Copernicus Publications on behalf of the European Geosciences Union.P. W. Kiedron and J. J. Michalsky: Non-parametric and least squares Langley plot methods τ i require verification and correction if the sun radiometer was poorly calibrated or its calibration has drifted.The BLB law implies that the slope of the straight line fitted to points τ i , 1/m i estimates the calibration constant correction factor (Cachorro et al., 2004(Cachorro et al., , 2008)).Finding this slope is a mathematically equivalent approach to finding the intercept in the Langley plot method.In both approaches the calibration success hinges on correctly locating the subset when the optical depth of atmosphere is constant.
It is often overlooked, that the presence of the straight line in the data set ln(V i ), m i or τ i , 1/m i does not imply that the actual τ is constant.The existence of a straight-line fit is the necessary condition but not a sufficient one for the constancy of the optical depth.Shaw (1976) may have been the first to point this out.He observed that in cases when the optical depth of aerosols is a parabolic function of time the pairs ln(V i ), m i create a perfect straight line, but its intercept is not the actual ln(V 0 ).Then the intercept is biased, and thus it cannot be used as a calibration constant.Several authors - Tanaka et al. (1986), Nieke et al. (1999), Harrison et al. (2003), and Campanelli et al. (2004) -mention this problem more or less explicitly.More recently, Marenco (2007) devoted his paper to this phenomenon.Equation (1) implies that when the actual τ A contains a varying component that is inversely proportional to the air mass the data align along the straight line with a slope τ , but the intercept is now ln(V 0 ) − ε.The hyperbolic dependence 1/m produces a straight line.
The time series of ln(V 0 )'s over several days are used to weed out cases when ε is not zero and estimate the true calibration.An individual Langley plot cannot identify the value of ε.The data (ln(V i ), m i ) do not contain information on the presence of a nonzero ε.
The process of removing outliers from the Langley plot may actually facilitate selecting a straight line from the data that will contain a spurious value of ε.Different methods of removing outliers may cause an inadvertent selection of a different value of ε, some larger and some smaller (positive or negative).
We will call a Langley plot with a nonzero ε an anomalous Langley plot.The anomalous Langley plot cannot be identified because data (ln(V i ), m i ) do not contain information on the presence of a nonzero ε.Statistical analysis of the time series of intercepts ln(V 0 )−ε leads to a better estimate of the calibration constant ln(V 0 ) = < ln(V 0 )−ε> = ln(V 0 )−<ε>, where < > denotes an average used in time-series analysis that usually is a combination of mean and median moving averaging windows.When the statistics of ε is unbiased, <ε> tends to zero as the number of samples in the time series increases.
Mountaintops like Mauna Loa in Hawaii or Izaña on Tenerife provide environments where the constancy of the optical depth (ε = 0) is frequent.With a small standard deviation SD(ε) at mountaintops a smaller number of Langleys is necessary to achieve the desired precision of calibration.A long history of measurements at sites like Mauna Loa also warrants the belief that the statistics of ε there is considered to be unbiased.Thus, within the range of validity of this belief an accurate calibration is possible.
In most places where sun photometers are deployed, periods of stable atmospheres are much less common, and they are frequently interrupted by cloud passages, changes in atmospheric conditions like varying humidity that promulgate aerosol size changes and by aerosol plume incursions.Large numbers of outliers call for special Langley plot analyses going beyond standard straight-line fitting procedures.The statistics of ε is likely to be biased with a large standard deviation SD(ε).This necessitates a larger number of points in the time series to achieve a desired precision while the frequency of Langley events at sites like these is low.There is no guarantee that the statistics of ε is unbiased.It should be emphasized that the difference between easy sites like mountaintops and difficult sites like Billings, Oklahoma, is quantitative, not qualitative.The same statistical analysis of time series must be applied in both cases; however, precision and accuracy of results will differ.
When there are no other available independent measurements, time-series analysis is the only option for in situ calibration of sun photometers.Photometers that also measure aureole radiance simultaneously with direct solar flux can be calibrated when optical depth is not constant (Tanaka et al., 1986;Nieke et al., 1999;Zieger et al., 2007).These methods can identify anomalous Langley events and estimate the value of ε.Then, in principle, a single Langley is sufficient.
The main objective of the paper is to analyze the efficacy of non-parametric and least squares methods of straight-line fitting to identify Langley plots useful for calibration.We use time-series analysis only to determine the impact of the methods on the estimated calibration constant.While we are not concerned with the identification of anomalous Langley plots, their presence is manifested in outliers of the time series.We do not deal with minor issues related to the differences in air mass among various air constituents: aerosols, ozone, and molecular scattering, in particular, and higherorder effects like atmospheric refraction's dependence on wavelength impact air mass.In other words, we presume that the BLB as given by Eq. ( 1) is valid.
Non-parametric statistics is the branch of statistics where no a priori assumptions on the statistical distribution of variables are made (Kendall, 1938).Non-parametric does not mean that there are no parameters but that the involved parameters have no assigned statistical properties a priori.For instance, a histogram is used as a substitute for the probability distribution function rather than, for instance, a fitted Gaussian to the histogram.For Langley plot processing the non-parametric methods are particularly appropriate.We do not know the statistical distribution of outliers in a Langley plot even if we understand their physical origins.We do not know the statistics of ln(V 0 )'s in a time series of them, but we still want to get the best estimate of the calibration constant of a sun photometer.
For simplicity of notation in the rest of the paper Eq. ( 1) is replaced with a linear equation: The organization of the paper is as follows.In Sect. 2 we define a Langley plot.In Sects. 3 to 7, 11 methods of finding a Langley plot are described: in Sect.3, two least squares methods; in Sect.4, the so-called objective algorithm method; in Sect.5, four non-parametric regressions methods following Theil (1950) and Siegel (1982); in Sect.6, a non-parametric method of identifying outliers and a modified Siegel (1982) method with sequential removal of outliers; and in Sect.7, a method of analyzing histograms of slopes and intercepts.In Sect.8 we describe the set of data used in the comparisons for all methods.Section 9 presents results of these analyses and comparisons.The final section summarizes the paper.

Our definition of a Langley plot
For any set of points P = {(x i , y i ) : i = 0, . .., n − 1} and for any nonnegative δ find a subset L ⊂ P for which a line y = α + βx can be defined such that one of the following metrics that measures the magnitude of residuals is smaller than δ( <δ) on the set L, where the r j = y j − α − βx j are residuals and k is the number of points in the subset L. We refer to the points of the subset L as δ-collinear.The subset L is not unique, or it may not exist when δ is too small, ignoring subsets consisting of two points only.Therefore, we add a requirement that L should be the largest subset with this property of residuals.In other words, a Langley plot is the most numerous δ-collinear subset of set P .The size of the set L (i.e., the number of points that "actually" define the line) is important in judging the viability of the resulting Langley plot (i.e., the subset L).The quotation marks around "actually" are justified for non-parametric methods, because they do not identify outliers explicitly.This problem is related to the pattern recognition problem.The human eye and mind are able to solve the problem in a qualitative way very quickly by identifying points that are approximately collinear.The human eye and mind can perform this task regardless of the plot orientation.The result is rotationally invariant: neither of the axes x or y is treated preferentially.Furthermore, the human eye and mind can almost instantaneously identify a data set P that has no potential of containing any subset L of a significant size and rejects this case as not providing a viable Langley plot.Mathematically this problem reduces to the straight-line fit and to a method of identifying and removing outliers.Any one of the criteria (Eq. 3) can be used to define the quality of the fit.
Some researchers (Augustine et al., 2003) seemingly avoid the issue of removing outliers altogether by selecting clearsky days based on the method of Long and Ackerman (2000) for which measurements from collocated broadband shaded and unshaded pyranometers are required.This approach, while effective, misses many Langley plots from partially clear days, so it does not fit the scope of this paper.Furthermore, many sites with sun radiometers do not have collocated shaded and unshaded pyranometers.

LSF with sequential removal of outliers (SRO)
This is the most straightforward and, probably, the most commonly used method in existence.The least-square fit (LSF) is applied to the set of points P , and the largest residual (negative or positive) is removed.Then the root mean square (rms) of residuals is calculated.The process is repeated until rms ≤ rms max and the number of remaining points k ≥ k min , with rms max and k min chosen based on experience.Usually, most of the outliers are negative (e.g., cloud passages), but there are less frequent cases when the atmosphere has periods of stability at larger τ that may be temporarily interrupted by a cleaner air mass.For this case the outliers at smaller τ are positive.For this reason the method must allow removal of the positive outliers.On our data set we note that we obtained good results (meaning that both the number of false and missed identifications of Langley plots were small) when about every fifth outlier that was removed was a positive one.However, we do not claim the value of five is a general rule.
Usually measurements with sun radiometers are performed at equal time steps.This means that the values of x, the air masses, are not evenly distributed.For 1-minute intervals x at x = 2 might be 10 times smaller than x at x = 6 at midlatitudes, where the air mass x is the Rayleigh air mass (Kasten, 1965).Some researchers recognized the bias introduced by the uneven distribution of x on α due to the larger number of points at low air masses.For example, Forgan (2000) performed Langley plots on (y/x, 1/x) for his sun photometric studies.In the Dobson ozone spectrophotometer community Langley plots for the ozone extraterrestrial constant (ETC) are performed in coordinates (y/x, 1/x) to give a smaller weight to points that are more sparse at large air masses (Dobson and Normand, 1958).Note, however, that for the Brewer UV ozone spectrophotometers the standard equation y = α + βx is used (Redondas, 2005;Ito et al., 2014).
The change of variable from x to 1/x leads to the linear equation y/x = β + α(1/x).This approach is equivalent to P. W. Kiedron and J. J. Michalsky: Non-parametric and least squares Langley plot methods applying weights w = 1/x 2 to the LSF of the original equation y = α + βx.Herman et al. (1981) considered applying other weighting methods.
A sequential removal of outliers can be applied and the method may have the same terminating criterion in terms of rms < rms max ; however, the residuals must be calculated for the equation y = α + βx.
We label these two methods LSF SRO − x and LSF SRO − 1/x, where the subscript SRO stands for "sequential removal of outliers".

The objective algorithm of Harrison and Michalsky
We describe some aspects of the objective algorithm (OA) because (a) its development is an excellent example of how a mathematical method was stimulated by the human eye-andmind approach, (b) it is based on physical phenomena that are responsible for the curve shape and the outliers, and (c) it is basically a non-parametric method despite the fact that LSF is used for the final filtering.
When Harrison and Michalsky (1994) developed the OA they tested it by comparing a set of cases from 384 days using 500 nm channel data where α and β were obtained by the eye-and-mind method of Michalsky, who disqualified the non-viable cases and identified the ones that, after the removal of outliers, produced Langley plots.Then he performed the LSF on the retained points.The OA did not try to produce "an artificial intelligence" emulating Michalsky's approach.Instead it identified several physical phenomena (like cloud passages, overcast skies, curvature in the plot, etc.) that were responsible for outliers and the non-linearity of the Langley plots.This justifies the term "objective" in the method's name.The method applies consecutive filters: each meant to deal with outliers produced by one of the identified physical phenomena that produced them.The last filter is LSF that shaves off outliers larger than 1.5 of the standard deviation of all residuals of points that survived the previous filtering.
The successful Langley plot for the OA is the one for which rms ≤ 0.006 (of retained residuals) and the k/n ≥ 1/3.The value of rms max = 0.006 was chosen to maximize the agreement with the eye-and-mind Michalsky method using 143 cases of successful Langley plots.The value of rms max is valid for the wavelength λ = 500 nm.For other wavelengths the rms max will be different because aerosols' impact on outliers is wavelength dependent.

Non-parametric fits (NPFs)
LSFs use means of {x i }, {y i }, {x 2 i, }, and {y i, x i }.The so-called breakdown point of the mean is 0 %: a single outlier can significantly change the value of a mean.On the other hand, the breakdown point of a median is 50 %.Theil (1950) opened the field of the so-called non-parametric regression fits that are based on medians and, thus, are much more robust.
The set P produces an nxn matrix of all possible slopes {b i,j }, where b i,j = (y i −y j )(x i −x j ).The matrix is symmetric with a diagonal that has indeterminate values.Its upper or lower triangles each have n(n−1)/2 points.They are used to calculate the slope: From the slope β the intercept is obtained also as a median: (5) Theil's (1950) algorithm robustness is 29.3 %, which means that when outliers exceed 29.3 % of all points the performance of the algorithm is not guaranteed; at this level it reaches its breakdown point.The increase of robustness to 50 % was achieved by Siegel (1982) with his method of "repeated medians" that uses two medians in Eq. ( 4) rather than one: one median along the rows of the matrix {b i,j } and then the median of this column of medians where all n(n − 1) values of the matrix {b i,j } are used.The Langley plot is chiefly concerned with obtaining the intercept and, unlike Theil's focus, the slope is secondary.Instead of obtaining the slope first, one can obtain the intercept first.From the matrix {a i,j } of intercepts a i,j = (y j x i − y i x j )/(x i − x j ), one gets the intercept α with Eqs.(4) or ( 6) and then gets the slope from The weighted median methods (Jaeckel, 1972) can also be applied.The uncertainty of slopes b i,j stems from the measurement errors of y i values.The uncertainty is inversely proportional to |x i − x j | if uncertainties for y i are the same.When |x i − x j | is small, the measurement errors have inversely proportional 1/|x i − x j | larger impact on the slope.For the intercepts the weights for a i,j are proportional to x i − x j /(x 2 i + x 2 j ) 1/2 .The exact formulas for using weighted medians for Theil (1950) methods can be found in Birknes and Dodge (1993).The weighted medians, one would expect, should offer an advantage when x i 's are not uniformly distributed, which is exactly the case for air masses.However, our simulations with weighted medians did not confirm this expectation.In fact, in our experience the weighted median methods introduced unacceptable biases in α and β.
The LSF methods have several advantages over NPF methods.They perform better when data have no outliers, and they have uniformly distributed noise among y values, particularly when the noise is Gaussian.The solutions for α and β in LSF methods are closed-form formulas derived from smooth analytic functions.The NPF methods are inherently discrete.They depend on medians.The "granulation effect" in a small data sample may lead to errors because of discontinuities.The residuals from the LSF methods are unbiased (sum of residuals equals zero), while this is not guaranteed for the methods based on the discreet processes.Theil (1950) and Siegel (1982) NPF methods do not identify outliers explicitly.α and β are generated for any set P , but the outliers must be identified to get the value of the metrics and reach the decision of whether these particular α and β define a Langley plot or not.The last problem we solved by using the following method of outlier identification.
For a given α and β we sort all points according to the ascending value of their residuals (r j ≤ r j +1 ).Then we calculate the root mean square of residuals (rms J ) of the first j = 0, . .., J − 1 points.We find J for which rms J ≤rms max and rms J +1 > rms max .All points with indices j ≥ J are considered outliers.Points with indices j < J are retained.They define the Langley plot, which is considered to be successful when J ≥ n/3.Residuals of retained points can be de-biased by performing a LSF on them.This reduces the value of rms J and slightly changes α and β.
We chose to use rms (the first criterion in Eq. 3) because the data set for OA, which is used in comparisons, is based on rms metrics.The method of identifying outliers in the set P , when α and β are given, we label OSM (outlier sorting method).In Sect.9 we also apply the OSM to the results obtained from the OA method.
We described four NPF methods of finding a Langley plot.We label them as T OSM -β, T OSM -α, S OSM -β, and S OSM -α, where T and S stand for Theil and Siegel, respectively, and α and β stand for the "intercept-first" and "slope-first" methods and OSM for the outlier sorting method with a final residual de-biasing LSF.The de-biasing on average increases the intercept α.At most (S OSM -β method) α increases by 0.0028 (0.28 % in terms of V 0 ).

Identifying outliers from the dispersions of slopes
Neither the Theil (1950) nor Siegel (1982) methods identify outliers explicitly.They produce the slope and intercept directly for any set P .In this section we describe a nonparametric method that identifies outliers without calculating the values of residuals.
For each row i of matrix {b i,j } we calculate d i , which is the measure of dispersion among the points of the ith row.d i can be the standard deviation of the row or its median absolute deviation (MAD).We used the latter.The "largest" outlier is the one for which d i is largest.Then we remove row i and column i from the matrix {b i,j } and calculate new values d i and find the one that is the largest, and so on.The largest dispersion D m , where m is the index of the iteration process, forms a descending sequence with a decreasing steepness.
Large drops in the sequence indicate a removal of a significantly "large" outlier.The largeness or smallness of outliers should be understood as their values of dispersion, though it may correlate very well with the value of residuals from the straight line y = α + βx.
One can analyze the sequence of {D m }.Once it flattens, this indicates that points that remain approximately a straight line and the process of outlier removal can be stopped, but if {D m } remains strongly decreasing it implies that there is no "collinear" subset in the data.We did not explore the potential of finding a criterion for stopping the iteration process from features and behavior of the {D m } sequence.Instead we used the remaining points to calculate rms and stopped the process when rms became smaller than rms max .
We applied the same method to the matrix of intercepts {a i,j } and obtained similar results: the sequences of removal of "large" outliers for both {a i,j } and {b i,j } matrices were the same but not identical when only "small" ones were left.
This approach of outlier identification and removal led us to a modified Siegel (1982) method.At each stage when a row and a column are removed we calculate new α m and β m with the Siegel (1982) method.Initially we were surprised that after a removal of an outlier the new values of α m+1 and β m+1 were not always changing significantly until we realized that this is a consequence of the robustness of the Siegel method.The method is stopped when rms m ≤ rms max , and the result is retained if (n − m)/n>1/3.This method of finding a Langley plot we label S SRO -β or S SRO -α, where SRO stands for a sequential removal of outliers.Keep in mind that SRO from this section is strictly a non-parametric method unlike SRO in Sect. 3 on LSF methods.

Histograms of slopes and intercepts
In this method we analyze a histogram of slopes to identify the subset of nearly collinear points.The histogram of slopes is constructed from elements of the matrix {b i, j }.We locate the cell [b, b + b] at which the histogram's counts are maximum, where b denotes the width of histogram cells.Next we identify all pairs of points [(y i , x i ), (y j , x j )] that produce slope b i,j from within the interval [b, b + b].We are interested in points y i , x i that create many such pairs.Let c i denote the number of such pairs that the ith point creates.The points for which c i = 1 are rejected, and the median of remaining c i is calculated.Then the points with c i less than the median are rejected.LSF is applied to the remaining y i , x i points and the rms is calculated.This method defines the Langley plot with a small number of points, but it is very efficient at detecting subsets of collinear points.For this reason we did not require that the number of retained points had to be larger than n/3 for this method.We labeled this method H-β.A similar process can be performed with the histogram of intercepts.The results, however, were not as good as with the histogram of slopes.
In Fig. 1a-d we show three cases in which the H-β method fails and one case with a large number of outliers for which H-β works correctly.The cases are shown to illustrate the usefulness of a histogram analysis in prescreening cases and possibly designing a more sophisticated method that could produce not one but several Langley plots from one set of points P .
In Fig. 1a there are two regions: 2 ≤ x ≤ 2.12, which contains 46 points, and 3.51 ≤ x ≤ 4.16, which contains 19 points.These two regions produced two mono-modal, narrow histograms, implying that there are many collinear points from within two regions.
In Fig. 1b there are two regions with collinear points.However, each region has a different slope, and it extrapolates to a different α.Both histograms are bimodal.Two distinct Langley lines could be produced in this case.Only one, if either, can be right.The question that one of them or both are anomalous Langley plots can be posed.
Figure 1c shows a very interesting case.The histogram of slopes is mono-modal, but the histogram of intercepts is bimodal.The height of the second mode is less than half of the dominant mode.Two regions 2.25 ≤ x ≤ 3 and 3 ≤ x ≤ 6 have similar slopes as they produce the mono-modal histogram of slopes.But at x = 3 there is a step change.It is not possible that it was produced by a change in the optical depth if one excludes the change in ε responsible for the anomalous Langley.It is possible that at x = 3 something affected the responsivity of the instrument.In this case we will get two Langley plots that are almost parallel with different α's.
Figure 1d depicts a case without major outliers.The points can be fairly well approximated with a third-degree polynomial (a thin line is depicted), which means that τ is the second-degree polynomial of air mass.Both histograms are bimodal and rather broad.One may pose the question of whether histograms could be used to detect nonlinearity.The method proposed by Kuester et al. (2003) has the potential for detecting the nonlinearities; however, it seems that the authors did not explore this possibility.

The data set
The comparison and analysis of Langley algorithms was performed on a data set produced by the rotating shadowband spectroradiometer (RSS) (Harrison et al., 1999) deployed at the Department of Energy Southern Great Plains We removed many overcast days and several corrupted files that, together with the instrument down times, reduced the data set to 1023 morning or afternoon sets P = {(y i , x i ) : 2 ≤ x i ≤ 6}.Data from one pixel, out of 1024, at approximately λ = 500 nm are used in this analysis.
The RSS was lamp calibrated every 2-3 weeks with calibrators that had lamp calibrations traceable to a NIST standard.The values of V were normalized by the responsivity obtained from each lamp calibration and interpolated between the calibration days.This reduced trends in V 0 due to the instrumental instability caused by optical elements aging, diffuser degradation, and CCD response changes.Nevertheless, the RSS displayed quasi-periodic instabilities due to what we later discovered was an outgassing problem that led to a deposition of a thin film on the cooled windowless CCD.This resulted in a wavelength-and time-dependent etalon effect that affected the CCD's response.Lamp calibrations mitigated the effect but did not remove it from the data completely.

Comparison of methods
In the previous sections we described 11 methods to identify a Langley plot.In this section we compare them at four different values of rms ma x = 0.010, 0.008, 0.006, and 0.004 using one data set.First we look at some statistical parameters concerning α between each two methods, and then we look at calibration constant time series derived from each method.
In Table 1 we collected information on the number of Langley plots for each method (in the diagonal of the table) and the number of common Langley plots between the two methods (above the diagonal).Then we included some statistical parameters on α between each of the two methods; there are 55 combinations.Above the diagonal is the standard deviation of α, and below the diagonal mean and median of α.The order of subtraction in α is as follows: α for the method from the row minus α for the method from the column.
The LSF methods produce the largest number of Langley plots (601 and 598) while the OA the lowest (284).The OA's α's are larger than any other method by 0.0026-0.0065,which translates to 0.26-0.65 % in V 0 .For most cases medians of α are 5 to 10 times smaller than means of α.This is because the main contribution to differences α among methods comes from the tails of α distributions.In other words, outliers are responsible for the main differences among the methods; however, some biases exist among them.A large median α indicates that the bias between the methods given by a mean is real and does not apply only to the outliers.In Fig. 2 we demonstrate the effect of the number of outliers on the differences among the methods.The average of the standard deviations of α for all methods is plotted against the percentage of points retained by Langley methods.When Langley plots have no more than 10 % outliers, the standard deviations between the methods are an order of magnitude smaller than when the number of outliers is up to 67 %.This implies that the main differences between methods are due to different handling of outliers by each method, and the outcome is more method dependent when the Langley plot consists of a smaller number of points.Also, we plotted the number of Langley plots vs. the number of points remaining in the Langley plot for the LSF SRO − x method to demonstrate how strongly the number of available Langleys diminishes with the number of outliers for the site in Oklahoma.Also, we plotted the number of Langley plots vs. the number of points remaining in the Langley plot for the LSF SRO − x method.
The Theil and Siegel methods (T OSM -β, T OSM -α, S OSM -β, S OSM -α) produce very similar results with some of the lowest means and medians of α.In some cases medians are zero.We could not discern a difference between the "intercept-first methods" (T OSM -α, S OSM -α) and the "slope-first methods" (T OSM -β, S OSM -β).
LSF SRO −x yields larger α than LSF SRO −1/x (by 0.0055).For rms max = 0.010, 0.008, and 0.004, it is 0.0083, 0.0065, and 0.0039, respectively.The standard deviation of 0.0198 is not exceptionally large or small in comparison with other methods.
The Siegel methods with sequential removal of outliers (S SRO -β, S SRO -α) yield significantly different numbers of Langley plots (364 vs. 475).But on the common set the mean and median of α are very small.
Table 1.The comparison among 11 Langley methods for rms ≤ 0.006.Number of successful Langley plots is on diagonal.Above the diagonal the number of common Langley plots for two methods and a standard deviation of differences α.Below the diagonal the mean and median of α.The order of subtraction in α is as follows: α for the method from row minus α for the method from the column.The outlier sorting method when applied to OA increases the number of Langley plots by 75 %.On the common set of data the OA OSM produces smaller α's (by 0.0028).This is the opposite effect of OSM compared to NPF methods (see Sect. 5).The extra Langleys produced by the OA OSM method do not necessarily indicate an improvement.Many of them are large outliers in the time series.We conclude that the OA, if it errs, it errs on being conservative: it has a fairly large missed detection rate (rejecting data sets with viable Langley plots), and at the same time the ones that are detected sometimes could be improved by a removal of few extra outliers.
When evaluating individual plots, and we looked at almost all 11 × 1023 of them, we found for each method cases when it went astray.There were cases of missed detection and false detection when judged by eye.However, we cannot quantify which of the algorithms has the most favorable missed and false detection rates.Out of all algorithms used, only OA deals explicitly with curvature.This perhaps might be a chief reason why it produces significantly fewer Langleys, which leads to a smaller number of large outliers in the time series.
The comparison of calibration constants that can be derived from α's obtained by each method gives us additional insight about each method as well as a strategy one should use when generating the calibration constants.We compare the behavior of time series of derived calibration constants α cc (d), where d indicates each day from 5 October 2003 to 30 March 2006.The α cc (d)'s are obtained from the time se-ries of α(d j ) independently for each Langley method, where d j denotes days at which Langleys were obtained.The α cc (d) might be considered "the best" estimate of the calibration constant for a given day, "the best" in the sense of the method that we use to remove outliers, interpolate and smooth the series α(d j ).
The method consists of a moving median window of width d med days that removes outliers and interpolates, which is followed by a moving boxcar filter.For each day d the median is calculated from d med number of α(d j ) values: d med /2 values for d j ≤ d and d med /2 values for d<d j .Then the 1055long series α cc (d) is smoothed with a boxcar filter of d smth days.By trial and error we decided on d med = 30 days and d smth = 25 days.
Prior to applying the method described above, the values of α are corrected for Earth-Sun distance a by a substitution α ← −α + 2ln(a), where a is in astronomical units.The Earth-Sun distance is calculated with the ephemeris program published by Michalsky (1988).
One of the methods of time-series smoothing of V 0 's to obtain the absolute calibration constants of a sun photometer (MFRSR) was validated against calibrations at Mauna Loa by Michalsky and LeBaron (2013).The discrepancy between the time-series smoothing and Mauna Loa calibration constants in terms of V 0 's was always smaller than 0.6 %.This at air mass m = 2 translated to aerosol optical depth uncertainty of less than 0.003 and at m = 6 to less than 0.001.This MFRSR was located in Boulder, Colorado, which in terms of number of sunny/clear sky days is somewhat superior to the SGP site.The OA method was used to identify the Langley plots and obtain individual V 0 's.In Fig. 3 we show all 11 α cc (d) curves and individual intercepts α's from OA and H-β methods for rms ≤ 0.006.In the course of 1055 days the RSS's calibration constants vary within ±3.5 % (in terms of V 0) band.All methods follow these changes; however, there are differences among them.Statistically, the differences are ±1.4 for 95 % of days from the calibration constant curve that is the average of all 11 curves.For other values of rms max 0.010, 0.008, and 0.004, the differences are ±1.6,±1.45, and ±1.4 %, respectively.The effect of the maximum rms on the differences between the α cc (d) curves is not very dramatic.This is because the parameter d med = 30 days of the median filter is relatively large.
In Fig. 4 we show nine calibration constant curves and intercepts α's for LSF SRO − x and S SRO -α methods for rms ≤ 0.010.In this case we used α's from Langley plots that had no more than 10 % outliers.The OA and H-β results did not pass the filter of 10 % outliers only.The differences between the calibration constant curves are ±0.8 %.For 40, 30, and 20 % outliers, the calibration constant curves are within ±1.2, ±0.9, and ±0.87 % bands, respectively.So, the effect of number of outliers removed to obtain a Langley plot has a larger effect on the spread among the methods than the effect of rms max .
We note that the OA and H-β calibration constant curves from Fig. 3 are marginally within the band defined by the curves in Fig. 4.
The majority of points in Fig. 4 are outliers, and they are defined by Langley plots with 90 % or more points.By the criterion rms ≤ 0.01 the points are collinear.Nevertheless, they are off and some by more than ±5 % (in terms of V 0 ).In our opinion, the majority of the outliers are cases of anomalous Langley plots.The topic of anomalous Langley plots will be pursued in another paper.

Conclusions and summary
Eleven Langley plot methods were compared.Two of them were the least square methods and nine were nonparametric methods which included the objective algorithm (OA) method by Harrison and Michalsky (1994).

P. W. Kiedron and J. J. Michalsky: Non-parametric and least squares Langley plot methods
We developed two methods to terminate the nonparametric methods in order to determine the existence of the Langley plot: the outlier sorting method (OSM) that was applied to two Theil (1950) and two Siegel (1982) methods, and a new non-parametric method of sequential removal of outliers (SRO) was applied to two Siegel methods resulting in two new iterative Siegel methods.
We found that analysis of histograms of slopes and intercepts can be an excellent tool to prescreen a data set for Langley-plot viability.The histogram of slopes was used to generate Langley plots that produce lines defined by a small number of points.The histogram method offers a possibility to extract Langley plots when outliers dominate and to find all subsets of collinear points.
The OA method turned out to be robust though conservative.It identifies the lowest number of Langley plots.It produces intercepts slightly larger than all other methods.
The Siegel (1982) and Theil (1950) methods with OSM produce very similar results.The two least square methods yield the largest number of Langley plots, with expected bias between them.
The largest differences among methods are on Langley plot cases that turn out to be outliers in terms of the calibration constant curve.Predominantly these are the cases that produce Langley plots but with a small number of points.In cases that are close to the calibration constant curve the differences are small, but there are systematic biases.
We have no way of determining which of the methods produces results closest to the truth.In fact, the answer may depend on the data set.When the number of outliers in a Langley plot is small, all methods tend to produce similar results.
The metrics used to define the Langley plot was rms of residuals.The effect of the value of rms, whether it was 0.10 or 0.06, had no great impact on the calibration constant curves: all methods produced calibration constant curves within a band between ±1.4 and ±1.6 % for 95 % of days.It is the number of outliers in the data set that has a greater impact.The calibration curves generated using a smaller number of Langley plots with each Langley defined by a larger number of points produce calibration constant curves that are less dependent on the method.For instance when Langley plots retain 80 % of the points all calibration constant curves are within ±0.9 % band for 95 % of days.
The outliers from the calibration constant curves are predominantly caused by anomalous Langley plots when the optical depth has a hyperbolic component as a function of air mass.This effect cannot be detected from the data set, and no Langley plot method can determine if this hyperbolic change with air mass is occurring.This effect at difficult sites like the SGP ARM site in Oklahoma sets the ultimate limit of accuracy of in situ calibrated sun photometers.

Figure 1 .
Figure 1.Four cases that illustrate how slope and intercept histograms can aid in analysis of points for the extraction of the Langley plot.In (b) and (d) both the slope and intercept histograms are bimodal.In case (c) only the histogram of intercepts is bimodal.Case (a) has many outliers, but both histograms are mono-modal indicating the existence of a single Langley plot.Case (d) has no large outliers, but y vs. x is nonlinear (third-degree polynomial).

Figure 2 .
Figure 2. The average over all methods of standard deviations of α as a function of number of points in a Langley plot.

Figure 3 .
Figure 3. Calibration constant curves for all 11 methods (rms < 0.006) and individual intercepts α for OA and H-β methods.

Figure 4 .
Figure 4. Calibration constants curves for nine methods (rms < 0.01) from Langley plots with no more than 10 % outliers and individual intercepts α's for LSF SRO − x and S SRO -α methods.(The OA and H-β results did not pass the filter of 10 % outliers only.) Its results are most similar to results produced by Siegel methods and LSF SRO − x.