Some problems in multifractal spectrum computation using a statistical method

In practical calculations of a multifractal spectrum only limited data are available due to data overflow, measurement errors and limited calculation time. An effective method to reduce the data overflow is proposed. Some parameters are introduced to evaluate the incomplete degree of a partial multifractal spectrum. Quantitative expressions of the evaluation parameters on the partial multifractal spectra calculated using a statistical method for the Cantor sets p/1−2p/p and p/0/1−p are derived. Approximate evaluation parameters of the partial multifractal spectra calculated using the statistical method are estimated for two examples with a random fractal character. The characteristic of the evaluation parameters for the above examples is discussed in detail.


Introduction
Since the concept of multifractals was introduced in 1970s by Mandelbrot, it has been used to describe and distinguish varieties of complicated figures, systems and processes in nature. Now, it has been applied to more and more fields [1] and the theory has been developed quickly.
Many algorithms have been developed to measure the multifractal dimensions [2]- [13]. Using these standard algorithms, we can derive the Hentschelt-Procaccia fractal dimension D(q) [2], the f(α) spectrum and the other related fractal measures for precise data.
These algorithms have been tested by introducing some 'noises' (or 'errors') on the data of a system with a regular fractal characteristic. The advantage of different algorithms has been discussed intensively, such as consequences of noise on the effectiveness of fractal analysis algorithms [3,4], the convergence of the box-counting fractal [3], the comparison of the numerical correlation integral algorithm with the Badii-Politi algorithm [13] by applying the algorithms to Euclidean point sets, Koch asymmetric snowflake and the imprecise data produced by perturbing the x and y coordinates of the fractal points by 'random values' and so on.
Because of the restriction on practical measurements, we can obtain an incomplete multifractal spectrum in the finite scaling range only. The fractal characteristic in a finite scaling range has been studied in detail. However, to the best of our knowledge, quantitative discussion on the incomplete characteristic of a multifractal spectrum has not been reported in the literature so far.
A method to calculate multifractal spectra was proposed by Halsey in 1986 [14]. By means of it, the multifractal spectra for quantities or states of a system with a random fractal character can be calculated from measured data. It is so easy and simple that it has been widely used to compute varieties of multifractal spectra. It is called the statistical method in the present paper for the sake of simpleness.
The wavelet transform modulus maxima (WTMM) method [15] is a well-known method to investigate the multifractal scaling properties of fractal and self-affine objects in the presence of nonstationarities.
The detrended fluctuation analysis (DFA) method is a widely used technique for the determination of (mono) fractal scaling properties and the detection of long-range correlations in noisy, nonstationary time series. The multifractal-DFA (MF-DFA) method [16], which is based on a generalization of the DFA, is produced for the calculated multifractal characterization of nonstationary time series.
Theoretically, it can be demonstrated that since the range of a special parameter in the practical computation cannot be taken as infinite, a complete multifractal characteristic cannot be gotten using any one among the three methods (the statistical, the WTMM and the MF-DFA methods).
In the statistical method, the special parameter q is used to calculate the function χ q (l) [14]; in the WTMM method the special parameter is used to calculate the function Z(q, s) [16]; and in the MF-DFA method the special parameter is used to calculate the function F q (s) [16], where l and s are scale parameters.
The aim of the present paper is an attempt to make the first step to solve the quantitative description of an incomplete multifractal spectrum. Obviously, the completeness of a calculated multifractal spectrum is important. The more complete the multifractal spectrum, the more actual the description. By using these more complete spectra, we can find out the rule among complicated systems more easily.

A method of computing multifractal spectra
The calculation steps of the statistical method are briefly described as follows [14,17]. It is well known that the introduced multifractal is to study the singular distribution properties of some quantities or states of a system. At first, a suitable space dimension is selected, which is 1, 2 or 3 according to the fractal character of the actual quantity. Then, the space is divided into many normalized boxes of size l (l 1 and the first one is l = 1). p i (l) is the distribution probability of a quantity or a state of a system studied in box i. Thus, the multifractal is described as where α is the singularity of the probability subset, N α (l) is the number of boxes of size l with the same probability and f(α) is the fractal dimension of subset α. The multifractal spectrum f(α) can be calculated by the following function: where q is the moment order, −∞ q ∞. If the object studied has a characteristic of multifractal, then we have where The χ q (l) with the same q but different l can be calculated and the curve of ln χ q (l) − ln l is plotted. For a fractal set, the curve is a straight line. Therefore, τ(q) can be obtained from the slope of the line In practical calculation, the range of l corresponding to the linear part of the curve is called as the scaling range of the multifractal. Then, the α subset can be obtained from equation (5) by solving it for different q. The expressions for α and f are α = dτ dq (7) and It is known from expression (4) that only one p in the sum plays the main role and its corresponding τ(q) is the smallest among all the p i . Consequently, we get the following results: • The subset when q = 0 corresponds to the one with the maximum fractal dimension f max , and the related q and α are denoted as q 0 and α 0 , respectively. • The subset as q < 0 corresponds to that with α > α 0 , and the smaller the q, the larger the α. • The subset for q > 0 corresponds to that with α < α 0 and the larger the q, the smaller the α.

Finite scaling range
In practical measurements, we cannot get sufficient data that distribute in many orders of magnitude of the scaling range. It is obvious that a phenomenon occurs only during a finite time and in a finite space, and the accuracy of data is also limited by the resolutions of the instruments used. Considering that real structures exhibit self-similarity only over a finite range of scales, Berntson and Stoll [18] proposed a technique by which they estimated only over a statistically identified region of self-similarity and introduced the finite scale-corrected dimension (FSCD). Avnir et al [19] surveyed a total of 96 papers published in all Physical Review journals (Phys. Rev. A to E and Phys. Rev. Lett.) from 1990 to 1996, in which the experimental data were reported to exhibit a fractal character. The results indicated that the scaling range of the experimentally declared fractality was extremely limited, which was centred in 1.3 orders of magnitude and spanned mainly between 0.5 and 2.0. As the distribution probability decreases due to the limit of experimental conditions, some unusual phenomena may appear in the multifractal spectra. Wang and Wu [20] analysed the problem in detail.

Incomplete multifractal spectra
According to the fractal theory, q should range from −∞ to +∞ to get a complete multifractal spectrum. However, this condition cannot be satisfied.

Data overflow
In fact, the data overflow restricts the range of q in the calculation. On a computer, the range of a double-precision number is from 10 −308 to 10 308 . Therefore, p q i in expression (4) may result in a data overflow, as q is 100. Now, a method is introduced to improve this situation greatly.
The key point is to indicate each positive real number x in a special structure described by two parameters as follows. Let x = a × 10 b , where b is an integer and 1 a < 10. A positive x can be transferred to the new structure as follows: • If x 1, b = int(log 10(x)) and a = x/pow (10, b).
The new structure can be obtained for different operations easily.
The standardization process is as follows. Let It should be noticed that the a q 1 may still be overflowed as |q| is large enough. The following method is helpful to reduce the overflow probability of the p q i in expression (4).
where A is an adjustable positive number and should be satisfied with the condition, 10 −308 < p A i < 10 308 , j is a natural number and q = q − j × A. Standardize p A i and p q i . The calculation range of q increases greatly by using j times of multiplication instead of one time of standardization p q i except that the p i is overflowed. However, it leads to the decreasing computation speed. Although the calculation range of q can be taken wider and wider, it is still finite. Therefore, the evaluation of the partial multifractal spectra is very important.

Effect of measurement errors
Let σ τ indicate the standard error of τ(q), which arises from the fitting line of the ln χ q (l) ∼ ln l curve. If q is very small, α = dτ/dq ≈ τ/ q. So, the standard error of α is approximate to that of τ. Let the errors of two near τ be approximately equivalent. On the basis of equations (7) and (8), we have Obviously, σ f increases rapidly with the increasing |q|. Therefore, the range of q should be restricted to get exact results. Obviously, only the data whose errors satisfy our requirements are useful. Consequently, the errors from measurement bring us the new restriction on the calculated multifractal spectra.

Evaluation parameters
For multifractal spectra with the shape of a bell or a hook, some parameters are defined as follows.
The α corresponding to the maximum fractal dimension f max is denoted as α 0 . The part with α < α 0 is denoted as section I and that with α > α 0 , section II. The f(α) corresponding to α min and α max is represented as f 01 and f 02 , respectively. Let For a complete multifractal spectrum, f 01 , f 02 , α 01 and α 02 are used to denote the changes of f(α) and α in sections I and II, respectively and α 0 is the total change of α.
The calculated minimum and maximum α are referred to as α 1 and α 2 , respectively; obviously, α 1 < α 0 < α 2 . The symbols, f 1 and f 2 , indicate the f(α) corresponding to α 1 and α 2 , respectively. Let For a partial multifractal spectrum, f 1 , f 2 , α 1 and α 2 are used to indicate the changes of f(α) and α in sections I and II, respectively, and α is referred to as the total change of α.
The parameters α 1 / α 01 , α 2 / α 02 , α/ α 0 , f 1 / f 01 and f 2 / f 02 are called as the evaluation parameters. A partial multifractal spectrum can be evaluated using these parameters. We shall discuss the evaluation parameters of the partial multifractal spectra for two Cantor sets and two examples of the random fractals calculated by the statistical method in detail.

Multifractal spectrum computed from generator
For this Cantor set, the original region is divided into three segments with the same size in each step, but the mass distribution probability from the left to the right segment is p, 1 − 2p and p, respectively. Hence, the Cantor set is expressed as the Cantor set p/1 − 2p/p. The method of computing the multifractal spectrum for different Cantor sets has been described in the literature [14,17,20]. We will briefly summarize as follows.
After k steps, the total segments with a size of l = ( 1 3 ) k are 3 k . The mass probabilities form a distribution set {p i }, The number of intervals with p i is Let m = kξ where 0 ξ 1. As k → ∞ (or l → 0), ξ will be continuous. Therefore, Taking the logarithm, we have Since According to the Sterling formula n! ≈ √ 2πn(n/e) n , we have From expression (15), it is known that f takes its maximum f max = 1 as ξ = 2 3 .
It can be seen from equation (13) that If ξ 01 , ξ 02 , ξ 1 and ξ 2 are known, then f 01 , f 02 , f 1 and f 2 can be found by expression (15), and q 1 and q 2 are given by expression (21). As f max = 1, then we have On the basis of expression (21), the following expressions are obtained: and 6. Cantor set p/0/1 − p For this Cantor set, the calculation process and the results of the multifractal spectrum have been given in detail in the literature. Hence, we shall just give the expressions for α and f(α) deduced using these two methods directly. Then, we shall derive the expressions for the evaluation parameters using our discussion mode. For this Cantor set, the original region is divided into three equivalent segments with p, 0 and 1 − p of the mass distribution probability from the left to the right segment, respectively, in each step. Similar to the Canter set p/1 − 2p/p, the following expressions are obtained from the generator: where 0 ξ 1. For the statistical method, we have and where −∞ q ∞. Let Thus, expression (27) is same as (29) and expression (28) is equal to expression (30).

The effect of q on the completeness of calculated results
For the two Cantor sets, the quantitative expressions for the evaluating parameters are derived. The expressions show that the range of q in calculating the multifractal spectra greatly affects the evaluation parameters.

The relationship between p and q
From expressions (22), (23), (15), (33), (34) and (28), it is seen that if ξ 1 and ξ 2 are definite, the values of the evaluation parameters should be definite for the above two Cantor sets. However, from expressions (24)-(26) and (35)-(37), it can be seen that the required calculation range of q is related not only to ξ 1 and ξ 2 but also to p.
Let p min and p max indicate the minimum and the maximum distribution probability in the generators of the Cantor sets, respectively. For the Cantor set p/0/1 − p,  For the Cantor set p/1 − 2p/p, It can be seen from expressions (21), (24)-(26), (32) and (35)-(37) that all q, q 1 , q 2 and q are proportional to |ln(p max /p min )| −1 . It means that when the partial multifractal spectra have the same evaluation parameters (or have the same described region relative to the complete spectrum), the smaller the difference of the distribution probability in the generators, the wider the range of q in calculating the multifractal spectra using the statistical method. |ln(p max /p min )| −1 is called the inhomogeneous effect factor. The two Cantor sets have the same expressions for α/ α 0 and q, i.e. α/ α 0 = ξ and q = ln p max p min Therefore, if ξ 1 , ξ 2 and p max /p min take the same values for the Cantor sets when calculating the partial multifractal spectrum using the statistical method, α/ α 0 and q will take the same values as well. The multifractal spectra for the Cantor sets 0.3/0.4/0.3 and 0.333/0.334/0.333 are shown in figures 2(a) and (b), respectively. The solid lines indicate the complete multifractal spectra obtained from the generators and the filled circles represent the results calculated using the statistical method. It can be seen from figure 2 that the required calculation range of q to get the complete multifractal spectrum is quite different for the two cases. There is little difference between the complete multifractal spectrum and the partial multifractal spectrum as |q| 20 in figure 2(a). In contrast, the partial multifractal spectrum for |q| 100 in figure 2(b) only describes a very small part of the complete multifractal spectrum. The large difference in the required calculation range of q comes from very different value of p max /p min for the two cases, which is 4 3 in figure 2(a) and 1.003003 in figure 2(b), respectively. The inhomogeneous effect factor is 3.476 and 333.5 for the two cases, respectively. The ratio of the inhomogeneous effect factors for the two cases is about two orders of magnitude. It leads to that the difference between two calculation ranges of q is also about two orders of magnitude. Sun et al [21] discussed the qualitative relationship between the inhomogeneous degree and q in the saturated regions using examples in the two Cantor sets. They found the smaller the inhomogeneous degree, the larger the q.

Approximate evaluation parameters
For a regular fractal, α min , α max , f 01 and f 02 are obtained from the generator. While for a random fractal, α min , α max , f 01 and f 02 are approximate. On the basis of the practical situation, two approaches can be used to estimate α min and α max . However, the approximate f 01 and f 02 cannot be found out if the errors of f(α) are too large.
Approach 1: In fact, α and the f(α) tend to be saturated when |q| is large enough. In the saturated regions, α is selected as the approximate α min and α max , and q relative to α min and α max are denoted as q max and q min , respectively.
Approach 2: Let l k represent the minimum of size l, and p k min and p k max be the minimum and the maximum distribution probability of the quantity studied, respectively as l = l k . Then, the approximate α min and α max can be written as α min = ln p k max /ln l k and α max = ln p k min /ln l k , respectively.
For a random fractal, we cannot get quantitative expressions to describe the range of α and f(α) as a function of q 1 and q 2 in a finite range of q using the statistical method. Moreover, we cannot find quantitative expressions about the required range of q related to the inhomogeneous effect factor also. However, the approximate evaluation parameters of a partial multifractal spectrum can be obtained. On the basis of expression (1), the approximate α 0 indicates the inhomogeneous degree of the fractal object studied. From expression (13) (for the Cantor set p/1 − 2p/p) and expression (27) (for the Cantor set p/0/1 − p), we get α 0 ∝ [ln(p max /p min )]. Therefore, the qualitative relationship between the inhomogeneous degree and the required range of q can still be studied.

Example 1: Hang Seng index in the Hong Kong stock market
The data are taken from the Hang Seng index in the Hong Kong stock market in a time period of 30 continuous trading days except holidays from 3 January 1994. There are totally 6240 indexes. Its multifractal spectrum is calculated using the statistical method. Figure 3(a) shows the α-q curve with |q| 300. In the two saturated regions shown in figure 3(a), q max = 150 and q min = −150 are selected using approach 1. The corresponding α is taken as the approximate minimum and maximum α of the complete multifractal spectrum, respectively. Figure 3(b) represents the f(α)-α curve (filled squares) and the σ τ -α curve (filled circles), as |q| 150.

Example 2: A photograph of a soybean leaf
The data are taken from the grey-scale values of a section of a Jindou 11 leaf's photograph. The Jindou 11 is one species of soybean, which originates from Shanxi province of China. There are 840 × 840 pixels. The multifractal spectrum is calculated using the statistical method. Figure 4(a) shows the α-q curve as |q| 30. q max = 10 and q min = −20 are selected from the two saturated regions shown in figure 4(a), respectively, according to approach 1. Figure 4(b) displays the f(α)-α curve (filled squares) and the σ τ -α curve (filled circles) as −20 q 10.

Discussion
Now, we arrive at the following conclusions by comparison.   (1) In the two examples, the evaluation parameters depend on errors. It can be seen from figures 3(b) and 4(b) that except for very small σ τ near q = 0, the σ τ increases rapidly with increasing |q|, and both σ α and σ f also increase with σ τ according to equation (9). As 2q 2 1, σ f ≈ √ 2qσ τ . σ f is proportional to |q| and increases significantly. For a partial multifractal spectrum, only the data whose error satisfies the requirement are useful. If the conditions, i.e. σ τ 2.0 × 10 −4 and σ α 2.8 × 10 −4 are fulfilled, q 1 and q 2 are restricted to -37 and 39 (−37 q 39) in figure 3(b) and −2.0 and 6.2 (−2.0 q 6.2) in figure 4(b), respectively. Therefore, using approach 1, the approximate evaluation parameters α 1 / α 01 and α 2 / α 02 are 53 and 50% for example 1 and 92 and 25% for example 2, respectively. For example 1, the corresponding σ f 1.5 × 10 −2 , but for example 2, σ f 2.5 × 10 −3 . Obviously, the calculated results are more reliable and give more information in section I (α < α 0 ) for example 2 when compared with those for example 1.
(2) The approximate α 0 for example 2 (the soybean leaf) is over 10 times as large as that for example 1 (the Hang Seng index), i.e. the inhomogeneous degree is smaller for the Hang Seng index than that for the soybean leaf. It is clear from figures 3(b) and 4(b) that if the evaluated parameters in the partial multifractal spectra are the same for the two examples, the required range of q is much wider for the Hang Seng index when compared with that for the soybean leaf. It is similar to the Cantor sets p/1 − 2p/p and p/0/1 − p.

Summary
A set of the evaluation parameters is introduced to describe partial multifractal spectra with the shape of a bell or a hook. The incompleteness of the multifractal spectra calculated using the statistical method comes from the data overflow in the computation, the measurement errors and the finite time of computation. An effective method to reduce the data overflow is proposed. Quantitative expressions of the evaluation parameters for the Cantor setsp/1 − 2p/p and p/0/1 − p are deduced and discussed. For the random fractals, two approaches of computing the approximate evaluation parameters are introduced and two examples are given. For the partial multifractal spectra with the same evaluation parameters using the statistical method, the required range of q in the computation is closely related to the inhomogeneous degree of the fractal object. The smaller the inhomogeneous degree, the larger the required range of q. The α-q and σ τ -α curves are helpful in selecting a reasonable range of q in the computation.