Skewness and kurtosis unbiased by Gaussian uncertainties

# Skewness and kurtosis unbiased by Gaussian uncertainties

Lorenzo Rimoldini
Observatoire astronomique de l’Université de Genève, chemin des Maillettes 51, CH-1290 Versoix, Switzerland
ISDC Data Centre for Astrophysics, Université de Genève, chemin d’Ecogia 16, CH-1290 Versoix, Switzerland
lorenzo@rimoldini.info
Draft version: April 29, 2013
###### Abstract

Noise is an unavoidable part of most measurements which can hinder a correct interpretation of the data. Uncertainties propagate in the data analysis and can lead to biased results even in basic descriptive statistics such as the central moments and cumulants. Expressions of noise-unbiased estimates of central moments and cumulants up to the fourth order are presented under the assumption of independent Gaussian uncertainties, for weighted and unweighted statistics. These results are expected to be relevant for applications of the skewness and kurtosis estimators such as outlier detections, normality tests and in automated classification procedures. The comparison of estimators corrected and not corrected for noise biases is illustrated with simulations as a function of signal-to-noise ratio, employing different sample sizes and weighting schemes.

Skewness and kurtosis unbiased by
Gaussian uncertainties

Lorenzo Rimoldini
Observatoire astronomique de l’Université de Genève, chemin des Maillettes 51, CH-1290 Versoix, Switzerland
ISDC Data Centre for Astrophysics, Université de Genève, chemin d’Ecogia 16, CH-1290 Versoix, Switzerland
lorenzo@rimoldini.info
Draft version: April 29, 2013

## 1 Introduction

Measurements generally provide an approximate description of real phenomena, because data acquisition compounds many processes which contribute, to a different degree, to instrumental errors (e.g., related to sensitivity or systematic biases) and uncertainties of statistical nature. While instrumental effects are addressed before data analysis, statistical uncertainties propagate in subsequent processing and can affect both precision and accuracy of results, especially at low signal-to-noise () ratios. Correcting for biases generated by noise can help the characterization and interpretation of weak signals, and in some cases improve a significant fraction of all data (e.g., the number of astronomical sources increases dramatically near the faint detection threshold, since there are many more sources far away than nearby).

In this paper, noise-unbiased estimates of central moments and cumulants up to the fourth order, which are often employed to characterize the shape of the distribution of data, are derived analytically. Some of the advantages of these estimators include the ease of computation and the ability to encapsulate important features in a few numbers. Skewness and kurtosis measure the degree of asymmetry and peakedness or weight of the tails of the distribution, respectively, and they are useful for the detection of outliers, the assessment of departures from normality of the data (D’Agostino, 1986), the classification of light variations of astronomical sources (Rimoldini, 2013a) and many other applications. Various estimators of skewness and kurtosis are available in the literature (e.g., Moors et al., 1996; Hosking, 1990; Groeneveld & Meeden, 1984; Bowley, 1920), some of which aim at mitigating the sensitivity to outliers of the conventional formulations. On the other hand, robust measures might miss important features of signals, especially when these are characterized by outliers (as in astronomical time series where stellar bursts or eclipses from binary systems represent rare events in the light curve) and weighting might help distinguish true outliers from spurious data (employing additional information such as the accuracy of each measurement), so the traditional forms of weighted central moments and cumulants are employed in this work.

Moments are usually computed on random variables. Herein, their application is extended to data generated from deterministic functions and randomized by the uneven sampling of a finite number of measurements and by their uncertainties, whereas the corresponding ‘population’ statistics are defined in the limit of an infinite regular sampling with no random or systematic errors. This scenario is common in astronomical time series, where measurements are typically non-regular due to observational constraints, they are unavoidably affected by noise, and sometimes also not very numerous: all of these aspects introduce some level of randomness in the characterization of the underlying signal of a star.

While the effects of sampling and sample size on time series are studied in Rimoldini (2013a, b), this work addresses the bias, precision and accuracy of estimators when measurements are affected (mostly) by Gaussian uncertainties. Bias is defined as the difference between expectation and population values and thus expresses a systematic deviation from the true value. Precision is described by the dispersion of measurements, while accuracy is related to the distance of an estimator from the true value and thus combines the bias and precision concepts (e.g., accuracy can be measured by the mean square error, defined by the sum of bias and uncertainty in quadrature).

Noise-unbiased expressions are provided for the variance, skewness and kurtosis (central moments and cumulants), weighted and unweighted, assuming Gaussian uncertainties and independent measurements. The dependence of noise-unbiased estimators on S/N is illustrated with simulations employing different sample sizes and two weighting schemes: the common inverse-squared uncertainties and interpolation-based weights as described in Rimoldini (2013a). The latter demonstrated a significant improvement in the precision of weighted estimators at the high end.

This paper is organized as follows. The notation employed throughout is defined in Sec. 2, followed by the description of the method to estimate Gaussian-noise unbiased moments in Sec. 3. Noise-unbiased estimates of moments and cumulants (biased and unbiased by sample-size) are presented in Sections 4 and 5, in both weighted and unweighted formulations, and the special case of error-weighted estimators is presented in Sec. 6. The noise-unbiased estimators are compared with the uncorrected (noise-biased) counterparts with simulated signals as a function of ratio in Sec. 7, including weighted and unweighted schemes and two different sample sizes. Conclusions are drawn in Sec. 8, followed by detailed derivations of the noise-unbiased estimators in App. A.

## 2 Notation

For a set of measurements , the following quantities are defined.

• Population central moments with mean , where denotes expectation, and cumulants , , (e.g., Stuart & Ord, 1969).

• The sum of the -th power of weights is defined as .

• The mean of a generic set of elements associated with weights is .

• Sample central moments and corresponding cumulants .

• Sample-size unbiased estimates of central moments and cumulants , i.e., and .

• The standardized skewness and kurtosis are defined as , and , with population values and . and satisfy consistency (for ) but are not unbiased in general (e.g., see Heijmans, 1999, for exceptions).

• Noise-unbiased estimates of central moments and cumulants are denoted by an asterisk superscript.

• No systematic errors are considered herein and random errors are simply referred to as errors or uncertainties.

• Statistics weighted by the inverse-squared uncertainties are called ‘error-weighted’ for brevity and interpolation-based weights computed in phase (Rimoldini, 2013a) are named ‘phase weights’.

## 3 Method

The goal is to derive an estimator as a function of observables (measurements with corresponding uncertainties ) which is unbiased by the noise in the data, i.e., such that the expectation equals the estimator in terms of the true (unknown) values aimed at by the measurements.

The noise-unbiased estimator is obtained with the following procedure and assumptions. If independent measurements are associated with independent Gaussian uncertainties , the expected value of the estimator is evaluated from measurements and the joint probability density , for given true values and measurement uncertainties :

 ⟨T(x)⟩=∫RnT(x′)p(x′|{\boldmathξ},{\boldmathϵ})dnx′, (1)

where

 p(x′|{\boldmathξ},{\boldmathϵ})=n∏i=11√2πϵiexp[−(x′i−ξi)22ϵ2i]. (2)

As shown in App. A, the expectation of the estimators considered herein can be decomposed as

 ⟨T(x)⟩=T({\boldmathξ})+f({\boldmathξ},{\boldmathϵ}). (3)

Thus, the noise-free estimator can be estimated in terms of measurements and uncertainties by the noise-unbiased estimator , where and, by definition, . The term is derived first by computing and then by replacing terms depending on in with terms as a function of which satisfy the requirement (see App. A). A property often used in the following sections is that a noise-unbiased linear combination of estimators is equivalent to the linear combination of noise-unbiased estimators:

 [N∑i=1ciTi(x)]∗=N∑i=1ciT∗i(x), (4)

where the coefficients are independent of the measurements .

## 4 Gaussian-noise unbiased sample moments and cumulants

Weighted sample central moments unbiased by Gaussian uncertainties, such as the variance , skewness , kurtosis and the respective cumulants about the weighted mean are derived assuming independent measurements , uncertainties and weights , as described in full detail in App. A. They are defined as follows:

 m∗2=m2−1V1n∑i=1wiϵ2i(1−wiV1)=k∗2 (5) m∗3=m3−3V1n∑i=1wiϵ2i(xi−¯x)(1−2wiV1)=k∗3 (6) m∗4=m4−6V1n∑i=1wiϵ2i[(xi−¯x)2(1−2wiV1)−ϵ2i2(1−2wiV1)2+m∗2wiV1]−3V41(n∑i=1w2iϵ2i)2 (7) (m22)∗=(m∗2)2−4V21n∑i=1w2iϵ2i[(xi−¯x)2−ϵ2i2(1−2wiV1)]+2V41(n∑i=1w2iϵ2i)2 (8) k∗4=m∗4−3(m22)∗. (9)

By definition, the above expressions satisfy

 ⟨m∗r⟩=1V1n∑i=1wi(ξi−%$¯ξ$)r. (10)

The unweighted forms can be obtained by substituting (for all ) and (for all ) in all terms, leading to:

 m∗2=m2−n−1n2n∑i=1ϵ2i=k∗2 (11) m∗3=m3−3(n−2)n2n∑i=1ϵ2i(xi−¯x)=k∗3 (12) m∗4=m4−6(n−2)n2n∑i=1ϵ2i(xi−¯x)2−6m∗2n2n∑i=1ϵ2i+3(n−2)2n3n∑i=1ϵ4i−3n4(n∑i=1ϵ2i)2 (13) (m22)∗=(m∗2)2−4n2n∑i=1ϵ2i(xi−¯x)2+2(n−2)n3n∑i=1ϵ4i+2n4(n∑i=1ϵ2i)2 (14) k∗4=m∗4−3(m22)∗. (15)

## 5 Gaussian-noise and sample-size unbiased moments and cumulants

The estimates of weighted central moments which are unbiased by both sample-size and Gaussian uncertainties, such as the variance , skewness , kurtosis and the respective cumulants, are defined in terms of the noise-unbiased sample estimators as follows:

 M∗2= V21V21−V2m∗2=K∗2 (16) M∗3= V31V31−3V1V2+2V3m∗3=K∗3 (17) M∗4= V21(V41−3V21V2+2V1V3+3V22−3V4)(V21−V2)(V41−6V21V2+8V1V3+3V22−6V4)m∗4+ −3V21(2V21V2−2V1V3−3V22+3V4)(V21−V2)(V41−6V21V2+8V1V3+3V22−6V4)(m22)∗ (18) K∗4= V21(V41−4V1V3+3V22)(V21−V2)(V41−6V21V2+8V1V3+3V22−6V4)m∗4+ −3V21(V41−2V21V2+4V1V3−3V22)(V21−V2)(V41−6V21V2+8V1V3+3V22−6V4)(m22)∗. (19)

The derivation of the sample-size unbiased weighted estimators is described in Rimoldini (2013b). The corresponding unweighted forms can be achieved by direct substitution for all , leading to:

 M∗2= nn−1m∗2=M2−1nn∑i=1ϵ2i=K∗2 (20) M∗3= n2(n−1)(n−2)m∗3=M3−3n−1n∑i=1ϵ2i(xi−¯x)=K∗3 (21) M∗4= n(n2−2n+3)(n−1)(n−2)(n−3)m∗4−3n(2n−3)(n−1)(n−2)(n−3)(m22)∗ (22) K∗4= n2(n+1)(n−1)(n−2)(n−3)m∗4−3n2(n−2)(n−3)(m22)∗. (23)

## 6 Special cases

If weights are related to measurement errors as , the noise-unbiased weighted sample moments and cumulants reduce to the following expressions:

 m∗2=m2−n−1V1=k∗2 (24) m∗3=m3−3V1n∑i=1(xi−¯x)=k∗3 (25) m∗4=m4−6V1n∑i=1[(xi−¯x)2−ϵ2i2]+6m∗2V1−3V21 (26) (m22)∗=(m∗2)2−2(m∗2+m2)V1 (27) k∗4=m∗4−3(m22)∗. (28)

In the case of constant errors, i.e., for all , some of the unweighted estimators are equivalent or similar to their noise-unbiased counterparts:

 Skewness: k∗3=k3~{}~{}and~{}~{}K∗3=K3~{}~{}~{}(also m∗3=m3~{}~{}and~{}~{}M∗3=M3), (29) Kurtosis: k∗4≈k4~{}~{}and~{}~{}K∗4=K4, (30)

where the approximation holds for large values of or ratios since

 k∗4−k4k22=6ϵ20(k∗2+k2)nk22≈6[1+2(S/N)2]n[1+(S/N)2]2, (31)

considering that, for constant errors, and . For sample cumulants up to the fourth order, only the variance depends strongly on noise. However, this is an important estimator because it is often involved in definitions of standardized skewness ( and ) and kurtosis ( and ) as follows:

 g1=k3/k3/22, G1=K3/K3/22, (32) g2=k4/k22, G2=K4/K22. (33)

For consistency with the above definitions, the noise-unbiased equivalents are defined as

 g∗1=k∗3/(k∗2)3/2, G∗1=K∗3/(K∗2)3/2, (34) g∗2=k∗4/(k∗2)2, G∗2=K∗4/(K∗2)2, (35)

although the truly noise-unbiased expressions should have been computed on the ratios in Eqs (32)–(33). The application of Eqs (34)–(35) should generally be restricted to larger samples (e.g., ) with ratios greater than a few, in order to avoid non-positive values of or .

## 7 Estimators as a function of signal-to-noise ratio

Noise-biased and unbiased estimators are compared as a function of signal-to-noise ratio with simulated data and different weighting schemes for specific signals, sampling and error laws. The values of the population moments of the continuous simulated periodic ‘true’ signal are computed averaging in phase as follows:

 μr=12π∫2π0[ξ(ϕ)−μ]rdϕ,      where   μ=12π∫2π0ξ(ϕ)dϕ. (36)

### 7.1 Simulation

Simulated signals are described by a sinusoidal function to the fourth power, which has a non-zero skewness and thus makes it possible to evaluate the precision and accuracy of the skewness standardized by the estimated variance without simply reflecting the accuracy of the variance. The level is evaluated by the ratio of the standard deviation of the true signal and the root of the mean of squared measurement uncertainties (assumed independent of the signal). The signal is sampled and 1000 times at phases randomly drawn from a uniform distribution, while the ratio varies from 1 to 1000 and determines the uncertainties of measurements as follows:

 ξ(ϕ)=Asin4ϕ (37) xi∼N(ξi,ϵ2i)         for~{}~{}ξi=ξ(ϕi)~{}~{}and~{}~{}ϕi∼U(0,2π) (38) ϵ2i=(1+ρi)μ2/(S/N)2           for~{}~{}ρi∼U(−0.8,0.8), (39)

where the -th measurement is drawn from a normal distribution of mean and variance . The latter is defined in terms of a variable randomly drawn from a uniform distribution so that measurement uncertainties vary by up to a factor of 3 for a given and ratio. Simulations were repeated times for each ratio (for and 1000).

The dependence of weighted estimators on sample size and the corresponding unbiased expressions were presented in Rimoldini (2013b). Herein, only large sample sizes are employed so that sample-size biases are negligible with respect to the ones resulting from small ratios. A sample signal and simulated data are illustrated in Fig. 1 for and . The reference population values of the mean, variance, skewness and kurtosis of the simulated signal are listed in table 1 of Rimoldini (2013b).

Error weights are defined by , while mixed error-phase weights follow Rimoldini (2013a), assuming phase-sorted data:

 wi=h(S/N|a,b)w′i∑nj=1w′j+[1−h(S/N|a,b)]ϵ−2i∑nj=1ϵ−2j       ∀i∈(1,n) (40) w′i=ϕi+1−ϕi−1                          ∀i∈(2,n−1) (41) w′1=ϕ2−ϕn+2π (42) w′n=ϕ1−ϕn−1+2π (43) h(S/N|a,b)=11+e−(S/N−a)/b        for  a,b>0. (44)

Weighting effectively decreases the sample size, since more importance is given to some data at the expense of other ones and results depend mostly on fewer ‘relevant’ measurements (e.g., weighting by the inverse-squared uncertainties can worsen precision at high levels). Weighted procedures are desirable when the dispersion and bias of estimators from an effectively reduced sample size are smaller than the improvements in precision and accuracy (e.g., weighting by inverse-squared uncertainties can improve both precision and accuracy at low ratios). Also, weighting might exploit correlations in the data to improve precision, as it is shown employing phase weights (Rimoldini, 2013a). Since correlated data do not satisfy the assumptions of the expressions derived herein, their application might return biased results. However, small biases could be justified if improvements in precision are significant and, depending on the extent of the application, larger biases could be mitigated with mixed weighting schemes, such as the one described by Eqs (40)–(44).

Estimators derived herein assume a single weighting scheme and combinations of estimators (like the variance and the mean in the standardized skewness and kurtosis) are expected to apply the same weights to terms associated with the same measurements. The function constitutes just an example to achieve a mixed weighting scheme: tuning parameters offer the possibility to control the transition from error-weighted to phase-weighted estimators (in the limits of low and high , respectively) and thus reach a compromise solution between precision and accuracy for all values of , according to the specific estimators, signals, sampling, errors, sample sizes and their distributions in the data.

### 7.2 Results

The results of simulations are illustrated for sample estimators, since the conclusions in Rimoldini (2013b) suggested that phase-weighted sample estimators can be more accurate and precise than the sample-size unbiased counterparts in most cases, especially for large sample sizes as considered herein.

Figure 2 illustrates the sample mean in the various scenarios considered in the simulations: sample sizes of and 1000, unweighted and with different weighting schemes (error-weighted, phase-weighted and combined error-phase weighted). While accuracy is the same in all cases, the best precision of the mean is achieved employing phase weights (including the low end, unlike other estimators).

Figures 316 compare noise-biased (‘uncorrected’ ) and noise-unbiased (‘corrected’ ) estimators as a function of , evaluating the following deviations from the population values:

 m2/μ2−1         vs  m∗2/μ2−1, (45) m3/μ3/22−γ1     vs  m∗3/μ3/22−γ1,      g1−γ1             vs  g∗1−γ1, (46) m4/μ22−3−γ2  vs  m∗4/μ22−3−γ2,   m4/m22−3−γ2  vs  m∗4/(m∗2)2−3−γ2, (47) k4/μ22−γ2        vs  k∗4/μ22−γ2,         g2−γ2              vs  g∗2−γ2, (48)

in both weighted and unweighted cases, for and 1000. The dependence on is described in more details in Rimoldini (2013b). Estimators standardized by both true and estimated variance are presented to help interpret the behaviour of the ratios from their components.

All figures confirm that ‘corrected’ and ‘uncorrected’ estimators have similar precision and accuracy at high levels (typically for ). Noise-unbiased estimators are found to be the most accurate in all cases and over the whole range tested. Their precision is generally similar to the noise-uncorrected counterparts, apart from estimators standardized by the estimated variance, such as , and , for which the uncorrected version can be much more precise (although biased) for , typically. As expected, the precision of estimators employing measurements per sample was greater than the one obtained with sample sizes of .

Weighting by the inverse of squared measurement errors made the estimators slightly less precise at high ratios, but more precise and accurate at low levels (except for the mean).

Weighting by phase intervals led to a significant improvement in precision of all estimators in the limit of large ratios and a reduction of precision at low (apart from the case of the mean). Tuning parameters such as and in Eq. (40) were able to mitigate the imprecision at low reducing to the error-weighted results, which appeared to be the most accurate and precise in the limit of low ratios (in these simulations). This solution might provide a reasonable compromise between precision and accuracy of all estimators, at least for .

Figures 58 show that the skewness moment is quite unbiased by noise, while the standardized version is underestimated at high because of the overestimated variance (as shown in Figs 34). While the accuracy of deteriorates at low , its precision is much less affected by noise.

The kurtosis moment (Figs 912) is less precise and accurate than the noise-unbiased equivalent, and its normalization by the squared variance reduces dramatically its inaccuracy and imprecision (since and exhibit a similar trend as a function of ). The kurtosis cumulant , instead, is much closer to its noise-unbiased counterpart, as shown in Figs 1316. The normalization of by the squared variance improves its precision at the cost of lower accuracy for : the bias of is similar to (greater than) the precision of for ().

The lower the level is, the less precise estimators are and the noise-unbiased variance can be underestimated (and even become non-positive). Thus, the skewness and kurtosis estimators standardized by or , as in Eqs (34)–(35), should be avoided in circumstances that combine small sample sizes (up to a few dozens of elements) and low S/N ratios (of the order of a few or less).

Figures related to moments and cumulants of irregularly sampled sinusoidal signals are very similar to the ones presented herein, with the exception of , which would have a similar precision but with no bias, as a consequence of the null skewness of a sinusoidal signal (since the mean of estimates is zero, they are not biased by the standardization with an overestimated noise-biased variance).

From the comparison of noise-biased and unbiased estimators with different weighting schemes, it appears that, for large sample sizes, noise-unbiased phase-weighted estimators are usually the most accurate for (apart from the special cases of standardized skewness and kurtosis when their true value is zero). For noisy signals (e.g., ), error weighting seems the most appropriate, at least with Gaussian uncertainties, thus noise-unbiased error-phase weighted estimators can provide a satisfactory compromise in general. Further improvements might be achieved by tuning parameters better fitted to estimators and signals of interest, in view of specific requirements of precision and accuracy.

## 8 Conclusions

Exact expressions of noise-unbiased skewness and kurtosis were provided in the unweighted and weighted formulations, under the assumption of independent data and Gaussian uncertainties. Such estimators can be particularly useful in the processing, interpretation and comparison of data characterized by low regimes.

Simulations of an irregularly sampled skewed periodic signal were employed to compare noise-biased and unbiased estimators as a function of in the unweighted, inverse-squared error weighted and phase-weighted schemes. While noise-unbiased estimators were found more accurate in general, they were less precise than the uncorrected counterparts at low ratios. The application of a mixed weighting scheme involving phase intervals and uncertainties was able to balance precision and accuracy on a wide range of levels. The effect of noise-unbiased estimators and different weighting schemes on the characterization and classification of astronomical time series is described in Rimoldini (2013a).

## Acknowledgments

The author thanks M. Süveges for many discussions and valuable comments on the original manuscript.

## References

• Bowley (1920) Bowley A.L., 1920, Elements of Statistics, Charles Scribner’s Sons, New York
• D’Agostino (1986) D’Agostino R.B., 1986, Goodness-of-fit techniques, D’Agostino & Stephens eds., Marcel Dekker, New York, p. 367
• Groeneveld & Meeden (1984) Groeneveld R.A., Meeden G., 1984, The Statistician, 33, 391
• Heijmans (1999) Heijmans R., 1999, Statistical Papers, 40, 107
• Hosking (1990) Hosking J.R.M., 1990, J. R. Statist. Soc. B, 52, 105
• Moors et al. (1996) Moors J.J.A., Wagemakers R.Th.A., Coenen V.M.J., Heuts R.M.J., Janssens M.J.B.T., 1996, Statistica Neerlandica, 50, 417
• Rimoldini (2013a) Rimoldini L., 2013a, preprint (arXiv:1304.6616)
• Rimoldini (2013b) Rimoldini L., 2013b, preprint (arXiv:1304.6564)
• Stuart & Ord (1969) Stuart A., Ord J., 1969, Kendall’s Advanced Theory of Statistics, Charles Griffin & Co. Ltd, London

## Appendix A Derivation of noise-unbiased moments

The derivations presented in this Appendix involve weighted estimators under the assumption of independent measurements, uncertainties and weights. Definitions and some of the relations often employed herein are listed below.

• For brevity, , and and are implied to involve all (from the 1-st to the -th) terms, unless explicitly stated otherwise.

• The following integral solutions are often employed:

 ⟨xsi⟩=∫∞−∞x′is√2πϵiexp[−(x′i−ξi)22ϵ2i]dx′i=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ξi for s=1ξ2i+ϵ2i for s=2ξ3i+3ξiϵ2i for s=3ξ4i+6ξ2iϵ2i+3ϵ4i for s=4. (49)
• The expected value of a generic estimator of independent data with Gaussian uncertainties is computed as follows

 ⟨m⟩ =∫Rnm(x′)n∏i=11√2πϵiexp[−(x′i−ξi)22ϵ2i]dnx′ (50) =∏h∫∞−∞m(x′)1√2πϵhexp[−(x′h−ξh)22ϵ2h]dx′h (51) =∑iai∫∞−∞x′is√2πϵiexp[−(x′i−ξi)22ϵ2i]dx′i ∑j≠ibj∫∞−∞x′jt√2πϵjexp[−(x′j−ξj)22ϵ2j]dx′j × ×∑k≠i,jck∫∞−∞x′ku√2πϵkexp[−(x′k−ξk)22ϵ2k]dx′k∑l≠i,j,kdl∫∞−∞x′lv√2πϵlexp