Fitting FFT-derived Spectra: Theory, Tool, and Application to Solar Radio Spike Decomposition
Gelu M. Nita, Gregory D. Fleishman, Dale E. Gary, William Marin, Kristine Boone
aa r X i v : . [ a s t r o - ph . S R ] J un Fitting FFT–derived Spectra: Theory, Tool, and Application toSolar Radio Spike Decomposition
Gelu M. Nita , Gregory D. Fleishman , Dale E. Gary , William Marin , and KristineBoone ABSTRACT
Spectra derived from fast Fourier transform (FFT) analysis of time-domaindata intrinsically contain statistical fluctuations whose distribution depends onthe number of accumulated spectra contributing to a measurement. The tail ofthis distribution, which is essential for separation of the true signal from the sta-tistical fluctuations, deviates noticeably from the normal distribution for a finitenumber of the accumulations. In this paper we develop a theory to properlyaccount for the statistical fluctuations when fitting a model to a given accu-mulated spectrum. The method is implemented in software for the purpose ofautomatically fitting a large body of such FFT-derived spectra. We apply thistool to analyze a portion of a dense cluster of spikes recorded by our FST instru-ment (Liu et al. 2007) during a record-breaking event (Cerruti et al. 2006) thatoccurred on 06 Dec 2006. The outcome of this analysis is briefly discussed.
Subject headings:
1. Introduction
Digitized time domain signals with a given instantaneous bandwidth are often ana-lyzed at high spectral resolution using the Fast Fourier Transform (FFT) (Benz et al. 2005a;Mannan et al. 2000; Heydt et al. 1999). As we show below, fitting FFT-derived spectra witha model function can be challenging because of their non-Gaussian statistical fluctuations,which depend on both the true signal and any sources of added noise. The exact propertiesof the fluctuations depend on the number of raw FFT spectra accumulated prior to analyzingthe spectrum. Center For Solar-Terrestrial Research, New Jersey Institute of Technology, Newark, NJ 07102, USA Binghamton University, Vestal, NY 13850, USA University of Calgary, Calgary, Canada n time-domain samples,which are converted to n -channel raw spectra and then accumulated over M such spectra toget a final output spectrum. This averaging reduces the amplitude of statistical fluctuations,but at the same time degrades the temporal or spectral resolutions. Therefore, there is atrade-off in the number of channels n and accumulations M for a given application. Fortu-nately, our FST observations offer a unique data set for studying the interplay between thespectral resolution and the level of statistical fluctuations, because the instrument directlyrecords digitized time-domain samples obtained with 1 ns resolution. These are organized in100 µ s contiguous blocks, which are separated by 19.9 ms time-domain gaps needed to readout and store the previously acquired contiguous block (Liu et al. 2007). Direct recordingof time-domain data allows one to post-process each contiguous time-domain block withdifferent FFT settings (different values of n and M ), to obtain accumulated spectra withadjustable spectral resolution and accumulation length.In this study we take advantage of this flexibility to identify the optimal spectral reso-lution needed to perform spike decomposition and obtain reliable statistical distributions ofspike parameters. To achieve this goal we developed a theory describing properties of thestatistical fluctuations for M -accumulate raw FFT spectra, tested it on simulated data, andcreated a software tool capable of automatically fitting the corresponding data. This toolhas been employed in the analysis of FST data in order to get the amplitude and bandwidthdistributions of the observed spikes, both of which turn out to be asymmetric distribu-tions with prominent power-law tails toward high amplitudes or large bandwidths. Theopen source code of this tool has been made publicly available as part of the SolarSoftWare(SSW- ) scientific software repository, and is being locatedin the SSW distribution tree at .. \\ ssw \ radio \ ovsa \ idl \ fitting \ spike explorer.pro . In sec-tion 2 we describe the general problem in using FFT data to measure reliable parametersof the observed spikes. The basic approach is to add spectral-peak model components rep-resenting individual spikes until an overall minimum χ criterion is met. We find that acorrect determination of χ requires a new theoretical treatment of the noise statistics forFFT-derived spectra, which we provide in section 3. To actually apply the theory to ourproblem requires not only estimating errors for a given model, but also determining whetherand where to place additional model components. A novel approach based on statistical 4 –likelihood of the run-length of regions of mismatch between data and model is developed insection 4. In section 5 we describe the model-fitting algorithm using simulated data. Fi-nally, the algorithm is applied to the 6 December 2006 event in section 6, and the results arediscussed in section 7.
2. Statement of the Problem
To illustrate the problem we face in this study, consider the simulated spike cluster shownin Figure 1. The problem is to construct a model that adequately fits the data, consisting ofone or more overlapping spectral components, each of which we take here as being adequatelydescribed by a Gaussian of adjustable amplitude and width. The justification for this choicewill be given in Section 5, but the general algorithm is independent of the specific form ofthe model. The figure presents three standard least-squares solutions obtained from fittingan FFT-derived spectrum with different spectral models assuming superpositions of a flatbackground and one, two, or three Gaussian peaks, i.e. s = ξ + N X k =1 α k e − (cid:16) x − βkγk (cid:17) , (1)where N = 1 , { ξ, α k , β k , γ k } (that is, background level, amplitude, positionand width), are the unknown free parameters corresponding to each component. It mayappear that there are dozens of spikes in this cluster. However, because the mean and noiseare related for FFT-derived spectra, the noise is largest in channels where the intensity isgreatest. A statistically-based algorithm designed for the specific case of FFT-derived noisestatistics is needed to assess whether a deviation of the model from the data can be explainedby noise, or instead requires placement of an additional component.The next two sections use the specialized language and notation of statistics theory toformally derive the relevant relationships necessary to address the problem outlined above.Readers interested only in the application and results can skip to section 6. 5 –Fig. 1.— Least-squares solution (thick line) for the same FFT-derived spectrum based onthree different spectral models assuming superposition (thin lines) of a flat background andone, two, and three Gaussian peaks (panels a, b, and c, respectively). 6 –
3. Statistical Properties of FFT-derived Spectra3.1. Statistical Fluctuations
Statistics of the FFT-derived spectra can be formulated in terms of the standard
Gamma probability distribution function (PDF) G ( x ; k, θ ) = x k − e − xθ θ k Γ( k ) , (2)where Γ( z ) = R ∞ t z − e − t dt is the Euler’s Gamma function, and k and θ are, respectively, theshape and scale parameters of the distribution, which determine its mean, µ = kθ , variance, σ = kθ , skewness, α = 2 / √ k , and kurtosis β = 3 + 6 /k .In particular, G ( y j , , µ j ) represents an exponential PDF of mean µ j and variance µ j ,which defines the statistical properties of a single ( j –indexed) channel of a raw FFT-derivedspectrum (Nita et al. 2007; Nita & Gary 2010a), while G ( S j , M, µ j ) represents the PDFcorresponding to the sum of M such random samples, which defines the statistical propertiesof each channel in an accumulation of M , i –indexed, consecutive raw FFT-derived spectra(Nita & Gary 2010a), S j = M X i =1 y ( i ) j . (3)Hence, the frequency–dependent population means and variances of such an accumula-tion are given by s j = M µ j and σ j = s j /M , respectively. The channel–dependent parentdistributions of the accumulated power S j may be expressed in terms of their individualstatistical means s j and their common accumulation number M as G ( S j ; M, s j /M ) = M M Γ( M ) 1 s j ρ M − j e − Mρ j , (4)where we introduce the sample-to-mean ratio ρ j = S j /s j sout to denote the ratios betweenthe random variables S j and their corresponding population means s j , which we are going toeventually derive from the data in the form of the Sample to Model Ratio (SMR) estimator, b ρ j = S j b s j , (5)with b s j being the most likely solution for s j . 7 – If a model function b s j , described by a set of yet to be determined free parameters { p k } , ( k = 1 , ν ), is a true solution for the parent population means s j corresponding to theFFT-derived spectral points S j , ( j = 1 , N ), the PDF given by equation (4) can be used tobuild the associated likelihood function (Bevington & Robinson 1992) L ( p , p , ..., p ν ) = N Y j =1 G [ S j ; M, M b s j ( p , p , ..., p ν )] , (6)which represents the conditional probability density function associated with the observation.The maximum likelihood estimates of the model parameters { p k } are obtained by max-imizing the likelihood function or, equivalently, by minimizing the more mathematicallyconvenient negative log–likelihood function (Bevington & Robinson 1992) λ ( p , p , ..., p ν ) ≡ − L ] = (7)2 M N X j =1 (cid:26) S j b s j + ln( b s j ) − (cid:18) − M (cid:19) ln ( S j ) + c ( M ) (cid:27) , where c ( M ) = ln [Γ( M )] /M + ln( M ) is a channel–independent constant. For Gaussianstatistics, λ is equivalent to the standard χ function χ = N X j =1 w j ( S j − b s j ) , (8)where w j are the statistical weights, related by χ = λ + constant. A major drawback of maximum likelihood is that it does not provide a standard methodfor assessing goodness of fit. The standard χ -based goodness-of-fit assessment may not bevalid unless the χ parameter associated with the problem follows a χ PDF with ν degreesof freedom, i.e. PDF [ χ ( ν )] = G (cid:16) χ ; ν , (cid:17) , (9)where ν = N − n is given by the difference between the number of the data samples N andthe number of free parameters n of the model function. 8 –The expectation of the random variable χ is E ( χ ) = ν or, equivalently, E ( χ ν ) = 1,where χ ν ≡ χ /ν is the reduced chi-square associated with the re-normalized distribution PDF [ χ ν ] = G (cid:18) χ ν ; ν , ν (cid:19) . (10)To investigate whether the standard χ function (8) works for FFT-derived spectra, weneglect all sources of noise other than the FFT statistical fluctuations, so that statisticalweights w i in equation (8) are described by the variance of the parent Gamma distributionof the random variables s j , i.e. w j = 1 /σ j = M/s j , which depend exactly on the quantitiesto be estimated. We estimate these statistical weights from the model itself as w j = M b s j . (11)Substituting these yet unknown weights to equation (8) we find χ I ) ≡ M N X j =1 (1 − b ρ j ) , (12)with b ρ j defined by equation (5). The superscript ( I ) is used to distinguish this expressionfor χ with another that we will develop shortly.Clearly, the solution provided by the minimization of equation (10) is located in thevicinity of the solution provided by least-squared minimization of equation (12). Indeed,both minimization problems seek a spectral model estimate b s j that globally minimizes thedeviation of b ρ j from unity. The SMR estimator b ρ j is a random variable belonging to the sameparent population as ρ j . The PDF of ρ j may be straightforwardly derived from equation (4)through a simple scale transformation, leading to G (cid:18) ρ j , M, M (cid:19) = M M e − Mρ j ρ M − j Γ( M ) , (13)with mean 1 and variance 1 /M , cf. equation (2).Making use of parent PDF (13) of ρ j , the expectation of χ I ) (12) can be computedunder the assumption that the model function b s j represents the true spectral shape s j : E (cid:2) χ I ) (cid:3) = N. (14) 9 –Using equation (14), we may define the normalized functional form with unity expecta-tion χ I ) N ≡ M N N X j =1 (1 − b ρ j ) , (15)or in the case of n model parameters to be determined from the data, we replace N with thenumber of degrees of freedom ν = N − n , as in the standard reduced-chi-square estimator(Bevington & Robinson 1992), to arrive at our final goodness-of-fit estimator χ ν ≡ M ν N X j =1 (1 − b ρ j ) . (16)As shown in Appendix A, the PDF of χ ν may be approximated by the analytical ex-pression PDF (cid:2) χ ν (cid:3) ≈ G " χ ν ; 1 (cid:0) M (cid:1) ν , (cid:18) M (cid:19) ν , (17)which, in the limit of large accumulation length M , reduces to a classic chi-squared distri-bution normalized by its degrees of freedom, equation (10). Using this approximation, theprobability to observe a given χ ν or larger may be computed as P ( χ ν ) = γ (cid:0) MM +3 ν , MM +3 ν χ ν (cid:1) Γ (cid:0) MM +3 ν (cid:1) . (18)where γ ( a, z ) = Z ∞ z t a − e − t dt (19)is the incomplete Gamma function.
In order to validate the theoretical derivations presented up to this point, and also forfurther reference, we have performed a Monte Carlo simulation in which we generated asequence of 2 million ρ j random variables distributed according to a G ( ρ j , M, /M ) PDF,which corresponds to a structureless flat spectrum of unity mean. A particular value of M = 48 has been chosen for illustration. Figure 2( a ) displays the density distribution ofthe numerically generated SMR data set, as well as a curve representing the correspondingtheoretical PDF given by equation (13). 10 –The 2 million values have been divided into 2000 contiguous blocks of length N = 1000,from which are computed a sequence of 2000 SMR means, η i = P Nj =1 ρ j . Figure 2( b ) displaysthe density distribution of the SMR means, as well as the corresponding theoretical PDFcurve G ( η, M N, /M N ).Finally, equation (16) has been used to compute a set of 2000 χ ν deviates, where ν ≡ N ,since there are no free parameters. Figure 2( c ) displays the density distribution of thecomputed χ ν values, as well as the corresponding PDF curve of equation (17).Visual inspection of Figure 2 shows good agreement between the numerically generateddistributions and their theoretical expectations, not only in the case of the exact ρ j and η i PDFs, but also in the case of the χ ν PDF approximation.
Another drawback of the maximum likelihood method is its tendency to converge towarda local extreme, rather than an absolute one. Instead, we use a least-squares minimization forthe purpose of making an initial guess of the unknown parameters. However, the weightingscheme adopted in equation 11 is no good for this purpose, since the weights come from themodel itself, which we do not initially know. To comply with the standard implementation ofthe least-squared minimization algorithm, which employs weighting coefficients w i derivedfrom measurement errors, we employ the unbiased variance estimator b σ j = S j / ( M + 1)derived in Appendix B to assign the sample-based weights w j → ( M + 1) /S j based on themeasured data samples. This is equivalent to minimizing the quantity χ II ) ≡ ( M + 1) N X j =1 (cid:18) − b ρ j (cid:19) , (20)which pursues in a reciprocal space the same goal as the model-based weighing scheme, i.e.finding the model that globally minimizes the deviations of the same local SMR parametersfrom unity. However, the uncertainties associated with the estimated parameters providedby the local curvature at the point of absolute χ minimum (Bevington & Robinson 1992), σ p k = 2 (cid:20) ∂ χ ∂p k (cid:21) − (21)are different for the two weighting schemes, which affects the most likely solution obtained,as quantitatively illustrated in Appendix C.Minimization of estimator equation (20) is the more practical of the alternatives. Con-sider the composite spectrum illustrated in Figure 3, which consists of superposition of two 11 –Fig. 2.— Monte Carlo simulation of a set of 2 million M = 48– ρ j random deviates dividedinto 2 , N = 1 , ρ j density distribution (histogram)and its theoretical PDF (solid-thick line) given by equation (13). b) Observed density dis-tribution of a sequence of 2 , η i = (1 /N ) P Nj =1 ρ j mean SMR values and its theoreticalPDF (solid-thick line), G ( η, M N, /M N ). c) Observed density distribution of 2 , χ ν val-ues computed for each contiguous bloc and its theoretical PDF (solid-thick line) given byequation (17). 12 –partially overlapping Gaussian peaks, same as in Figure 1, shown by the red curves. Inorder to fit the simulated, noise-contaminated spectrum (thin black line), calculated for anaccumulation length M = 12, one must first estimate the number of Gaussian peaks fromthe data. Figure 3(a) illustrates the fitting results with an initial guess that consists of onlya single peak. The χ I ) (blue), χ II ) (green) and λ (yellow) minimizations lead to differ-ent solutions characterized by different goodness of fit parameters χ ν and different negativelog–likelihoods λ . The λ and χ I ) minimizations lead to solutions that have goodness-of-fit and likelihoods comparable with the true model function, while the χ II ) minimizationleads to a solution with poorer goodness-of-fit and likelihood. However, it has the advan-tage of being closer to the true shape of one of the overlapping peaks, unlike the alternativesolutions that provide estimates that are closer to the envelope of the overlapping peaks .The same conclusion is supported by Figure 3(b), which displays the local deviations ρ j andtheir averages over the entire spectrum η . Now, the poorer χ II ) minimization signals theinadequacy of the model while the non-uniform pattern of the SMR deviations indicates themost appropriate location to make a model adjustment. By adding a new spectral peak tothe model close to the suggested location, the results displayed in Figures 3(c) and (d) areobtained. These panels reveal that in the case of a good model fit, χ I ) and λ minimizationsprovide more accurate estimates of the true spectral peaks than the χ II ) minimization,as indicated by their χ ν and η estimators being closer to unity. This suggests a strategywhere values from χ II ) minimization are used to identify where new model components areneeded, and χ I ) and λ minimizations are used to assess the final model parameters oncethe model components are determined. In practice, we find that λ minimization yields thebest results.
4. Fitting Composite Spectra: Estimation of the Most Likely Spectral Model
As has been shown, χ II ) -minimization is the most sensitive to a local data–modelmismatch. The presence of local SMR deviations greater than statistically expected cansuggest the spectral location at which the spectral model should be amended (by adding onemore peak, for example). This offers a starting point for designing an adaptive, self-consistentspectral-fitting algorithm that should be capable of adding complexity to an initially crudespectral model until the remaining SMR fluctuations around unity are solely due to statistical This reflects typical performance of each of these estimation methods: the χ II ) functional form gives,by design, more weight to the low amplitude data points, which are less affected by FFT noise than thepeaks, while the χ I ) and λ functional forms favor global solutions that are more consistent with the entiredata set under the assumption that the true spectral model contains only one spectral peak.
13 –Fig. 3.— Fitting results of a composite spectrum formed by superposition of two partiallyoverlapped peaks under the assumptions of a one peak model (panels a and b) and two peakmodel (panels c and d). The exact (red), χ I ) (blue), χ II ) (green) and λ (yellow) solutionsare shown in the upper panels, and their corresponding SMR deviations in the lower panels.The associated χ ν , λ and η parameters, as well as the probabilities to have observed greater χ ν under the assumption of a correct spectral model are also displayed in correspondingcolors. 14 –fluctuations, within an acceptable confidence interval. The probability to observe a local SMR exceeding an observed value ρ due to statisticalfluctuations can be computed from equation (13) to yield p ( b ρ j > ρ ) = γ ( M, M ρ )Γ( M ) . (22)Consequently, the probabilities to observe a local SMR deviation above or below the mean, p a ≡ p ( b ρ j >
1) and p b ≡ p ( b ρ j < p a = γ ( M, M )Γ( M ) (23) p b = 1 − γ ( M, M )Γ( M )To identify a strong deviation from the model, equation (22) alone may be sufficient.In a more general case, however, especially for a smaller-amplitude but spectrally-resolvedsignal, we expect that the data will deviate from an imperfect model in some compactspectral region that will result in systematic deviations of the SMR estimator either aboveor below unity.To quantitatively address this aspect of the fitting problem, we derive the PDF describ-ing the probability of observing a compact region of a given length whose SMR deviates fromunity solely as the result of statistical fluctuations. This will allow a decision on where andwhen a local improvement to the model function is needed, based on some desired probabilityof false alarms.We first derive the discrete probability mass functions (PMF) for observing a SMRcompact region of a certain size n that is located above or below unity. The process ofsuch compact regions of a given size occurring randomly can be described as a memorylessBernoulli process in which any individual SMR in a random sequence may be above orbelow unity with the unequal probabilities p a and, respectively p b (equation 23). The twosize distributions should obey complementary geometrical distributions given by f ( n a ) = (cid:20) − Γ( M ) γ ( M, M ) (cid:21) (cid:20) Γ( M ) γ ( M, M ) (cid:21) n a (24) f ( n b ) = Γ( M ) γ ( M, M ) (cid:20) − Γ( M ) γ ( M, M ) (cid:21) n b ,
15 –which represent the conditional probabilities of observing a run of SMR deviations, all on thesame side of unity, of length n a,b , terminated by any SMR deviation of the opposite sense.We have confirmed this expectation with simulated FFT-derived spectral data.In addition to the run-length criterion, we can examine the probability for an identifiedregion to deviate in amplitude by more than statistically expected from the population meanof b ρ j which we define as the mean of the n contiguous individual SMR deviations that areall above, or all below, unity, η a = 1 n j + n − X j = j [ b ρ j ≥ , (25) η b = 1 n j + n − X j = j [ b ρ j ≤ . To derive the PDFs of such means, we first obtain the conditional PDFs (CPDF) as-sociated with each type of compact deviate by splitting the unconditional PDF defined byequation (13) into two truncated
Gamma distributions, above and below unity, and normal-izing them by the factors provided by equation (23). Hence,CPDF( ρ j >
1) = Γ( M ) H ( ρ j − γ ( M, M ) M M e − Mρ j ρ M − j (26)CPDF( ρ j <
1) = Γ( M ) H (1 − ρ j )Γ( M ) − γ ( M, M ) M M e − Mρ j ρ M − j , where H ( x ) represents the Heaviside unit step function.Although the desired PDF cannot be obtained in closed form, we can neverthelesscompute its moments in a standard way, which we do in Appendix D. Equations (53) and(54) give the means and central moments µ i of the random variables η a and η b , respectively.These central moments may be used to compute, for a given set of parameters M and n ,the more commonly used standard moments, i.e. standard deviation σ = √ µ , the skewness α = µ /µ / , and kurtosis β = µ /µ . For convenience, we will also use in the subsequentderivations the alternative standard parameter β = α = µ /µ . Having determined the exact first four standard moments of the Mean SMR probabilitydistribution functions, it may be shown (see Appendix E for a detailed quantitative analysis) 16 –that the Pearson Type I PDF (Pearson 1895) is a generally valid approximation, at least forlarge run-lengths n ≫
1, for the true PDF describing the parent populations of the η a and η b random variables.This is suggested by Figure 4, which displays the relationship between the kurtosis( β ) and squared skewness ( β ) as computed from equation (53) (red + symbols) and equa-tion (54) (blue + symbols) for n = 1 , , and 9, and M ranging from 1 to 1000. The symbolsare overlaid on a color-coded diagram that indicates the parameter regions correspondingto each of the 7 probability distribution types that form the Pearson system. The PearsonType I region is upper bounded by the line, 2 β − β − β − β − n have been chosen to illustrate the influence of the two free parameters n and M on the shape of the η a and η b distributions. The three parallel symbol stripes, whichcorrespond to different values of n , indicate that the skewness range decreases as n increases,while, for a fixed n , the η a and η b distributions converge toward the same skewness as M increases. However, unless the length of a compact SMR region is fairly large, n ≫
1, itsparent PDF cannot be accurately approximated by a purely Gaussian distribution (whichcorresponds to the point on the y -axis { β = 0 , β = 3 } ).In its most general form, the Pearson Type I PDF is a generalized Beta distributionhaving the functional form f ( x, a, b, p, q ) = Γ( p + q )( x − a ) p − ( b − x ) q − Γ( p )Γ( q )( b − a ) p + q − ; (27) a ≤ x ≤ b, where a and b define the endpoints of the limited domain of definition, and p and q are twoparameters that determine the shape of the distribution. These four free parameters haveto be determined such that the first four moments of the true distributions are matched bythe first four moments of their respective Pearson Type I approximations.The solution to this problem is given by (Bowman & Shenton 2007), a = µ − σ p √ dp + qb = µ + σ q √ dp + qp = r − ( r + 2) r β d ! (28) q = r r + 2) r β d !
17 –Fig. 4.— Kurtosis ( β ) and squared skewness ( β ) as computed from Eqn. 53 (red symbols)and Eqn. 54 (blue symbols) for n = 1 , , and 9, and M ranging from 1 to 1000. The symbolsare overlaid on a color-coded diagram which indicates the parameter regions correspondingto each of the 7 probability distribution types that form the Pearson system. All Pearsonlines and regions converge in the { β = 0 , β = 3 } point, which corresponds to the normaldistribution. As M increases, the red and blue symbols corresponding to the same value of n converge from both sides of the same line toward a central point { β , β } . As n increases, thecentral point corresponding to each such line approaches the point of the diagram definingthe normal distribution. 18 – r = 6 β − β −
16 + 3 β − β d = ( r + 2) β + 16( r + 1) , where all quantities are expressed in terms of the known standard moments µ, σ, β and β .The probability for the mean of a compact SMR region, located above or below unity,to lie above or, respectively, below a certain threshold, t , may be expressed in terms of thecumulative probability function of the Beta distribution, p ( t ) = R ∞ t f ( x, a, b, p, q ) dx , whichmay be conveniently written as p ( t ) = t ≤ aI (cid:16) (1 + t σt ) , t − t t , t + t t (cid:17) a < t < b t ≥ b, (29)where I ( z, α, β ) = R z ξ α − (1 − ξ ) β − dξ R ξ α − (1 − ξ ) β − dξ (30)is the regularized incomplete Beta function, and t = (4 β − α − t − µ ) + σα ( β + 3) (31) t = q α ( β + 78 β − − β − β − α t = 5 β − α − α − β + 6 t = 3 α ( α − β + 1)( β + 3)3 α − β + 6 , where the skewness α and not √ β has been explicitly used for properly taking its sign intoconsideration.To determine whether the deviation of an observed SMR mean is statistically signifi-cant, one may use equation (29) for deviations above unity, or its complement 1 − p ( t ), fordeviations below unity. Exceeding an acceptable probability of false alarm (PFA), such asthe standard 3 σ PFA of 0 . η b = a and η a = b , respectively. 19 –In conclusion, we can combine both conditions as a multivariate probability function,the product of the probability to have a region with a given length n from equation (24) andthe probability of exceeding threshold t from equation (29), i.e. p ( t, n ) = t ≤ a < h Γ( M ) γ ( M,M ) i h − Γ( M ) γ ( M,M ) i n h − I (cid:16) (1 + t σt ) , t − t t , t + t t (cid:17)i a < t < h − Γ( M ) γ ( M,M ) i h Γ( M ) γ ( M,M ) i n h I (cid:16) (1 + t σt ) , t − t t , t + t t (cid:17)i < t < b t ≥ b, (32)where the parameters a, b, t , t , t , t and σ depend only on M , n , and the position of t relative to unity. This expression provides accurate tail probabilities associated with theSMR mean t over a compact region of length n and therefore may be used in a fittingalgorithm to decide whether to add another component to the fit model. We note that thisis a general result that does not depend in any way on the choice of the fit components. Inwhat follows, we will decompose the spectrum into Gaussian components, but one could aseasily choose wavelets, lorentzians, Voigt, or other forms.
5. An algorithm for fitting spectral data with Gaussian components
This section describes a self-adaptive algorithm specifically tailored for the problemstated in section 2—that is, estimating the most probable superposition of Gaussians thatfit an observed spectrum obtained by accumulating a known number of raw FFT-derivedspectra. Our model choice (superposition of positive-amplitude Gaussians) necessarily meansthat the algorithm described below flags only positive SMR deviations to indicate where toadd components. The theory above, however, does not make any assumption about theform of the model, and for some models (e.g. superposition of wavelets, or spectra whereboth absorption and emission are expected) it could be appropriate to flag both positive andnegative SMR deviations for adding components. The selection of positive Gaussian spectralcomponents is not arbitrary for the case of radio spike bursts. Indeed, as has been said inthe Introduction, a coherent emission mechanism of the radio spikes implies the necessity ofmany e-folding amplifications of the radiation as it propagates through the source region ofa given spike. It is apparent that the spike spectrum has a peak where the correspondingspatial growth rate κ σ has a maximum as a function of frequency (Fleishman 2004b). Sincefor the required strong amplification, only a small vicinity of this maximum plays a role, the 20 –true frequency dependence of the growth rate can be reliably approximated by a parabolicone (Fleishman 2004a) κ σ ≃ κ mσ (1 − α ( f − f ) ) , (33)where κ mσ and α are constant parameters of this approximation. The intensity of the ampli-fied waves is proportional to exp( κ σ L ), where L is the spike source size along the line of sight;thus, the parabolic dependence of the growth rate translates to the Gaussian dependence ofthe wave energy density. The ECM mechanism, assumed here, is direct amplification of theescaping electromagnetic waves, which immediately leads to a Gaussian spectral shape ofelementary spike emission. Figure 5 offers a step-by-step illustration of the algorithm for the case of a syntheticspectrum that contains a superposition of 7 peaks with different amplitudes, bandwidths,and degrees of overlap. A noise level corresponding to M = 96 has been applied to match thetypical accumulation length of the FST instrument (Liu et al. 2007) and N = 61 frequencychannels have been chosen to provide an appropriate frequency scale visualizing the spectraldetails. The first step of the algorithm performs an initial guess of a minimal number of peakspresent in the spectrum. The algorithm starts with a simple flat spectrum indicated bythe blue horizontal line in Figure 5a, which results in the parameters listed in blue on theright side of the plot. Based on this initial model function, the local SMR deviations arecomputed and displayed in Figure 5b. At this stage, the SMR sequence is divided intoadjacent compact regions deviating below or above unity, and their mean SMR deviations η b , or η a , are computed. The probabilities of observing these averaged SMR deviationsare computed according to Eqn. 29 and all SMR regions above the adopted probabilitythreshold are flagged. An initial guess of Gaussian parameters are computed based on thenumber, location, width, and average amplitude in the region. For the particular exampleused for illustration, Figure 5b indicates three such regions delimited by solid green verticallines. Their averaged SMR values are shown as blue horizontal segments, while the red solidsegments indicate the maximum range of the deviations allowed by the Pearson Type I PDFapproximations in each case. The two red horizontal lines spanning the entire frequencyrange indicate the standard 0 . n = 1 region. 21 –Fig. 5.— Fitting algorithm steps in the case of a synthetic Gaussian superposition data(panel a–red curves) contaminated by M = 96 Gamma noise. The final solution, and itsratio to the true data are shown in panels ( m ) and ( n ), respectively. At each intermediatestep, the fitted Gaussian peaks (blue lines), their corresponding envelope (yellow filling), andthe new added component (green lines) are indicated on the left-hand panels. The right-hand panels display the SMR estimator (black lines), the range of the compact regions aboveunity (vertical green dashed lines), the observed SMR means (horizontal blue segments), theirPearson Type I upper limits (red horizontal segments), and the selected least-probable regionwhere a model adjustment will be attempted in the next step (vertical green solid lines). Forreference, the horizontal red lines indicate the Pearson Type I approximation limited rangecorresponding to an n = 1; M = 96 SMR region. The inset text on the left-hand panelsgive goodness-of-fit χ ν , the negative log–likelihood λ , the averaged SMR, η , and the currentnumber of spikes n for the evolving fitting solution. 22 –The set of these three Gaussian peaks shown in Figure 5c in solid green along with the flatbackground are then used as the initial guess of the true spectral shape. χ minimization and model adjustment The next stage of the fitting algorithm involves a standard χ II ) minimization that, inour example, results in the set of three Gaussian peaks shown in Figure 5c as blue solid lines.In this and all subsequent left-hand panels, the χ ν , λ , and η parameters from the fit in thatstep are indicated in blue. Note that for overlapping peaks the χ II ) minimization findsa solution close to the low-amplitude spectral components. This built-in bias of the χ II ) minimization, which results in the recomputed SMR deviations shown in Figure 5d, offersgreater sensitivity for finding additional components as discussed earlier. Now five regions(vertical green lines) are identified, where new spectral components should be added to themodel function.At this stage of the algorithm, the flagged SMR regions are sorted in order of theirchance probability, equation (32), and the most improbable one, marked by solid verticallines in our example, is chosen as the location at which an additional spectral peak is addedto the model. Panels e - l display a sequence of identical steps that are repeated up to thepoint that, as seen in panel l , no compact SMR region deviates above the upper limit ofthe Pearson Type I PDF, and no individual SMR deviation, whether isolated or part of acompact region, crosses the n = 1 standard probability threshold. χ minimization The last stage of the fitting algorithm consists of using the final solution of the χ II ) minimization as a starting point for finding the absolute minimum of the negative log-likelihood function λ , which, as shown in Figure 3 (see also Figure 14), should be situatedin the neighborhood of the χ II ) solution. The spectral shape solution provided by the λ minimization is displayed in Figure 5m and its corresponding SMR deviations are displayedin Figure 5n. One should note that the final λ minimization results not only, as expected, ina smaller value of λ , but also in a smaller value of the goodness-of-fit parameter χ ν (whichcorresponds to the χ I ) measure).Although, in the case of the particular example illustrated in Figure 5, this final λ minimization stage of the algorithm results only in a minimal change in the final χ II ) solution, this correction has the merit of having a proper goodness-of-fit estimate and also 23 –being fully consistent with a maximum likelihood criterion. Each spectral model component has its own amplitude-dependent level of statisticalfluctuations; thus, the successful fitting algorithm must fit both the true signal (determin-istically) and the associated fluctuations (statistically). Adding spectral components to analready well-determined model will reduce the fluctuations below the expected level, whichoffers an objective means for rejecting such an excessive model. To avoid adding an exces-sive spectral component, our fitting algorithm employs a gradual augmentation of the modelfunction, which except for the starting point, adds only one spectral component at eachstep up to the point where no SMR region exceeds the preset PFA thresholds. Therefore,by design, the algorithm is not expected to continue adding spectral components above theminimally needed number.However, the χ II ) minimization might not find a valid solution for a given set of initialguess parameters, or might settle on a local rather than a global minimum. To avoid this, weimplemented a mechanism that checks the validity of solution and, in case of failure, allowsan early termination of the process before all SMR residuals are brought within the desiredlimits. For this purpose, after each new minimization, the goodness-of-fit χ ν is computedand the new solution is rejected if its goodness of fit is larger than the one of the previoussolution. If this happens, while more than one SMR region had been flagged at the previousstep, the next in line is considered, and a new χ II ) minimization is attempted based onthe new set of initial parameters until the χ ν is successfully decreased or all the availableflagged regions have been tried. If the attempt to add a new spectral component fails for allSMR suggested locations, the χ II ) minimization stage is abnormally exited, and the partialsolution may be either rejected or flagged as possibly unreliable. As noted earlier, the FST instrument directly records time-domain data for 100 µ s (10 samples) at a time. These can be split into M accumulations of N -channel spectra, withthe constraint 2 M N = 10 , allowing us to trade spectral resolution, 500/ N MHz, for areduced level of statistical fluctuations, which depend on M . This motivates us to considerthe effect of this trade-off, using simulated data, on the algorithm’s ability to recover the truemodel. Figure 6 illustrates the algorithm’s performance in the case of a synthetic spectrum 24 –composed of the same underlying components as in Figure 5, but sampled with lower andhigher spectral resolutions. To mimic the design of an FFT-based spectrometer like FST,which digitizes the time domain signal at a fixed sampling rate, the changes in spectralresolution are accompanied by inversely proportional changes in the accumulation lengths M . The spectrum shown in Figure 6a,c, although less affected by noise fluctuations due toits increased accumulation length, M = 240, is negatively impacted by the proportionatelylower spectral resolution. It results in under-sampling the spectral shape, as quantified by thelarge χ ν = 10 .
45. In these conditions, the algorithm, while succeeding in reducing the SMRdeviations, fails to properly resolve the overlapping spectral peaks, and underestimates theamplitudes of the isolated ones, even though the final χ ν = 1 .
45 suggests a rather successfulfit. On the other hand, the higher frequency resolution of the spectrum shown in Figure 6e,g,which should facilitate the spectral separation of the overlapping peaks, is counterbalancedby the shorter accumulation length M = 24, which results in a higher level of statisticalfluctuations. This level of fluctuations prevents the algorithm from resolving three closelyoverlapping peaks. However, the local systematic SMR deviations revealed by Figure 6h,suggests the need for further model adjustment. To programmatically achieve this, thedetection threshold would have to be decreased at the cost of a higher false-alarm rate.Figure 6i presents the last stage of the fitting process of the same spectrum as in Figure 6g,but with the choice of a tiny increase in the PFA threshold of 10 − %. The algorithm nowsucceeds in lowering the local SMR deviations and identifying the correct number of spectralpeaks.
6. Analysis of the 06 December 2006 19:41:00 UT data recorded by FST
We now apply the fitting algorithm to a 60-s data segment recorded by the FSTinstrument (Liu et al. 2007) during the record-breaking solar flare of 06 December 2006(Cerruti et al. 2006). During this exceptional solar flare, which lasted more than an hour,the radio data recorded by the FST instrument in the 1 . − . > sfu) ever recorded during a solarradio burst. The spikes became less numerous during this late decay phase, but neverthelessoccurred in such high numbers that many overlapping spikes are seen, which the algorithmdescribed above must handle.FST observes the Sun interferometrically with three antennas. Although both total 25 –Fig. 6.— (a-d) First and last steps of the fitting algorithm for the same synthetic Gaussiansuperposition as in Figure 5, as would be seen with a larger number of accumulations ( M =240), which is achieved at the cost of decreasing the frequency resolution. The noise varianceis reduced by a factor of 2 .
5, but the lower frequency resolution does not fully resolve all ofthe spectral features in the model. (e-h) First and last steps of the fitting algorithm for thesame synthetic Gaussian superposition, as would be seen with higher frequency resolutionand a smaller number of accumulations ( M = 24), which increases the noise variance by afactor of 4. (i-j) The last step of the fitting algorithm corresponding to (g-h), but for a tinyincrease in PFA threshold of 10 − %, which now succeeds in separating the spike cluster nearfrequency bin 50 into its three components. 26 –power (signals from each antenna separately) and correlated data (cross-correlations of am-plitude and phase between each pair of antennas) are available, the FFT-based statistics wedescribe in this paper apply only to total power. As shown in Figure 7, the total powerrecords are essentially identical, which confirms that the fluctuations are dominated by thestatistics of the solar signal and not by uncorrelated thermal fluctuations in each antenna.Incidentally, it also places a constraint on the anisotropy in the spike emission on the scale ofthe spatial separation between the antennas (baselines of 136, 244, and 280 meters): there isclear one-to-one correspondence between spikes recorded by the three different antennas incontrast to what has been reported by D¸abrowski et al. (2011) for two different instrumentslocated in Toru´n and Ondˇrejov and so separated by 450 km.Figure 8a displays a 3-second detail of the M = 96 dynamic spectrum showing the spikeemission. The vertical dotted line indicates the time of a representative fit solution shown inFigure 8b. The χ ν goodness-of-fit parameter, and the average SMR deviation η are printedin the plot. Figure 8c displays the local SMR deviations. All but two local SMR deviationslie within the 0% PFA thresholds shown by the horizontal lines, but the algorithm did notadd spikes in these locations because doing so did not improve χ ν .We can check to what extent the model derived from the data is consistent with ourexpectation that the statistical noise is dependent solely on the amplitude of the signal. Todo this, we generate a synthetic dynamic spectrum by taking the superposition of Gaussiansfound for each spectrum of real data and adding numerically generated statistical noisecalculated for M = 96. The resulting synthetic spectrum for the 3-s period is shown inFigure 8d, which cannot be distinguished visually from Figure 8a. We then run the algorithmon these synthetic data and examine the fit results and the residuals. For example, the fittingsolution and corresponding SMR deviations for the same time as in Figure 8b,c are shownin Figure 8e,f, respectively. Although qualitatively similar to the real-data solution, thesynthetic-data solution in Figure 8e has only n = 12 spikes and a better χ ν = 1 .
29. TheSMR residuals in Figure 8c appear to be about twice those of Figure 8f. This quantitativemismatch between the data-based and simulation-based solution indicates that not all SMRresiduals in Figure 8c are due to M = 96 noise.This is further confirmed by Figure 9, which presents the distribution of χ ν for the150 times in Figure 8a,d and the distributions of the SMR deviations corresponding to eachspectral point of these times. The analysis is performed for three different accumulationlengths, M = 48 (left), M = 96 (middle), and M = 194 (right). Here we include only thosetimes for which the fitting algorithm succeeded to find a solution where all SMR deviationsare within their allowed range (0% PFA). Thus, for example, the fit in Figure 8b is rejected.The χ ν distributions for the synthetic spectra in Figure 9a-c cluster around unity, while 27 –those based on the original data are higher. This suggests that there is another source offluctuations than those based on accumulation number M , which becomes masked only when M is small enough that M -based fluctuations dominate. These fluctuations could be due tounresolved spikes. The distributions of the individual SMR deviations shown in the bottomrow of Figure 9 complement the above conclusion by revealing that the SMR distributions ofboth real (solid yellow) and synthetic (blue) data fits are not only closest to the theoreticaldistribution (red curve) in the M = 48, but also consistent to each other, i.e. they havesimilar excess values ǫ = 15% and ǫ = 7%, respectively, relative to the expected theoreticalSMR distribution given by equation (13).For the remainder of this section we continue to compare the results derived from actualdata in yellow with those derived from model-based synthetic spectra in blue, as a means ofhighlighting potential limitations. Figure 10 compares the data-derived distributions (solid yellow histograms) of spikenormalized amplitude (upper row), location frequency (middle row), and bandwidth (bottomrow), with the distributions derived from the corresponding synthetic spectra (overlaid bluehistograms).The amplitude normalization factor has been chosen based on the average FST back-ground in absence of any spike emission. For reference, the unit normalized amplitude, whichtherefore correspond to the average noise level, is indicated on each of the upper panel plotsby vertical red lines. The number of spikes extracted from real data, i.e. N = 3194( M = 48), N = 3426( M = 96), and N = 3386( M = 194), as well as the number and relative percent-age of the extracted spikes from the synthetic spectra, i.e. 78%( M = 48), 56%( M = 96),and 67%( M = 194), are shown in the plots. The amplitude distributions display high-endpower-laws, as well as a low-end rollover clearly associated with the background noise level.While the power-law slopes and intercepts of the real- and synthetic-data amplitude distri-butions are consistent in each case, it is evident that most of the missing spikes from thesynthetic-data distributions come from below the noise level. We conclude that the distri-butions above unit normalized amplitude are reliable, and further, as we concluded fromFigure 9, the M = 48 spectra provide the best combination of statistical noise and frequencyresolution. The spike extraction algorithm performance can be quantified by an overall 78%real-spike validation rate, as well as by only a minor difference between the lower ends ofthe real and synthetic data distributions. 28 –The center-frequency distributions of the spikes, shown in the middle row of Fig-ure 10, have a remarkably pure exponential distribution in the available spectral domain,1—1.5 GHz. The shape of the distribution, however, suggests that there are spikes at f > . f < − . GHz frequency range of the instrument, the pair of solid red horizontal lines represent a fixed1 − M Hz bandwidth range, and the dotted red horizontal line indicates the frequencyresolution corresponding to the accumulation length M . There is a clear separation of fittedspikes into two families, one within the physically justifiable confines of the inner red box,and another well outside the region of interest. Therefore, we conclude that the parameter-space boundaries drawn by the red solid lines, along with a minimum normalized-amplitudelimit of unity, can be used as filtering criteria to reliably separate the true solar spikes fromnon-solar artifacts.Consequently, Figure 12 presents the results of fitting the power-law indexes of thefiltered amplitude and bandwidth distributions. However, for illustration purposes, the dis-tributions shown in this figure were only partially filtered as follows: to build the densitydistributions shown on the upper row, only the bandwidth and location frequency filters wereapplied, while the bandwidth distributions shown on the bottom row were obtained fromdata filtered only by the location frequency and amplitude criteria. In contrast to Figure10, the distributions shown in Figure 12 are the density distributions, i.e. the counts ineach logarithmic bin have been divided by the corresponding variable bin width. On eachplot, the amplitude or bandwidth ranges used to perform the linear fits, which complete in 29 –each case the full filtering of data, are indicated by vertical red dotted lines. The resultingpower-law indexes, and their fit uncertainties, are indicated in corresponding colors for thereal (yellow) and synthetic (blue) distributions. On the bottom row plots, the vertical redsolid lines indicate the spectral resolution corresponding to each accumulation length M .The results displayed in Figure 12 confirm quantitatively that the best agreement be-tween the slopes of the real and synthetic distributions is reached in the case M = 48, themost favorable tradeoff between spectral resolution and the expected level of statistical noisefluctuations for the solar spikes observed in this event.As one final refinement, we noted earlier (Fig. 6) that the finding of spikes might beimproved by lowering the spike detection thresholds from the P F A = 0% level to a lessconservative, but still close to zero,
P F A = 10 − %. We therefore repeat the fitting based onthis detection threshold for M = 48 and summarize the results in Figure 13. In Figure 13awe compare the goodness of fit estimators of the original-data fits, h χ ν i = 1 . ± .
07, withthose from fitting the corresponding synthesized spectra, h χ ν i = 1 . ± .
05. Changing thePFA level to a small, but a nonzero value resulted in more than doubling the number of spikesextracted from data, i.e. N = 7538 vs N = 3194, while the percentage of self-consistentlyvalidated spikes increased slightly to 82% from the 79% figure indicated in Figure10a. Weinterpret this result as an indication that using a slightly less conservative PFA thresholdallows extraction of significantly more spikes, and consequently a better spike model, withouta significant increase in misidentification of statistical noise fluctuations as true solar spikes.To quantitatively support this assertion, we calculate that, when directly applied to thenumber of extracted spikes, N = 7538, the particular value we chose for the probability fora given detected peak to be created by a statistical fluctuation, P F A = 10 − %, indicatesthat the probability to have any one misidentified spike from the detected 7538 peaks is7 . × − ; hence, from a statistical perspective, less than one additionally detected spikecould have resulted from a false alarm flagging of the spectral model. We thus conclude that,for practical purposes, the criterion of keeping the absolute number of potentially false spikesbelow 1 may provide an objective means for adoption of a certain non-zero PFA suitable forthe amount of data being investigated.The distribution of the local SMR deviations shown in Figure 13b, when compared withFigure 9d, also indicates an improvement in algorithm’s performance, which is quantized bya reduction of the distribution excess parameters from ǫ = 15% to ǫ = 7 .
5% for the real datafits, and from ǫ = 7% to ǫ = 2 .
8% for the synthesized data fits.Similarly, the two-dimensional distributions of spike bandwidth versus spike locationshown in Figures 13c,d, confirm that the rectangular parameter space bounded by the limitsof the observed frequency range (vertical red solid lines), the spectral resolution (horizontal 30 –red dashed line) and the empirical high bandwidth limit of 100 MHz (horizontal red solidline), clearly confines the bulk of the true solar spike population and provides a reliablemeans for filtering out the non-solar artifacts.The very good visual agreement of the spike normalized-amplitude and bandwidth dis-tributions obtained from the original data (solid yellow) and synthesized data (blue) in Fig-ure 13e,f is quantitatively confirmed by their identical power law slopes, α = − . ± . β = − . ± .
04. The low-amplitude roll-over in Figure 13e is well below the noise (unitnormalized amplitude) limit, and hence the shape of the true distribution in this amplituderegion is not well recovered, as seen from the mismatch between the yellow and blue distribu-tions. In contrast, the low-bandwidth roll-over (peaking around 3 MHz) in Figure 13f is wellabove the frequency resolution (0.5 MHz) and is undoubtedly a real feature. The power-lawtail of this distribution results in a mean bandwidth (7.5 MHz) that is higher than the peak(most probable) bandwidth.
7. Discussion
The proposed spike decomposition algorithm allowed us to extract more than 7500individual spikes from a relatively sparse fragment of the powerful spike cluster observedfrom the 06 Dec 2006 flare; the obtained number of spikes is sufficient for a detailed statisticalanalysis of spike properties. The fitted distributions of the spike bandwidth have a (real)primary peak at a few MHz and two secondary peaks that we have concluded are artifacts—one below 1 MHz and the other around 1000 MHz. The first of them originates from verynarrowband (unresolved) fluctuations (possibly including interference) misidentified by thealgorithm as true spikes while the second is due slight background variations. Filtering thesetwo artifacts results in sharp edges in the distribution that prevent us from straightforwardlyfinding the moments of the true spike bandwidth distribution. The reliable portion of themeasured bandwidth distribution obeys a power-law with index around − .
25, mean of7.5 MHz, and mode of about 3 MHz.Similar asymmetric shapes of the bandwidth distribution are typical and have been re-ported for both decimeter and microwave spikes (Elgarøy & Sveen 1973; Csillaghy & Benz1993; Messmer & Benz 2000; Rozhansky et al. 2008; Nita et al. 2008). Rozhansky et al.(2008) employed a theory of the ECM bandwidth formation in a source with random in-homogeneities of the magnetic field developed by Fleishman (2004a) and demonstrated thatthe observed distribution was in almost perfect agreement with this theory. The comparisonof the moments of the observed and modeled distributions yielded an estimate of the corre-sponding magnetic irregularities (relatively small-scale turbulence) to be h δB i /B ∼ − . 31 –In our case, because of the aforementioned high-end and low-end artifacts, we cannot con-fidently estimate the higher moments of the bandwidth distribution, which is needed tocompute the magnetic irregularity level. We, nevertheless, can tell that a comparably lowlevel of the magnetic irregularities would be sufficient to yield the observed distribution.Although it may seem surprising that such a small value can in fact be measured through itseffect on the spike spectrum shape, one must keep in mind that these small magnetic fieldvariations have to be compared with another very small value—the ECM natural bandwidth,which is typically as small as ∼ .
2% (Fleishman 2004b).The central frequency, f , distribution carries some information of the global parameterrange at the source of the spike cluster, presumably, one (or more) coronal loops involved inthe flaring process. The available spectral range, 1–1.5 GHz implies the magnetic field rangeof 357–535 G if the ECM emission is produced at the fundamental of the gyrofrequency,or of 179–268 G in case of the second harmonic. Clearly, the actual range of the magneticfield is much broader because the central frequency distribution continues well outside ourspectral ’window’. The ability to produce either fundamental or harmonic ECM emissionalso constrains the plasma density, because the growth rates at these gyroharmonics havea strong dependence on the ‘plasma parameter’ Y = ω pe /ω Be (Fleishman & Melnikov 1998;Stupp 2000, and references therein). In particular, fundamental extraordinary-mode ECMemission requires Y < .
25 ( n e . × cm − ), fundamental ordinary-mode ECM emissionrequires 0 . < Y < n e . . × cm − ), and second harmonic extraordinary ECMemission requires 1 < Y < . n e . × cm − ). Identification of the emission mode(recall the spike emission is 100% polarized in this event) through the spatially resolvedmeasurements along with other context data can help to identify the emission mode and soimprove the use of the spike properties as a probe of source parameters.The derived distribution of spike amplitude does not display any apparent artifact andso is more conclusive. Indeed, in the meaningful range, roughly between 1 and 100, thedistribution is well fitted by a power-law with an index around − .
6, while it displays low-and high- amplitude rollovers outside this range. The low-amplitude rollover results frominability of the algorithm to fully recover the peaks having the amplitudes comparable orlower than the mean level of the signal. We emphasize that there must be plenty of the low-amplitude spike for a reason related to the FST recording mode. The FST time resolution is20 ms; however, the signal accumulation time is only 100 µ s within the 20 ms time interval.We checked that the spike signal level does not change measurably during this 100 µ s timeinterval implying that the typical spike duration is noticeably longer than 100 µ s. On theother hand, we did not find clear cases when the same spike would be recorded over twoconsequent 20 ms intervals, implying that the spike duration is shorter than 20 ms; notethat according a phenomenologically established regression law, the typical spike duration 32 –at 1 GHz is only marginally less than 20 ms (see Fig. 3 in Rozhansky et al. 2008). This meansthat during one measurement the given spike comes and goes, so the 100 µ s snapshot cancorrespond to any point of the rise, peak, or decay phase. Even if all the spikes were to havethe same peak amplitude, the described measurements would result in a broad distributionof the observed amplitude whose exact shape depends on the exact shape of the spike lightcurve; for example, for an exponential growth and a similar exponential decay this will bea power-law distribution with index −
1. The observed index, − .
6, deviates from − − . A. Moment-based Approximation for the PDF of the Goodness of FitEstimator
Although finding an exact closed form for the PDF of the random variable χ N maybe a difficult, if not impossible, mathematical task, the less challenging task of findinga approximation sufficient for practical applications may be straightforwardly accomplishedby evaluating the expectation E ( χ N ) = 1 and the higher statistical moments of the goodnessof fit parameter, i.e. E (cid:0) χ kN (cid:1) .For this purpose, we write χ N = 1 N N X j =1 ξ j , (34)where ξ ≡ √ M ( ρ j −
1) (35)is a random variable that, being linearly related to ρ j , follows a re-scaled and translated G PDF, i.e. G ∗ ( ξ ) = (cid:16) √ M (cid:17) M (cid:16) ξ + √ M (cid:17) M − Γ( M ) e −√ M ( ξ + √ M ) . (36)From this perspective, χ N represents the sample mean of the squared random variable ξ , which allows us to use the PDF of ξ for computing the expected variance of χ N from asimple formula that relates it to the expectations of ξ and ξ (Nita & Gary 2010a), σ χ N ≡ N (cid:2) E ( ξ ) − E ( ξ ) (cid:3) = 2 N (cid:18) M (cid:19) . (37)This shows that, for a fixed N , the variance of the goodness of fit estimator asymptoticallydecreases as the accumulation length M increases, from a maximum value of 8 /N toward2 /N , which is the variance of a standard χ N PDF.Having determined the mean and variance of the random variable equation (15) and,implicitly, of the goodness of fit estimator defined in equation (16), one may define anapproximation for the yet unknown PDF by choosing a functional form that matches thesefirst two standard moments.A natural choice for such approximation is equation (17), which in the limit of largeaccumulation length M reduces to a classic chi-squared distribution normalized by its degreesof freedom (equation 10). 34 –Within the limits of this approximation, the expected skewness, α and kurtosis, β , ofthe χ ν goodness of fit estimator are α ≈ √ √ N r M (38) β ≈ N (cid:18) M (cid:19) , and the probability to observe a given χ ν or larger is given by Eq. (18) B. An Unbiased Estimator for the Variance of an Accumulated FFT Spectrum If b σ j is an unbiased estimator for the variance of the G ( S j ; M, s j /M ) parent distributionof an accumulation of M raw FFT spectra, S j = P Mi =1 y j , then E [ b σ j ] = s j /M = M µ j , (39)where µ j and µ j represent the mean and, respectively, the variance of the parent distribution G ( S j ; 1 , µ j ) corresponding to a raw FFT spectrum.If b µ j represents the sample-based unbiased estimator of the variance of the G ( S j ; 1 , µ j )distribution, i.e. E [ b µ j ] = µ j . , immediately follows that b σ j = M b µ j . (40)Hence, using the generally valid expression for the sample-based unbiased variance estimator(Kendall & Stuart 1958), b µ j may be written as b µ j = M S ( j )2 − S j M ( M − , (41)where the S ( j )2 represents the accumulated square of the power, S ( j )2 = P Mi =1 y j , a quantitythat is not routinely available from a standard spectrum analyzer, but due to its theoreticallyproven practical benefits (Nita et al. 2007; Nita & Gary 2010a,b), it has recently become astandard output for a new generation of spectral instruments (Dou et al. 2009; Gary et al.2010; Deller et al. 2011).However, even in the absence of directly measured S ( j )2 , by using an intermediate resultof Nita & Gary (2010a), who showed that the statistical expectation of the accumulation S ( j )2 is related to the expectation of the accumulation S j , E [ S ( j )2 ] = 2[ E ( S j )] M + 1 , (42) 35 –one may express the unbiased variance estimator b µ j solely in terms of the known accumula-tion S j .Indeed, by taking the expectation of the unbiased variance estimator given by Equation41, one gets µ j = E h b µ j i = M E [ S ( j )2 ] − [ E ( S j )] M ( M −
1) = (43) E ( S j ) M ( M + 1) = E (cid:20) S j M ( M + 1) (cid:21) , which proves that b µ j = S j M ( M + 1) (44)is an unbiased sample–based estimator for the variance of the G ( S j ; 1 , µ j ) distribution.Consequently, from equation (40) it immediately follows that b σ j = S j ( M + 1) (45)is an unbiased sample–based estimator for the variance of the G ( S j ; M, s j /M ) distribution,which describes the statistical properties of an accumulation of M FFT raw spectra. C. A Comparison between the Least-Squares and Maximum LikelihoodSolutions
The vanishing conditions of the first order partial derivatives of the spectral log-likelihoodfunction with respect to its model parameters provide a system of ν equations with ν un-knowns p k , ∂λ ( p , p ..., p ν ) ∂p k ≡ M N X j =1 (cid:18) − S j b s j (cid:19) ∂ ln( b s j ) ∂p k = 0 . (46)The formal uncertainties of the solution parameters may be expressed in terms of the secondorder partial derivatives of the log-likelihood function (Bevington & Robinson 1992) as σ p k = 2 (cid:20) ∂ λ∂p k (cid:21) − = (47) 36 –1 M n N X j =1 h (cid:18) − S j b s j (cid:19) b s j ∂ b s j ∂p k − (cid:18) − S j b s j (cid:19) (cid:18) ∂ ln( b s j ) ∂p k (cid:19) io − . Since the system given by equation (46) can be solved analytically only in some simplecases (Bevington & Robinson 1992), the minimization of the spectral log-likelihood functiongenerally requires a numerical method (Barret & Vaughan 2012, and references therein).At this point, one may ask whether the functional form χ II ) defined in equation (20)may provide an alternative choice for a goodness-of-fit estimator. Indeed, using the proba-bility density function G (cid:0) ρ j , M, M (cid:1) , one may compute the expectation E (cid:2) χ II ) (cid:3) = ( M + 1)( M + 2)( M − M − N (48)and define the alternative statistical estimator χ II ) N ≡ ( M − M − M + 2) 1 N N X j =1 (cid:18) − b ρ j (cid:19) , (49)which has unity expectation. However, when compared with χ I ) N , its more cumbersomemathematical dependence on b ρ j makes χ II ) N a less suitable choice for a convenient goodness-of-fit estimator.To check the validity of our expectation that the least-squares minimization may lead toa point in the parameter space that is located in the vicinity of the true maximum likelihoodsolution, we compare in Figure 14 the results obtained by minimizing the negative log–likelihood function (equation 7) and the two alternative χ functions (equations 12 and 20)in the case of a simulated M = 12 FFT spectrum containing a single Gaussian spectral peak, s = ξ + αe − ( x − βγ ) , (50)where α , β , and γ , i.e. the amplitude, location, and dispersion, respectively, of the Gaussianpeak, and ξ represents a flat spectral background. The true signal, and the peaks estimatedby these three alterative methods are shown in panel (a).Based on Figure 14b we may conclude, by visual inspection, that the log–likelihoodminimization provides a more accurate representation of the hidden true signal than the χ I ) or χ II ) estimations, despite the fact that they correspond to smaller least–squareddeviations. The fact that the minimization of χ I ) results in an overestimation of the true 37 –amplitude α , while the minimization of χ II ) underestimates it, may be understood as adirect consequence of their reciprocal mathematical definitions, which result in differentweighting of the data points. At the same time, the minimized λ parameter and the non-minimized χ I ) and χ II ) values corresponding to the negative log–likelihood estimation arecloser to the true values produced by the random realization used for this test.Figure 14b displays the local SMR deviations corresponding to the noise–contaminatedtrue signal and the three alternative estimations, as well as the average SMR deviations, b η = (1 /N ) P Nj =1 b ρ j . While the lower and, respectively, higher than unity averaged SMRvalues resulting from the two χ minimizations are related to the same weighting bias dis-cussed above, the perfect b η = 1 .
00 from the log–likelihood minimization is consistent withthe known fact that the log-likelihood approach may not necessarily provide the most accu-rate description of an individual random realization, but does provide the statistically mostfavorable set of parameters that could have produced the observed samples.To provide the means for a more quantitative assessment of this test, we display in theother four panels of the Figures 14 the shapes of the minimized functions in the vicinity oftheir minima obtained by keeping three out of the four parameters fixed at their estimatedvalues, while the remaining one has been slightly varied. To facilitate a direct comparison,the χ curves have been normalized by their corresponding degrees of freedom, while theircorresponding minimum values have been subtracted from the λ curves. The true andestimated values of the spectral peak parameters are indicated in solid colors, while the3 σ standard ranges of uncertainty, computed according to the equations (47) and (21), areindicated by color-coded dashed lines. D. Standard Moments of the Mean of an SMR Compact Region
One may compute the characteristic functions corresponding to the probability dis-tribution functions describing the means η a and η b by taking the n th power of the Fouriertransforms of the conditional PDFs defined by equation (26), followed by the standard changeof variable t → t/n in their arguments (Kendall & Stuart 1958). This procedure leads toΦ( t ) (cid:12)(cid:12)(cid:12) { ρ j > } = (cid:20) M M E (1 − M, M − itn ) γ ( M, M ) (cid:21) n (51)Φ( t ) (cid:12)(cid:12)(cid:12) { ρ j < } = (cid:20) M M ( M − itn ) − M [Γ( M ) − γ ( M, M − itn )] γ ( M, M ) (cid:21) n , where E ( a, z ) = Z ∞ e − zt t − a dt, (52) 38 –is the standard exponential integral function.Although closed forms for the inverse Fourier transforms of equations (51) have notbeen found, which would have directly provided the PDF we are interested in, by following astandard procedure (Kendall & Stuart 1958), one may use these characteristic functions tocompute the infinite sets of statistical moments associated with the random variables underinvestigation.Using the characteristic functions given by equation (51), one may straightforwardlycompute the first four central moments of the random variables η a and η b (equation 25),which are given by equations (53) and (54), respectively. µ = Γ( M + 1 , M ) M γ ( M, M ) (53) µ = γ ( M, M )Γ( M + 2 , M ) − Γ( M + 1 , M ) nM γ ( M, M ) µ = 1 n M γ ( M, M ) h M + 1 , M ) − γ ( M, M )Γ( M + 1 , M )Γ( M + 2 , M )+ γ ( M, M ) Γ(3 +
M, M ) i µ = 1 n M γ ( M, M ) h n − M + 1 , M ) − n − γ ( M, M )Γ( M + 1 , M ) Γ( M + 2 , M ) − γ ( M, M ) Γ( M + 1 , M )Γ( M + 3 , M )+3( n − γ ( M, M ) Γ( M + 2 , M ) + γ ( M, M ) Γ( M + 4 , M ) i µ = 1 − e − M M M − Γ( M ) − γ ( M, M ) (54) µ = e − M nM [Γ( M ) − γ ( M, M )] n e M M Γ( M ) + e M M M γ ( M, M ) + e M M γ ( M, M ) − e M Γ( M )[ M M + 2 e M M γ ( M, M )] − M M o µ = e − M n M [Γ( M ] − γ ( M, M )] × n e M M Γ( M ) + 3 e M M M γ ( M, M ) 39 –+ e M ( M − M M γ ( M, M ) − e M M γ ( M, M ) − M M − e M Γ( M ) h e M M γ ( M, M ) − ( M − M M i + e M Γ( M ) h e M M γ ( M, M ) − e M ( M − M M γ ( M, M ) − M M io µ = e − M n M [Γ( M ] − γ ( M, M )] × n n − M M + 3 e M M ( M n + 2)Γ( M ) − n − e M M M γ ( M, M ) − e M M M [ M (6 n − − n + 11] γ ( M, M ) + e M M M [(6 n − M + 6] γ ( M, M ) +3 e M M ( M n + 2) γ ( M, M ) − e M Γ( M ) h M M [(6 n − M + 6]+12 e M ( M n + 2) γ ( M, M ) i + e M Γ( M ) h M M [(6 n − M + 3 n − e M M M [(6 n − M + 6] γ ( M, M )+18 e M M ( M n + 2) γ ( M, M ) i − e M Γ( M ) h − n ) M M − e M M M × [(6 n − M − n + 11) γ ( M, M )3 e M M M [(6 n − M + 6]Γ[ M, M ] e M M ( M n + 2) γ ( M, M ) io E. Pearson Type I Approximations vs. Monte Carlo Simulations
One should always take in consideration that the Pearson Type I are only approxima-tions to the true distributions, therefore their performance and limitations should be alwaystested on a case-by-case basis.To facilitate such critical assessment, we present in Table 2 the inferred Pearson Type Iparameters for the particular case M = 48 used in the Monte Carlo experiment described in 40 – § α , the allowed range { a, b } ,and the 0 . t a an t b , as defined in § tunneling probabilities for η b and η a to begreater and, respectively, lesser the unity are given. The fact that these probabilities are non-zero is the direct consequence of the fact that both η b and η a Pearson I approximations haveranges that actually cross above and respectively below unity, although it is mathematicallyimpossible to observe in practice means above or below unity while averaging quantities thatare all below and respectively above these borders.However, except for the special case n = 1, these pseudo tunneling probabilities arebelow unity and their negligibility increases as fast as one order of magnitude as n increasesby 1, despite the fact that neither the distribution limits nor the relative tunneling depthschanges noticeably with n as M is kept fixed. These results suggests that the accuracy ofthese approximations may increase with n .The special case n = 1 may be easily explained by recalling the fact that for n = 1 wealready found the true distributions, which are the truncated Gamma distributions givenby equation (26), which we have used as basis for our derivation. The n = 1 distributions,despite having the firsts standard moments satisfying the Pearson Type I region criterion,as clearly illustrated in Figure 4, do not not satisfy the basic Pearson family criterion ofhaving a continuous derivative. Therefore, the available true distributions must be used forthe case n = 1 and, since in all other cases the true distributions are not known, we comparethem with the outcome of the Monte Carlo simulations in order to validate the PearsonType I approximations derived for n >
1. The result of this comparison, which we presentin this Appendix, allows us to conclude that, at least from a practical point of view, the truetail probabilities of the SMR distributions are accurately estimated by the Pearson Type IApproximations in all cases including the special case n = 1.Figure 15a shows the randomly generated distributions of the mean SMR s b and η a for n = 1 ,
6, while Figure 15b shows the distributions for n = 7 ,
12. For comparison,their corresponding Pearson Type I PDF approximations, characterized by the parametersshown in Table 2 are drawn with solid lines. A very good agreement between simulationsand theoretical approximations is evident except for the special case n = 1, for which thecorresponding true theoretical distributions are shown as dashed-dotted lines. Although,for reasons already discussed in the precedent section, the differences between the n = 1approximations and the true distributions are significant especially in the vicinity of theirpeaks, Figure 15a suggests that, even in the case n = 1, the Pearson Type I curves maystill offer good approximations for the true tail probabilities. To support this assertion, wepresent in Table 1 a comparison between some key parameters of the approximative and 41 –exact distributions for M = 48 and n = 1. Table 1 reveals that the standard 0 . t a and t b , computed using the Beta approximations are practically identicalwith the thresholds computed using the exact distributions. Moreover, the true probabilitiesof observing SMR regions beyond the limited ranges of the
Beta approximations are alsopractically negligible.From a practical point of view, to address the fitting problem, we are interested onlyin the tail probabilities. Therefore, we may assert that even for the case n = 1 it would besafe to assume that any outlier located outside the range of the Pearson Type I distributionswould be due to the model function rather than the result of some expected statisticalfluctuations. This result may greatly simplify and speed up the fitting algorithm since,instead of paying the computational cost of having to evaluate an incomplete Beta functionfor checking the statistical significance of each observed systematic deviation above or belowunity, and compare it with an arbitrary chosen non zero probability threshold, one may choseinstead to compute only the limited ranges of the Pearson Type I approximations, such asthose listed by Table 2 , and use them to reject any systematic region that is situated beyondthese limits.To conclude this section, we present in Figure 16 two data sequences corresponding tothe simulated data regions, that contain the largest SMR deviations below (0.554–panel a)and above (0.897–panel b)unity, the largest systematic SMR regions below ( n = 20–panelc) and above ( n = 23–panel d) unity, as well as the least probable SMR regions randomlygenerated below (2 . × − %–panel e) and above (2 . × − %–panel f) unity.The ranges of the selected regions are indicated by blue vertical lines and their meandeviations are shown as horizontal blue solid lines on each panel. The red dashed horizon-tal lines indicate the corresponding standard 0 . n , the size probability, p ( n ), and the SMRmean probability p ( s ) to be observed amongst its own size class are also indicated.Note that, in a real situation, both n = 1 SMR regions shown in panels ( a ) and b wouldfail the standard 0 . a region, the panel b regionwould also survive the 0% PFA test, being still inside the limited boundaries of its PearsonType I approximation, despite its comparatively lower probability to be observed amongstits own class. At the same time, despite being very unlikely to be randomly generated, asresulted from their size distributions (Eqn. 24), i.e. 9 . × − % and 2 . × − %, respec-tively, the regions shown in panels ( c ) and ( d ) would survive both standard and boundarytests due to their high probabilities to be observed amongst their own class, i.e. 51 .
9% and67 . T o t a l P o w e r [ a r b . un i t s ] Ant 5Ant 6Ant 7
Fig. 7.— Comparison of total power records for each of the three antennas (5, 6 and 7) for asingle time as a function of frequency, after normalizing to the background continuum. Thesignal for antenna 6 is shifted by 2, and the signal for antenna 7 is shifted by 4 in amplitudefor clarity. The signals are essentially identical down to even minor fluctuations, especiallyat the lower end of the frequency range where the continuum is stronger.Table 1. Exact PDFs vs. Pearson Type I Approximations for M = 48 and n = 1Type p ( s < a )% a t b (0 . p ( s > · · · p ( s < a )% a t a (0 . p ( s > · · · − . N = 511 frequency channels( ∼ M = 96 accumulation length. The abundant spike emissiondisplaying different degrees of overlapping is evident. The dotted vertical line indicates aselected time frame used for illustration. b) Data and estimated spectral components for theselected time frame. The estimated number of spikes, n = 13, the goodness of fit χ ν = 2 . η = 1 .
01 are indicated in the figure inset. c) Local SMRdeviations corresponding to the solution shown in panel (b). All but two SMR deviationslie within the maximum allowed range (0% PFA) of the n = 1 Pearson Type I PDF. Panels( d, e, f ): The same as in panels ( a, b, c ), but for the synthetic data set built by contaminatingthe solution estimated from real data with pure M = 96 statistical noise. 44 –Fig. 9.— Goodness-of-fit analysis for the 150 spectra in Figure 8a,d, obtained with threedifferent accumulation lengths, M = 48 ,
96 and 194. Upper panels: χ ν distributions for theoriginal data (solid yellow) and for synthetic data (blue lines). The means and standarddeviations of the χ ν parameters are indicated on each plot. Lower panels: Distributionsof the SMR deviations for the original (solid yellow) and synthetic (blue lines) data. Thedistributions theoretically expected according to equation (13) are plotted in each panel (redlines) and the percentage ǫ of points falling outside the theoretical distribution is shown oneach plot. Only those time frames for which the algorithm has found solutions with allSMR deviations ranging within the expected theoretical limits, i.e. the 0% PFA thresholdsindicated by vertical red lines, were selected for this analysis. 45 –Fig. 10.— Upper row: Normalized amplitude distributions. Middle row: Central frequencydistributions. Bottom row: Bandwidth distributions. The first, second, and third rowdisplay the distributions corresponding to M = 46, M = 96 and M = 194, respectively.In all panels, the solid yellow histograms correspond to the spikes extracted from real data,while the blue histograms correspond to the spikes extracted from synthetic data generatedbased on the previously extracted spikes. For reference, we indicate in the upper row plotsthe normalized unity amplitude by vertical red solid lines. In the middle row plots, the pairsof vertical red solid lines mark the physical boundaries of the observed 1 − . − . GHz frequency range of the instrument, the pair of solidred horizontal lines represent a fixed 1 − M Hz bandwidth range, and the dotted redhorizontal line indicates the frequency resolution corresponding to the accumulation length M . 47 –Fig. 12.— Power-law fitting of amplitude and bandwidth distributions derived from real(yellow) and synthetic (blue) data. The amplitude distributions displayed on the upper rowwere filtered based on the bandwidth and location frequency boundaries shown on Figure11, while the bandwidth distributions shown on the bottom row were filtered based on thelocation frequency criterion, as well as by the minimum unit relative amplitude criterion.The fitting ranges, which complete the 3–criterion filtering process are indicated by verticalred dash-dotted lines and the power-law indexes and their associated fit uncertainties ofreal and synthetic distributions are indicated in corresponding colors on each plot. On thebottom row plots, the vertical red solid lines indicate the spectral resolution correspondingto each accumulation length M . 48 –Fig. 13.— Characteristics of the extracted spikes from the M = 48 spectrum with p = 10 − %PFA as obtained from the original (yellow) and corresponding synthetic (blue) data. Upperrow: Goodness of fit analysis presented in panels (a) and (b) in the same format as in Figure9. Middle row: Two-dimensional distributions of spike locations versus spike bandwidth, asin Figure 11. Panels (c) and (d) present the same information with alternate order in whichspikes extracted from the original (yellow) and synthetic (blue) data are plotted on top ofeach other. Bottom row: The amplitude (panel e) and bandwidth (panel f) distributions,and their corresponding power-law fits, as in Figure 12. 49 –Fig. 14.— Fit results in the case of a simulated M = 12 FFT spectrum containing a singleGaussian spectral peak superimposed on a flat background. The color-coded lines represent:red–true signal; blue– χ minimization according to Equation 12; green– χ minimization ac-cording to Equation 20; yellow– λ minimization according to Equation 7. The solutions, andtheir corresponding SMR values are shown on panels (a) and (b) respectively. The associ-ated χ ν (panel a), λ (panel a), and η (panel b) parameters are displayed in correspondingcolors. The shapes of the minimized functions in the vicinity of their minima obtained bykeeping three out of the four parameters fixed at their estimated values, while the remainingone has been slightly varied are shown in panels c–f. The χ curves have been normalizedby their corresponding degrees of freedom, while their corresponding minimum values havebeen subtracted from the λ curves. The true and estimated values of the spectral peak pa-rameters are indicated in solid colors, while the 3 σ standard ranges of uncertainty, computedaccording to the equations (47) and (21), are indicated by color-coded dashed lines. 50 –Fig. 15.— The observed simulated distributions of compact SMR regions of different sizes n and fixed accumulation length M . The good agreement with the Pearson type I approxima-tions is evident,(note the logarithmic scale on panel ( b )), for all cases except n = 1 which isexactly described by the truncated Gamma distributions (Eqn. 26) indicated on panel a bythe dashed-dotted lines. 51 –Fig. 16.— Simulation highlights: a) the region containing the largest SMR deviation belowunity; b) the region containing the largest SMR deviation above unity; c) the largest regiondisplaying a systematic deviation below unity; d) the largest region displayed a systematicdeviation above unity; e) the least probable region below unity; f) the least probable regionabove unity. The ranges of the selected regions are indicated by vertical blue solid lines andtheir corresponding mean deviations are shown by horizontal blue solid lines on each panel.The red dashed horizontal lines indicate the standard 0 . n , the size probability, p ( n ), and the SMR mean probability p ( s ) to be observed amongst itsown size class are also indicated. 52 –Hence, combining the probabilities provided by Eqn. 24 and Eqn. 29, one may computea relative ranking in the increasing order of their absolute probabilities to being randomlygenerated, i.e. ( d ) : 1 . × − % , ( f ) : 2 . × − % , ( e ) : 2 . × − % , ( a ) :2 . × − %, and ( c ) : 5 . × − %. This reveals that the largest ( n = 23) SMRregion shown in panel ( d ) is indeed the least likely region to be randomly generated, whilethe n = 20 region shown in panel c is actually the one that is most likely to be randomlygenerated, despite being much larger than the n = 2 regions.We consider as a significant outcome of our simulation the fact that out of a total of498 ,
611 SMR regions having n ≥
2, none has been found to cross the limited boundaries ofits corresponding Pearson Type I PDF Approximation. Moreover, out of a total of 500 , n = 1.This work was supported in part by NSF grants AGS-1250374 and NASA grants NNX11AB49Gand NNX13AE41G to New Jersey Institute of Technology. This work also benefited fromworkshop support from the International Space Science Institute (ISSI). 53 –Table 2. Parameters of the Pearson Type I Approximations for M = 48 and n = 1 , α a t b b p ( s > p ( s <
1) a t a b α @0 .
13% % % @0 . REFERENCES
Altyntsev, A. T., Lesovoi, S. V., Meshalkina, N. S., Sych, R. A., & Yan, Y. 2003, A&A, 400,337Aschwanden, M. J. 1990, A&AS, 85, 1141Aschwanden, M. J., Dennis, B. R., & Benz, A. O. 1998, ApJ, 497, 972Barret, D., & Vaughan, S. 2012, ApJ, 746, 131Benz, A. O. 1985, Sol. Phys., 96, 357—. 1986, Sol. Phys., 104, 99Benz, A. O., Grigis, P. C., Csillaghy, A., & Saint-Hilaire, P. 2005a, Sol. Phys., 226, 121—. 2005b, Sol. Phys., 226, 121Benz, A. O., Monstein, C., Beverland, M., Meyer, H., & Stuber, B. 2009, Sol. Phys., 260,375Benz, A. O., Saint-Hilaire, P., & Vilmer, N. 2002, A&A, 383, 678Bevington, P. R., & Robinson, D. K. 1992, Data reduction and error analysis for the physicalsciences; 2nd ed. (New York, NY: McGraw-Hill)Bowman, K., & Shenton, L. 2007, Far East Journal of Theoretical Statistics, 23, 133Cerruti, A. P., Kintner, P. M., Gary, D. E., et al. 2006, Space Weather, 4, 10006Chernov, G. P., Sych, R. A., Yan, Y., et al. 2006, Sol. Phys., 237, 397Csillaghy, A., & Benz, A. O. 1993, A&A, 274, 487D¸abrowski, B. P., Rudawy, P., Falewicz, R., Siarkowski, M., & Kus, A. J. 2005, A&A, 434,1139D¸abrowski, B. P., Rudawy, P., & Karlick´y, M. 2011, Sol. Phys., 273, 377Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP, 123, 275Dou, Y., Gary, D. E., Liu, Z., et al. 2009, Publications of the Astronomical Society of thePacific, 121, 512Dr¨oge, F. 1977, A&A, 57, 285 56 –Elgarøy, Ø., & Sveen, O. P. 1973, Sol. Phys., 32, 231Fleishman, G. D. 2004a, ApJ, 601, 559—. 2004b, Astronomy Letters, 30, 603Fleishman, G. D., Gary, D. E., & Nita, G. M. 2003, ApJ, 593, 571Fleishman, G. D., & Melnikov, V. F. 1998, Uspekhi Fizicheskikh Nauk, 41, 1157Gary, D. E., Hurford, G. J., & Flees, D. J. 1991, ApJ, 369, 255Gary, D. E., Liu, Z., & Nita, G. M. 2010, PASP, 122, 560G¨udel, M., & Benz, A. O. 1988, A&AS, 75, 243G¨udel, M., & Zlobec, P. 1991, A&A, 245, 299Heydt, G., Fjeld, P. S., Liu, C., et al. 1999, Power Delivery, IEEE Transactions on, 14, 1411Holman, G. D., Eichler, D., & Kundu, M. R. 1980, in IAU Symposium, Vol. 86, RadioPhysics of the Sun, ed. M. R. Kundu & T. E. Gergely, 457–459Isliker, H., & Benz, A. O. 1994, A&AS, 104, 145Kendall, M. G., & Stuart, A. 1958, The Advanced Theory of Statistics, Vol. 1, London:Charles Griffin and CoLiu, Z., Gary, D. E., Nita, G. M., White, S. M., & Hurford, G. J. 2007, PASP, 119, 303Magdaleni´c, J., Vrˇsnak, B., Zlobec, P., Hillaris, A., & Messerotti, M. 2006, ApJ, 642, L77Mannan, M. A., Kassim, A. A., & Jing, M. 2000, in Application of image and sound analysistechniques to monitor the condition of cutting tools., 969–979Melrose, D. B., & Dulk, G. A. 1982, ApJ, 259, 844Meshalkina, N. S., Altyntsev, A. T., Sych, R. A., Chernov, G. P., & Yihua, Y. 2004,Sol. Phys., 221, 85Messmer, P., & Benz, A. O. 2000, A&A, 354, 287M´esz´arosov´a, H., Karlick´y, M., Veronig, A., Zlobec, P., & Messerotti, M. 2000, A&A, 360,1126Nita, G. M., Fleishman, G. D., & Gary, D. E. 2008, ApJ, 689, 545 57 –Nita, G. M., & Gary, D. E. 2010a, PASP, 122, 595—. 2010b, MNRAS, 406, L60Nita, G. M., Gary, D. E., Liu, Z., Hurford, G. J., & White, S. M. 2007, PASP, 119, 805Pearson, K. 1895, Royal Society of London Philosophical Transactions Series A, 186, 343Rozhansky, I. V., Fleishman, G. D., & Huang, G.-L. 2008, ApJ, 681, 1688Slottje, C. 1978, Nature, 275, 520—. 1981, Atlas of Fine Structures of Dynamic Spectra of Solar Type IV-dm and Some TypeII Radio Bursts (The Netherlands: Dwingeloo, 1981)St¨ahli, M., & Magun, A. 1986, Sol. Phys., 104, 117Stepanov, A. V. 1978, Pisma v Astronomicheskii Zhurnal, 4, 193Stupp, A. 2000, MNRAS, 311, 251Treumann, R. A. 2006, A&A Rev., 13, 229