Uncertainty Quantification in Sunspot Counts

On white-light images, sunspots are visible as dark areas. They correspond to regions of locally enhanced magnetic field and act as an indicator of changing solar activity over time. They have been observed and counted since the invention of the telescope at the beginning of the seventeenth century. As such, the counting of sunspots constitutes one of the "longest-running scientific experiment[s]" (Owens 2013). In 1848, J. R. Wolf from Zürich Observatory created an index, denoted N_c, of solar activity by summing up the total number of sunspots N_s with 10 times the total number of sunspot groups N_g on a daily basis:

$$\begin{eqnarray}&&{N}_{c}=10{N}_{g}+{N}_{s}.\end{eqnarray} \tag{ 1 }$$

Figure 1 displays smoothed averages of the median value of these three quantities across a set of 21 observatories (also called "stations") chosen for the present study (see Section 2). Modeling the statistics of N_c is far from trivial, as this quantity jumps from 0 to 11 when a sunspot appears on the Sun (N_s = 1, N_g = 1, and thus N_c = 11). By construction, each active region appears thus twice in N_c. The multiplication factor in Equation (1) was introduced by J. R. Wolf to put the number of groups on the same scale as the number of spots. Indeed, during this historical period, a group contained on average 10 spots (Izenman 1985). Note that in recent solar cycles the average number of spots per group is rather around six.

Zoom In Zoom Out Reset image size

The index N_c, or rather the formula behind it, is at the basis of the ISN. The ISN is distributed through the World Data Center Sunspots Index and Long-term Solar Observations (WDC-SILSO).³ The N_c values from each observing station in the SILSO network are collected and rescaled, i.e., multiplied by a factor k, to compensate for their differing observational qualities. The N_c values are then combined on a monthly basis to produce the ISN (Clette et al. 2007), which constitutes the international reference for modeling solar activity over the long term. Despite the fact that it is arguably the most intensely used times series in all of astrophysics (Hathaway 2010), its historical part suffers from a number of errors and inconsistencies that have been partly addressed by the recalibration of the ISN in 2015 (Clette & Lefèvre 2016). Even the most recent part (1981–now) lacks proper error modeling and uncertainty quantification; see Section 1.2 below.

In order to place our work in context, we first describe how the ISN is currently obtained at the WDC-SILSO center. N_s and N_g are entered through an interface (www.sidc.be/WOLF) and stored in a database. The main processing is described in Clette et al. (2007), and we summarize a more recent version of it in Figure 2. The processing uses a pilot station, here the Locarno station, as a reference. It compares the values obtained by a station i to the pilot station via a scaling factor k_i, often referred to as the "k-coefficient":

$$\begin{eqnarray}&&{k}_{i}(t)=\displaystyle \frac{\mathrm{pilot}(t)}{{Y}_{i}(t)},\end{eqnarray} \tag{ 2 }$$

where Y_i(t) is the composite index of station i, observed at time t (expressed in days), and pilot(t) is the value of the pilot station. The monthly scaling factors are computed from a sigma-clipping mean of Equation (2), i.e., values differing by more than two standard deviations from the mean are eliminated from the computation process. This processing still suffers from its historical heritage, summarized in Table 1 of Dudok de Wit et al. (2016). For example, this table shows that between 1926 and 1981 (when the sunspot collection center was in Zurich) there were several standard observers, and no pilot station. As an index derived from count data, N_c (or N_s and N_g) does not necessarily follow a Gaussian distribution (Vigouroux & Delache 1994; Usoskin et al. 2003; Dudok de Wit et al. 2016). A processing based on sigma-clipping is thus not fully adapted, but still undoubtedly better than what was done during the Zurich era. Finally, some steps in the processing date back from the mid-nineteenth century (when J. R. Wolf introduced the sunspot index) and have not been upgraded when the collection and preservation center was moved from Zurich to Brussels in 1981. There were two reasons for this: (1) the new curators of the ISN wanted to keep the uniformity of the series, and (2) the numerical tools available at that time were limited.

Zoom In Zoom Out Reset image size

The WDC-SILSO team is currently working on improving the ISN computation and coordinates an important community effort to correct past errors. Such effort includes, among others, work by a team from the International Space Science Institute (ISSI) on recalibration of the SN,⁴ organization of sunspot workshops,⁵ and editorial work for a Solar Physics topical issue on SN recalibration (Clette et al. 2016).

Long-term analyses started with models of the shape of the sunspot number time series (Stewart & Panofsky 1938; Stewart & Eggleston 1940). They pursued the works by M. Waldmeier himself (Waldmeier 1939), who tried to understand the solar cycle and predict upcoming cycles. Later on, Morfill et al. (1991) investigated the short-term dynamical properties of the SN series using a Poisson noise distribution superimposed on a mean cycle variation. Vigouroux & Delache (1994) also use a Poisson distribution to approximate the dispersion of daily values of the SN at different regimes of solar activity. Usoskin et al. (2003) develop a reconstruction method for sparse daily values of the SN and model the monthly number of groups corresponding to a certain level of daily values by a Poisson distribution. Schaefer (1997) emphasizes the need for error bars on the AAVSO sunspot series⁶ (Foster 1999), and more recent results in Dudok de Wit et al. (2016) present a first uncertainty analysis of the short-term error, through time domain errors and dispersion errors among observing stations, still assuming a Poisson distribution. In Dudok de Wit et al. (2016), however, the authors uncover the presence of overdispersion in the SN and approximate the SN by a mix of a Poisson and a Gaussian distribution in an additive framework. Although non-Poissonian, this additive model fails to capture some of the characteristics of sunspot data. Chang & Oh (2012), on the other hand, use a multiplicative model to simulate sunspot counts in view of assessing the dependency of correction factors on the solar cycle.

Our goal in this work is to go beyond the above-mentioned historical heritage by developing a comprehensive uncertainty quantification model for the count data N_s, N_g, and N_c. These quantities are subject to different types of errors and do not behave exactly like Poisson random variables: (1) they experience more dispersion than the Poisson distribution, (2) they are not independent from one day to another (since sunspots can last from several minutes to several months on the Sun), and (3) they exhibit a large number of zeros owing to periods of minimal solar activity.

Our contribution is twofold. First, we develop robust estimators for the physical solar signal, denoted "true" signal in this paper, underlying N_s, N_g, and N_c. We propose a model for their densities that takes into account characteristics such as overdispersion and large number of zero counts. Our processing and estimators are robust to missing values and do not require filling in missing observations, contrarily to previous studies.

Second, we propose an uncertainty model that is motivated by first studies in Dudok de Wit et al. (2016) and that works within a multiplicative framework. Our model distinguishes three error types. The short-term error accounts for counting errors and variable seeing conditions from one station to another (e.g., weather, atmospheric turbulence), whereas the long-term error provides an overall bias in the number of spots (e.g., gradual ageing of the instrument or observer). Finally, a third error type aims at modeling inaccuracies occurring at solar minima and helps differentiating true from false zero counts. As an illustration, the short-term variations coming from the solar variability and the observational errors are clearly visible in Figure 1, superimposed on the approximate seasonality of the 11 yr solar cycle.

In future prospects, our work paves the way for a more robust definition of the ISN. Indeed, the analysis of the different error types allows studying the stability of observatories involved in the computation of ISN. Our study lays the ground for a future monitoring of all active stations within the SILSO network in quasi-real time, with the aim of (1) defining a stable reference of the network, (2) alerting the observers when they start deviating from the network, and (3) preventing large drifts from occurring. A new ISN could then be defined from a stable reference rather than a single pilot station and could benefit from the robust estimators and procedures (including rescaling of observing stations; see Section 4) developed in this work.

Our paper is structured as follows. Section 2 introduces the data set considered. The uncertainty model is presented in Section 3, while Section 4 details the preprocessing of the data. Section 5 provides the estimators (or proxy) for the "true" solar signal underlying N_s, N_g, and N_c, as well as their densities. Finally, Section 6 displays our results on quantification of the different error types, as well as a first stability analysis that takes into account both short- and long-term variability.

The observed counts are studied in two distinct situations: when there are sunspots (s > 0) and when there are none (s = 0). This separation is motivated by the idea that the absence or the presence of sunspots is led by a series of phenomena involving complex dynamo processes in the solar interior, and which can be modeled by a latent variable with two states. This analysis will lead to a better understanding of the observations and allows differentiating the "true" zeros of the counting process from the "false" zeros that occur when a station reports zero sunspot count in the presence of one or more spots on the Sun.

Let Y_i(t) represent either the number of spots, groups, or composite actually observed (raw, unprocessed data). The index 1 ≤ i ≤ N denotes the station, and 1 ≤ t ≤ T is the time. The probability density function (pdf) of $Y:= {Y}_{i}(t)$ may be decomposed as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{P}}(Y=0) & = & \mathop{\overbrace{P(Y=0| s(t)\gt 0)P(s(t)\gt 0)}}\limits^{1}\\ & & +\mathop{\underbrace{P(Y=0| s(t)=0)P(s(t)=0)}}\limits_{2}\\ {\boldsymbol{P}}(Y\geqslant y) & = & \mathop{\overbrace{P(Y\geqslant y| s(t)\gt 0)P(s(t)\gt 0)}}\limits^{3}\\ & & +\mathop{\underbrace{P(Y\geqslant y| s(t)=0)P(s(t)=0)}}\limits_{4}\ \mathrm{for}\ y\gt 0.\end{array}\end{eqnarray} \tag{ 3 }$$

Terms "1" and "3" in Equation (3) represent the short-term error in the presence of one or more sunspots. Term "1" reflects a situation where no sunspots are reported while there are actually some spots on the Sun ("false" zeros or observational errors due, e.g., to a bad seeing) and leads to an excess of zeros in short-term error distribution.

Term "2" captures the "true" zeros (no sunspot and no sunspot reported), while term "4" reflects a situation where the station reports some sunspots when there are no sunspot on the Sun. Term "4" is neglected outside of solar minima periods. Together, these two terms form the distribution of the error at minima, which has an excess of "true" zeros and a tail modeling the errors of the stations and the short-duration sunspots.

Results in Dudok de Wit et al. (2016) evidence a short-term, rapidly evolving, dispersion error across the stations that accounts for counting errors and variable seeing conditions. We define a similar term allowing a possible station dependence, and we denote it ₁(i, t). Assuming ${\mathbb{E}}({\epsilon }_{1}(i,t))=0$, where ${\mathbb{E}}$ is the expectation sign, our interest lies in modeling its variance and its tail to study the short-term variability of the stations.

Next, we introduce ₂(i, t) to handle station-specific long-term errors such as systematic biases in the sunspot counting process. We want to estimate its mean, μ₂(i, t), and detect whether this mean experiences sudden jumps or drifts on longer timescales.

Both ₁(i, t) and ₂(i, t) are multiplicative errors, as an observer typically makes larger errors when s(t) is higher (Chang & Oh 2012). Assembling these two types of errors, we propose the following noise model outside of solar minima:

$$\begin{eqnarray}&&{Y}_{i}(t)=({\epsilon }_{1}(i,t)+{\epsilon }_{2}(i,t))s(t)\ \mathrm{when}\ s(t)\gt 0.\end{eqnarray} \tag{ 4 }$$

Let ₃ denote the error occurring during minima of solar activity, when there exist extended periods with no or few sunspots. We assume the error ₃ to be significant when there are no sunspots (s(t) = 0) and otherwise negligible in order to not interfere with the errors ₁ and ₂. ₃ captures effects like short-duration sunspots and nonsimultaneity of observations between the stations. At solar minima, the model becomes

$$\begin{eqnarray}&&{Y}_{i}(t)={\epsilon }_{3}(i,t)\ \mathrm{when}\ s(t)=0.\end{eqnarray} \tag{ 5 }$$

Combining the three error types, we may write our uncertainty model in a compact and generic way as follows:

$$\begin{eqnarray}{Y}_{i}(t)=\left\{\begin{array}{ll}({\epsilon }_{1}(i,t)+{\epsilon }_{2}(i,t))s(t) & \mathrm{if}\ s(t)\gt 0\\ {\epsilon }_{3}(i,t) & \mathrm{if}\ s(t)=0.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

We assume the random variables (r.v.) ₁, ₂, and ₃ to be continuous, and the r.v. s, ₁, ₂, and ₃ to be jointly independent. Although the "true" number of counts s(t) is discontinuous, its product with a continuous r.v. (₁ + ₂) remains continuous. This is consistent with the fact that, after the preprocessing, the observed data Y_i(t) may be modeled by a continuous r.v.

All terms in Equation (6) exhibit an excess of zeros, that is, an unusual local peak in the density at zero, due to solar minima periods. As the solar minimum is an important part of a solar cycle, the zeros must be properly treated. Specific models such as the zero-altered (ZA) or the zero-inflated (ZI) two-part distributions may be used for this purpose (Zuur et al. 2009; Colin Cameron & Trivedi 2013). The main difference between both models is that the ZI distribution allows the zeros to be generated by two different mechanisms, contrarily to the ZA model, which treats all zeros in the same way.

As "true" and "false" zeros do not appear together in a single term of Equation (6), we find it appropriate to work with the ZA two-part model (also called the "Hurdle" model) and denote its density by f(x). In this model, the zero values are modeled by a Bernoulli distribution ${f}_{0}{(x)={b}^{1-x}(1-b)}^{x}$ of parameter b. Nonzero values follow a distribution described generically by f₁(x), either another discrete distribution (in case of modeling the counts μ_s(t) in Section 5) or a continuous distribution for ₁ and ₃ in Section 6:

$$\begin{eqnarray}f(x)=\left\{\begin{array}{ll}{f}_{0}(0)=b & \mathrm{if}\ x=0\\ (1-{f}_{0}(0))\displaystyle \frac{{f}_{1}(x)}{1-{f}_{1}(0)}=(1-b)\ \displaystyle \frac{{f}_{1}(x)}{1-{f}_{1}(0)} & \mathrm{if}\ x\gt 0.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

The ZA distribution will be used to model the estimated densities of μ_s(t) in Section 5 and those of ₁(i, t) and ₃(i, t) in Section 6.

To use Equation (6), we need an estimate of a proxy for s(t). We choose this proxy to be the mean of s(t), denoted μ_s(t). We propose as a robust estimator for μ_s(t) a transformed version of the median of the network:

$$\begin{eqnarray}&&{\hat{\mu }}_{s}(t)=T({M}_{t}),\end{eqnarray} \tag{ 9 }$$

where ${M}_{t}=\mathop{\mathrm{med}}\limits_{1\leqslant i\leqslant N}{Z}_{i}(t)$ represents the median of the network and T denotes a transformation composed of an Anscombe transform and a Wiener filtering (Davenport & Root 1968). This filtering is applied in order to clean the data from very high frequencies, which can lead to instabilities in the subsequent analysis. The generalized Anscombe transform stabilizes the variance (Murtagh et al. 1995; Makitalo & Foi 2013). It is written as

$$\begin{eqnarray}&&A(x)=\displaystyle \frac{2}{\alpha }\sqrt{\alpha x+\displaystyle \frac{3}{8}{\alpha }^{2}}.\end{eqnarray} \tag{ 10 }$$

It is commonly applied in the literature to Gaussianize near-Poissonian variables. It is needed here, as the Wiener filtering performs better on Gaussian data. Similarly to Dudok de Wit et al. (2016), pp. 14–15, we fix α = 4.2 in Equation (10). This is the optimal value found for the composite N_c. Before applying the Wiener filtering, missing values of the median of the network are imputed using the algorithm described in Dudok de Wit (2011). Only 49 values are imputed, which represents 0.2% of the total number of values on the period studied. The Wiener filtering is then applied on the transformed and complete set of median values and suppresses the highest frequencies of the signal. Finally, the imputed missing values were reset to NaN ("not a number") in ${\hat{\mu }}_{s}(t)$. The block diagram of the procedure is described in Figure 3.

Zoom In Zoom Out Reset image size

Among other tested estimators (based on the mean, the median of the network, or a subset of stations), with or without application of T, the estimator proposed in Equation (9) turns out to be the most robust to outliers.

To test the quality of our estimator ${\hat{\mu }}_{s}$, we compare it with a sunspot number extracted from satellite images of the Sun. We expect less variability when N_s, N_g, and N_c are retrieved from satellite images using automated algorithms, as the rules to count spots and groups are clearly defined. Nevertheless, the measurements are biased by these rules, and the most complex cases, e.g., at maxima, most often require either human intervention or a specific procedure in the algorithm. In any case, a measure of the "true" number of spots and groups does not exist.

As exercise for this comparison, we use the sunspot Tracking And Recognition Algorithm (STARA) sunspot catalog (Watson & Fletcher 2010), regrouping observations from 1996 May to 2010 October (solar cycle 23). This number is extracted using an automated detection algorithm from the images obtained by the MDI instrument on Solar and Heliospheric Observatory. It has a lower scaling than our estimator for the number of spots, ${\hat{\mu }}_{{N}_{s}}$, as expected since the definition of a spot in the detection algorithm excludes the pores (spots without penumbra).

We compare three quantities on the period where STARA data are available (1996–2010): N_s (STARA), ${\hat{\mu }}_{{N}_{s}}$, and N_s as recorded by the Locarno station. These are shown in Figure 4. We test the level of variability by computing the mean value of a moving standard deviation over a window of 81 days. It is equal to 14.07 for N_s (STARA) rescaled on ${\hat{\mu }}_{{N}_{s}}$ (5.38 for N_s (STARA) without scaling) against 15.68 for ${\hat{\mu }}_{{N}_{s}}$ and 27.13 for Locarno. As expected, N_s (STARA) experiences less variability than a single station, but its variability is comparable to the one of our estimator.

Zoom In Zoom Out Reset image size

Nevertheless, satellite images of the Sun have only been available since 1980, and data extracted from those images cannot be traced back until the seventeenth century. Gathering space observations during several decades also requires the use of different satellites and instruments, as instruments age in space. These instruments need calibrations that create additional inaccuracies to the extracted numbers. We thus conclude that ${\widehat{\mu }}_{s}(t)$ is a more robust estimator of the solar activity, and it will be used as a proxy for s(t) in the sequel.

We present here the statistical modeling of the number of spots N_s and the number of groups N_g. We do this separately for each component since their physical origins are driven by different phenomena: the groups convey information about the dynamo-generated magnetic field in the solar interior, whereas the emergence of individual spots would rather come from fragmented surface flux and agglomeration of small magnetic fields (Thomas & Weiss 2008). Together, the analysis of the spots and groups helps us to better understand the composite N_c and the solar activity, which is not satisfactorily described by only one of the two numbers (Dudok de Wit et al. 2016). In the remainder of this paper, we define a specific notation for the generic ${\hat{\mu }}_{s}(t)$ from Equation (9): ${\hat{\mu }}_{{N}_{s}}(t)$ for the number of spots, ${\hat{\mu }}_{{N}_{g}}(t)$ for the number of groups, and ${\hat{\mu }}_{{N}_{c}}(t)$ for the composite.

The authors in Dudok de Wit et al. (2016) showed that the numbers of spots and groups experience more overdispersion than actual Poisson variables. In order to estimate how far the distribution of the "true" s(t) departs from a Poisson distribution, we regress the conditional variance $\mathrm{Var}({{\rm{Z}}}_{i}({\rm{t}})| {\hat{\mu }}_{s}({\rm{t}})=\mu )$ versus the conditional mean ${\mathbb{E}}({Z}_{i}(t)| {\hat{\mu }}_{s}(t)=\mu )$ by OLS; see Figure 5. Whereas in a Poisson context the slope of the fit should be close to 1, for N_s > 10, our fit shows a slope of 1.5, with confidence interval (CI) CI_95% = [1.48, 1.51]. This points to overdispersion and the need for a generalization of a Poisson pdf. On the contrary, the same plot for N_g > 0 displays a slope of 0.96, with CI_95% = [0.93, 0.99], confirming the validity of a Poisson process assumption. Note that values <11 are excluded from the fit of N_s, as they seem to indicate a different regime. This change in the alignment may indicate the presence of a multimodal distribution; see Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Count data with overdispersion are widely modeled by the negative binomial (NB) distribution in the literature (Colin Cameron & Trivedi 2013; Rodriguez 2013) or by its generalization (Jain & Consul 1970):

$$\begin{eqnarray}&&g(x,r,p)=\displaystyle \frac{{\rm{\Gamma }}(r+x)}{{\rm{\Gamma }}(r){\rm{\Gamma }}(x+1)}{p}^{r}{q}^{x},\end{eqnarray} \tag{ 11 }$$

where r > 0, p ∈ (0,1), q = (1 − p) and Γ is the gamma function.

A visual inspection of the histogram of estimated values ${\widehat{\mu }}_{{N}_{s}}(t)$ in the left panel of Figure 6 reveals a local maxima in the distribution around 20–40. We refer to these local maxima as modes in the remainder of the article. The underlying density of ${\hat{\mu }}_{{N}_{s}}(t)$ may thus be multimodal, as suspected from the left panel of Figure 5. Such pdf's are classically modeled by a mixture model. As the density shows a typical excess of zeros as well, it requires the use of a ZA distribution defined in Equation (7). We thus fit the complete pdf of the estimated number of spots, ${\hat{\mu }}_{{N}_{s}}(t)$, by a ZA mixture of generalized NB distributions. The density at zero, f₀(x), is represented by a Bernoulli distribution, whereas the density outside zero, f₁(x) in Equation (7), is identified by a mixture of NB distributions:

$$\begin{eqnarray}&&{f}_{1}(x,{r}_{1},{r}_{2},{p}_{1},{p}_{2})={w}_{1}{g}_{1}(x,{r}_{1},{p}_{1})+(1-{w}_{1}){g}_{2}(x,{r}_{2},{p}_{2}),\end{eqnarray} \tag{ 12 }$$

where g₁, g₂ are NB densities and w₁ is the mixture weight.

Similarly, the histogram of ${\hat{\mu }}_{{N}_{g}}(t)$ exhibits a clear excess in the range of 1–3 compared to a Poisson-like distribution centered between 5 and 8. The pdf of ${\widehat{\mu }}_{{N}_{g}}(t)$ shows thus two modes: one around 1–3, and one around 5–8. Such a pdf may be modeled by a mixture of NB and Poisson distributions:

$$\begin{eqnarray}&&f(x,{r}_{1},{p}_{1},{\mu }_{2})={w}_{1}g(x,{r}_{1},{p}_{1})+(1-{w}_{1})\displaystyle \frac{{\mu }_{2}^{x}}{x!}{e}^{{\mu }_{2}},\end{eqnarray} \tag{ 13 }$$

where μ₂ > 0 and, as above, w₁ is the mixture weight.

The fit of these parametric densities is shown in Figure 6 by a black line superimposed on the histograms. All the fits in this article are computed using the maximum likelihood estimation (MLE). The nature of the different modes in the pdf of ${\hat{\mu }}_{{N}_{s}}$ and ${\hat{\mu }}_{{N}_{g}}$ will be discussed in Section 5.5.

We now use Equation (9) to estimate the ${\mu }_{{N}_{c}}$, the "true" value of the composite N_c = N_s + 10N_g. Again looking at the conditional mean–variance relationship, we observe in the left panel of Figure 7 an overdispersion with a slope of 1.29 and CI_95% = [1.27, 1.31] for N_c > 20. As a compound of both quantities, N_c experiences less overdispersion than N_s and more than N_g. A visual inspection of the histogram of ${\hat{\mu }}_{{N}_{c}}(t)$ values in the right panel of Figure 7 indicates an excess of zeros and several modes, most probably coming from the modes observed in the pdf's of ${\hat{\mu }}_{{N}_{g}}$ and ${\hat{\mu }}_{{N}_{s}}$. We thus find it appropriate to approximate the density of ${\hat{\mu }}_{{N}_{c}}(t)$ by a ZA mixture of three NB distributions, where the density outside zero values, f₁(x) in Equation (7), is identified with

$$\begin{eqnarray}&&{f}_{1}(x,{r}_{1},\ldots ,{r}_{3},{p}_{1},..,{p}_{3})=\displaystyle \sum _{i=1}^{3}{w}_{i}{g}_{i}(x,{r}_{i},{p}_{i}),\end{eqnarray} \tag{ 14 }$$

where w_i are the mixture weights and ${\sum }_{i=1}^{3}{w}_{i}=1$. The fit of f₁ is represented in the right panel of Figure 7 by a black line.

Zoom In Zoom Out Reset image size

A statistical analysis (not presented here) shows that the distribution of the ISN is statistically close to the distribution of ${\hat{\mu }}_{{N}_{c}}$. The uncertainty analysis for ${\hat{\mu }}_{{N}_{c}}$, presented in the remainder of the article, remains thus valid for the ISN.

Due to the physical nature of the data, the local maxima for the densities of ${\widehat{\mu }}_{{N}_{s}}$ and ${\widehat{\mu }}_{{N}_{g}}$ are not independent. We therefore look at the conditional correlation $\mathrm{Corr}({N}_{s},{N}_{g}| {\widehat{\mu }}_{{N}_{c}}=s)$ with the goal to better understand the nature of the modes observed in these two densities, and thus also in the density of ${\hat{\mu }}_{{N}_{c}}$.

Figure 8(a) shows the conditional correlation for different values of ${\hat{\mu }}_{{N}_{c}}$ between 0 and 400. Note that even when ${\hat{\mu }}_{{N}_{c}}=s$, the value of the composite N_s + 10N_g for a particular station may be larger (resp. smaller) than s. Our analysis highlights three regimes of activity:

• Minima: ${\hat{\mu }}_{{N}_{c}}\in [0,11]$. Here, the number of spots and groups oscillates between 0 and 1. As the number of spots equals exactly the number of groups, the correlation is high.
• Medium activity: ${\hat{\mu }}_{{N}_{c}}\in [12,60]$. The correlation progressively decreases, because the number of spots increases faster than the number of groups and then stabilizes. This regime is characterized by the development of smaller spots without penumbra or with a small penumbra. Figure 8(b) shows the bivariate boxplot of N_s and N_g when ${\hat{\mu }}_{{N}_{c}}=40$. For N_g = 1 or N_g = 2, we observe values of N_s as high as 40. We clearly observe groups containing a large number of spots, as well as groups, composed of fewer spots, that appear progressively as the penumbra grows and that indicate a transition toward groups with fewer but larger spots. The effect of this transition from small to larger spots is observed in Figure 5 from Clette & Lefèvre (2016).
• High activity: ${\hat{\mu }}_{{N}_{c}}\gt 60$. Figure 8(c) shows the bivariate boxplot of N_s and N_g when ${\hat{\mu }}_{{N}_{c}}=70$. The plot has a potato shape around (N_s = 30, N_g = 4). We now observe all kinds of groups. The correlation between groups and spots slightly increases as the number of groups begins to grow as well.

Zoom In Zoom Out Reset image size

The medium and high regimes are reflected in the estimated densities of ${\widehat{\mu }}_{{N}_{s}}$ and ${\widehat{\mu }}_{{N}_{g}}$ in Figure 6. The first mode of ${\hat{\mu }}_{{N}_{g}}$, ranging from 1 to 3, corresponds to the medium regime, while the second mode, ranging from 5 to 8, reflects the high regime. The two distinct regimes provide another justification for the use of a multimodal distribution to characterize the pdf of ${\hat{\mu }}_{{N}_{g}}$. Similarly, there is also a mode in the distribution of ${\widehat{\mu }}_{{N}_{s}}$ around 20–40 that comes from the transition between the medium and the high regime. The mode is correctly represented by a mixture model. The study of the conditional correlation constitutes the first step toward retrieving the distribution of N_c from its composites N_s and N_g. However, this task is challenging and goes beyond the scope of the article because (1) the distributions of N_s and N_g are complex mixtures and (2) the number of spots is nontrivially correlated to the number of groups.

The study of solar minima periods is complex, as the data show a large variability and dichotomy. Observed values of the error at minima, ₃, are defined as counts made by the stations when the proxy for s(t), defined in Equation (9), is equal to zero:

$$\begin{eqnarray}&&{\hat{\epsilon }}_{3}(i,t)={Z}_{i}(t)\ \mathrm{when}\ {\hat{\mu }}_{s}(t)=0,\end{eqnarray} \tag{ 15 }$$

where Z_i(t) corresponds to N_s, N_g, or N_c, and where the generic ${\widehat{\mu }}_{s}(t)$ has to be replaced by ${\widehat{\mu }}_{{N}_{s}}(t)$ for N_s, ${\widehat{\mu }}_{{N}_{g}}(t)$ for N_g, and ${\widehat{\mu }}_{{N}_{c}}(t)$ for N_c.

A visual inspection of the histogram of N_s (resp. N_g) in the left (resp. middle) panel of Figure 9 shows an important amount of "true" zeros together with two modes around one and two. Similar modes occur around 11 and 22 in the distribution of N_c in the right panel of Figure 9, as expected. These modes represent short-duration sunspots. Due to the nonsimultaneity of the observations between stations, the proxy for s(t) might be equal to zero even if some spots appear shortly (from several minutes to several hours) on the Sun. These modes can be represented by a t-location-scale (t-LS) distribution, which is a generalization of the Student t-distribution. This distribution has three parameters to accommodate for asymmetry and heavy tails: the location μ, scale σ > 0, and shape ν > 0 (see Taylor & Verbyla 2004; Evans et al. 2000). Its pdf is defined as

$$\begin{eqnarray}&&g{(x,\mu ,\sigma ,\nu )}_{t-\mathrm{LS}}=\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{\nu +1}{2}\right)}{\sigma \sqrt{\nu \pi }{\rm{\Gamma }}\left(\tfrac{\nu }{2}\right)}{\left(\displaystyle \frac{\nu +\tfrac{{\left(x-\mu \right)}^{2}}{{\sigma }^{2}}}{\nu }\right)}^{-\left(\tfrac{\nu +1}{2}\right)}.\end{eqnarray} \tag{ 16 }$$

The large proportion of zero values for ${\hat{\epsilon }}_{3}$ requires the use of a ZA model as in Equation (7). We choose a ZA mixture of t-LS for the complete distribution of ${\hat{\epsilon }}_{3}$. The density outside of zero, f₁(x) in Equation (7), is thus identified by such a mixture of t-LS distributions:

$$\begin{eqnarray}&&\begin{array}{l}{f}_{1}(x,{\mu }_{1},{\sigma }_{1},{\nu }_{1},{\mu }_{2},{\sigma }_{2},{\nu }_{2})\\ \quad =\,{w}_{1}g{\left(x,{\mu }_{1},{\sigma }_{1},{\nu }_{1}\right)}_{t-\mathrm{LS}}+(1-{w}_{1})g{\left(x,{\mu }_{2},{\sigma }_{2},{\nu }_{2}\right)}_{t-\mathrm{LS}},\end{array}\end{eqnarray} \tag{ 17 }$$

where, as before, w₁ is the mixture weight. The histograms and fitted distributions for ${\widehat{\epsilon }}_{3}$ are shown in Figure 9. The visual closeness between the histogram and the fitted distribution was used as a criterion to select the best pdf among a few intuitive candidates, while the parameters of the distribution are estimated via MLE.

Zoom In Zoom Out Reset image size

In the previous figures, where the error at minima is represented for all stations combined, outliers defined as ${\hat{\epsilon }}_{3}(i,t)\gt 2$ are not visible for N_g and N_s. A separate analysis (not presented here) shows that the percentage of outliers in each station is low (inferior to 0.5% for N_s). Some stations also observed a high maximal value at minima (e.g., a value of 35 was recorded in QU [Quezon, Philippines] for N_s). This extreme value for minima may correspond to a transcription error that might be verified in the future, before being encoded in the SILSO database.

When the proxy for s(t), defined in Equation (9), is different from zero, the short-term error $\widetilde{\epsilon }$ may be estimated using

$$\begin{eqnarray}&&\widehat{\widetilde{\epsilon }}(i,t)=\displaystyle \frac{{Z}_{i}(t)}{{\hat{\mu }}_{s}(t)}\,\ \mathrm{when}\ {\hat{\mu }}_{s}(t)\gt 0.\end{eqnarray} \tag{ 18 }$$

To select the best distribution, we proceed as follows. Different densities are fitted to the values of $\widehat{\widetilde{\epsilon }}$, outside of zero, using MLE.⁸ Then, the AIC criterion is used to choose the best pdf, which in this case is a t-LS distribution.

As we observe an excess of zero, we need a ZA t-LS distribution to represent the complete distribution of $\widehat{\widetilde{\epsilon }}$. Figure 10 shows the histogram and the fitted pdf of $\widehat{\widetilde{\epsilon }}$ outside of zero. For the latter, the mean is close to 1, indicating that on average the stations are aligned with ${\hat{\mu }}_{s}(t)$. The histogram exhibits a probability mass at zero representative of "false" zeros, that is, of stations that do not observe any sunspot when there are actually some on the Sun. The histogram also shows a tail on the right-hand side, caused by outliers. This asymmetry requires a t-LS rather than a Gaussian distribution to be fitted.

Zoom In Zoom Out Reset image size

The violin plots of four different stations are shown in Figure 11 for the number of spots N_s, where the differences between the stations are the most visible. The mean of the Locarno station (LO), the current reference of the network, is slightly higher than the three other means (and higher than the means of all other stations), around 1.19. This results from its particular way of counting: large spots (with penumbra) count for more than small spots without penumbra.

Zoom In Zoom Out Reset image size

Another characteristic feature is how the error is distributed around the mean. A violin plot may be seen as a pdf with the x-axis of the density drawn along the vertical line of the boxplot. For example, the pdf of the short-term error of LO is concentrated around the mean, but the entire distribution is shifted upward, unlike the pdf of Uccle (UC), which has much lower values. UC is a professional observatory. Different observers record from one week to another the number of spots, groups, and composite on the Sun. As their experience and methodology slightly change, the shift of observers probably increases the short-term variability of the station. Usually a team of observers experience more variability than a single person, like in FU (Fujimori, Japan). This station has remarkable short-term stability.

Similarly, the San Miguel (SM) station shows the typical shape of a professional observatory. On the other hand, the LO station shows an $\widehat{\widetilde{\epsilon }}$ distribution almost characteristic of a single observer: that is because until recently LO had one dominant main observer.

A generic estimator for the long-term error ₂(i, t) may be defined by

$$\begin{eqnarray}&&{\hat{\mu }}_{2}(i,t)={\left(\displaystyle \frac{{Y}_{i}(t)}{{M}_{t}}\right)}^{\star }\,\ \mathrm{when}\ {M}_{t}\gt 0,\end{eqnarray} \tag{ 19 }$$

where the ⋆ denotes the smoothing process, Y_i(t) are the raw observations, and ${M}_{t}=\mathop{\mathrm{med}}\limits_{1\leqslant i\leqslant N}{Z}_{i}(t)$ is the median of the network. The T transform from Equation (9) is not required here, as we apply a moving average (MA) of length L defined below.

This length L should be larger than what is considered as short-term, that is, periods inferior to one solar rotation (27 days). Long-term, on the other hand, is usually defined as periods above 81 days (Dudok de Wit 2011), beyond which the effects of the solar rotation and of the sunspot's lifetime are negligible. The midterm temporal regime corresponds to periods between 27 and 81 days. To select the long-term scale for a given station i, we make the assumption that, for all t belonging to a window of length L, we have

$$\begin{eqnarray}&&{\mathbb{E}}({\epsilon }_{2}(i,t))={\mu }_{2}(i,t)\,\simeq {C}_{i},\end{eqnarray} \tag{ 20 }$$

where C_i is a constant, which might differ from 1. Having C_i = 1 means that the station i is at the same level as the median of the network. We test different lengths L (L > 27 days) for the MA window and select the long-term regime as the shortest length for which the above assumption is valid. We consider thus ${\hat{\mu }}_{2}(i,t){\rm{s}}$ of Equation (19) in sliding windows of length L over the total period (1947–2013). We apply a nonparametric equivalent of the t-test (the Wilcoxon rank sum test; Bridge & Sawilowsky 1999) on the ${\hat{\mu }}_{2}(i,t){\rm{s}}$ to test whether Equation (20) is verified within each window. Longer windows correspond thus to a stronger smoothing but also contain more values to test. As a result of this procedure, we define the long-term regime as all scales above 81 days, as we found this to be the shortest length such that the constant assumption on μ₂(i, t) is not violated more than roughly 10% of the time. This ties in with what solar physicists consider as the long-term regime.

Depending on our interest in detecting long-term drifts or jumps, different window lengths may be chosen in Equation (19) (some well above 81 days). Indeed, drifts require long smoothing periods (several months, or even years) to be observed, whereas jumps might be oversmoothed by such long smoothing and hence need a smaller MA window.

Figures 12(a)–(d) represent the long-term drifts associated with N_s in four stations starting from 1960. Figures 12(e)–(h) show the scaling factors κ_i(t₂)s for the same stations used at short-term and minima regimes. We do not represent years before 1960 because FU and SM show too few observations in that period. FU and UC appear relatively stable, unlike stations LO and SM, which display severe drifts. Bias in the counting process is also larger during solar minima when there are short-duration sunspots. This effect is clearly visible in LO. Indeed, erroneous encoding of counts leads to much higher relative errors during minima than during the remaining part of the solar cycle. Some jumps are also visible on the graphs with the smallest MA length (81 days) in green. This scale is more appropriate to observe the jumps, while longer scales only highlight the long-term drifts of the stations.

Zoom In Zoom Out Reset image size

We emphasize here the strong link between the preprocessing and the long-term analysis. Indeed, the scaling factors presented in Section 4 are a rough estimate of the long-term error, inspired by the historical procedure of J. R. Wolf. This rough estimation is required to rescale the stations to the same level. This rescaling is used to compute the median M_t of Equation (19). Contrarily to the piecewise constant κ_is computed in Section 4, the ${\hat{\mu }}_{2}(i,t){\rm{s}}$ are smooth over time and hence are more adapted to a future monitoring of the stations.

In previous sections, we presented separately the estimations of the short-term error $\hat{\widetilde{\epsilon }}(i,t)$, the long-term error ${\hat{\mu }}_{2}(i,t)$, and the error at minima ${\hat{\epsilon }}_{3}(i,t)$. All three types of errors are needed to assess the quality and stability of one station. It is more important for a station to have a low variability (low interquantile range) than to be aligned on the mean on the network. Indeed, as seen in Section 4, it is easy to rescale a station on the mean of the network.

Figure 13 displays a visual representation of long-term against short-term error for each station. It shows the long-term versus short-term empirical interquantile range on a 2D plot and thus characterizes the stability of the stations outside of minima. Stations in red are the teams of observers. They usually experience more short-term variability than an individual. We see that MO (Mochizuki, Japan), FU, and KOm (Koyama, Japan) have low variability in both the short and long term. They correspond to long individual observers with stable observation practices. On the other hand, the LO station shows a poor long-term stability, while its short-term variability is remarkably low for a professional observatory. As mentioned earlier, this is due to the fact that there is a main observer. UC shows a large variability in the short term (due to many observers) but an interesting long-term stability, as already noticed in Figure 11. SM experiences the most severe long-term variability of the network. It also has a large short-term variability, characteristic of a team of observers. QU shows a large short-term variability and a low long-term variability level. Although it seems that it is a single observer, it appears there was a move from one place to another during the observing period, and maybe a change of instrument that would impact the short-term variability. This surprising behavior will prompt SILSO to ask for more metadata.

Zoom In Zoom Out Reset image size