Bio-Medical Materials and Engineering 24 (2014) 1041–1051 DOI 10.3233/BME-130902 IOS Press
1041
A Hurst exponent estimator based on autoregressive power spectrum estimation with order selection Yen-Ching Chang a,*, Li-Chun Lai b, Liang-Hwa Chenc, Chun-Ming Chang d and Chin-Chen Chuehe a
Department of Medical Informatics, Chung Shan Medical University and Department of Medical Imaging, Chung Shan Medical University Hospital, Taichung, 40201, Taiwan, ROC b Bachelor Program in Robotics, National Pingtung University of Education, Pingtung, 90003, Taiwan, ROC c Department of Computer Information and Network Engineering, Lunghwa University of Science and Technology, Taoyuan, 33306, Taiwan, ROC d Department of Applied Informatics and Multimedia, Asia University, Taichung, 41354, Taiwan, ROC e Department of Medical Informatics, Chung Shan Medical University, Taichung, 40201, Taiwan, ROC
Abstract. The discrete-time fractional Gaussian noise (DFGN) has been proven to be a regular process. According to Wold and Kolmogorov theorems, this process can be described as an autoregressive (AR) model of an infinite order. An estimator for the Hurst exponent based on autoregressive power spectrum estimation has been proposed, but without considering order selection. In this paper, six common order selection methods for the AR model were used to select appropriate orders of the AR model in order to raise the accuracy of estimating the Hurst exponent. Experimental results show that these six AR methods with considering order selection are more accurate than the original AR method without considering order selection. Keywords: Fractional Brownian motion, fractional Gaussian noise, Hurst exponent, autoregressive, order selection
1. Introduction In the fields of nature and biology, signals usually exhibit a strong long-term correlation and selfsimilarity [1–26]. Among them, Pentland addressed the problems of representing natural shapes and computing their description from image data [4]. Jennane et al. studied the fractal analysis of bone Xray tomographic microscopy projections [16]. Peréz et al. used the estimation of Hurst parameter and fractal dimension to analyze the hemogram components [17]. Esgiar et al. used the fractal dimension to separate the normal and cancerous images [18]. Chang et al. applied the fractal dimension to analyze the dynamics and synchronization of rhythms urodynamics of female Wistar rats [20]. Lee et al. used the fractal dimension of a skeletonized cerebral surface to analyze the hemispheric asymmetry *
Corresponding author. E-mail: [email protected].
0959-2989/14/$27.50 © 2014 – IOS Press and the authors. All rights reserved
1042
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
[21]. For the purpose of analysis, it is often necessary to model these signals. Among all models that are suitable for description, the fractional Brownian motion (FBM) [7, 10] is well known and has been invoked by engineers who engage in the field of statistical signal processing. The simplicity of FBM lies in the fact that a simple parameter, i.e., the Hurst exponent, with its values ranging from 0 to 1, can characterize a wide choice of biomedical and natural signals. Another model is called the fractional Gaussian noise (FGN), which is an increment process of FBM. On the contrary, FBM is an accumulative process of FGN. In practical applications, data are usually sampled from a continuous-time process. The sampled FBM is referred to as the discrete-time FBM (DFBM), and the sampled increment process of FBM as the discrete-time FGN (DFGN). There are many different methods estimating the Hurst exponent of DFBM and DFGN, including the variance method [1, 5, 27], box-counting method [8–9, 11–12], maximum likelihood estimator (MLE) [5], moving-average (MA) method [14], and autoregressive (AR) method [19]. The variance method needs to consider the order of linear regression, or its accuracy is weaker and even less reliable. Chang et al. has proved that the variance estimator can use the autocorrelation function (ACF) to estimate the Hurst exponent, and that 4-point linear regression is almost optimal in most cases [27]. Box-counting methods are rough and need to deal with the problems of scales or box sizes. Lundahl et al. [5] proposed an MLE method and applied it to the texture of the X-ray images of bones. This method possesses the optimal accuracy but requires substantial repeated matrix operations during estimation. Since its computational cost is very high, the MLE is not recommended to systems that need quick response. To improve the execution efficiency of the MLE, Liu and Chang [14] proposed a moving-average method and further applied it to the classification problem of image texture. The MA method exploits Wold decomposition theorem to decompose the DFGN into two uncorrelated processes: one is regular and the other is singular [28–30]. After the ACFs are calculated, Levinson algorithm [31–32] is used to estimate the parameters of the AR model, and then the algorithm [33–34] of Wold decomposition theorem is utilized to estimate the parameters of the MA model. This MA method is indeed faster than the MLE, but less accurate. If DFGN can be shown to be a regular process, then the transformation of the AR model into the MA one is no longer necessary, which will significantly reduce the computational cost. For the purpose, Chang and Chang employed a comparative method to show that the DFGN is indeed a regular process [19]. This comparative method is based on the fact that the ACFs between DFGN and fractionally differenced Gaussian noise (fdGn) [35] are similar and their approximate behavior on the ACFs is nearly the same. Thus, a more accurate and efficient estimator than the MA method, called the AR method, was proposed. The AR method was proved to be an “approximate” MLE [14, 31, 36]. In theory, an AR model requires an infinite order to describe a regular process. Unfortunately, the given data length is always finite and much larger errors will be induced by higher order estimates through the estimated ACFs. Therefore, it is essential for the AR method to take order selection into consideration. In this paper, six common order selection methods are considered to select the order of the AR model. The AR model with order selection can be used to derive a more accurate algorithm than the original AR model without considering order selection [19]. The rest of this paper is organized as follows. Section 2 briefly describes mathematical preliminaries. Section 3 provides the procedures of estimating the Hurst exponent and its application. Results are discussed in Section 4. Finally, the conclusions are drawn in Section 5.
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1043
2. Mathematical preliminaries In this section, we will introduce some related terminologies and notations for the rest of the paper. Besides, Levinson algorithm and common model order selection methods will be briefly introduced. The notation {x(t ), t ∈ R} is used to denote a continuous-time random process, and {x[n ], n ∈ Z } a discrete-time process. 2.1. FBM and its variants The FBM is a generalized Brownian motion, first officially proposed by Mandelbrot and Van Ness [29] and denoted as the random function BH (t , ω ) , where ω belongs to a sample space Ω , simplified as BH (t ) . It can be defined as follows: BH (0 ) = b0 , BH (t2 ) − BH (t1 ) =
1 Γ (H + 1 2 )
{³
t2
−∞
}
(t2 − s )H −1 2 dB(s ) − ³−∞ (t1 − s )H −1 2 dB(s ) , t1
(1)
where H is the Hurst exponent, the only parameter with a value lying between 0 and 1. It is easy to understand that FBM will become an ordinary Brownian motion at H = 0.5 . Unfortunately, FBM is not a stationary process, i.e., its spectrum is time-dependent. For obtaining a time-independent spectrum, its increment of FBM, denoted as BH′ (t ) and called FGN, is considered. FGN is a stationary process and its spectrum is expressed as follows: S BH′ (Ω ) =
1 Ω
2 H −1
.
(2)
Furthermore, the discrete version of FBM, called the DFBM, is denoted as BH [n] = BH (nTs ) , where Ts denotes the sampling time. The DFGN is denoted as x H [n] = BH [n] − BH [n − 1] and its corresponding ACF as rH [k ] =
(
)
VH 2H 2H 2H , k +1 − 2 k + k −1 2
(3)
where VH = Γ(1 − 2 H ) cos(Hπ ) Hπ [1, 6]. The ACF rH [k ] behaves asymptotically as k 2 H −2 = k −α , α ∈ (0, 2 ) [10]. Chang and Chang [19] have proved the DFGN is a regular process using comparisons between two types of ACFs of DFGN and fdGn. Because of the regular property, an AR model can be used to describe the behavior of DFGN. Furthermore, Levinson algorithm is utilized to efficiently obtain the coefficients of an AR model. Some of its characteristics have been summarized by Lundahl et al. [5] Moreover, Samorodnitsky and Taqqu [10] have a complete introduction to these topics. 2.2. Levinson algorithm The Levinson algorithm [31], or Levinson-Durbin recursion [32], is an efficient method for solving the Yule-Walker equation or the Wiener-Hopf equation. Its computational cost can further be reduced by using the Schur algorithm [32], in which inner products are no longer required. This algorithm can provide the AR parameters for all lower-order AR models fitting to the ACFs. This property is useful, especially when users do not know a prior the correct model order. If the AR model is finite and free
1044
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
of noise, one can use this algorithm to derive the correct model order and its corresponding parameters according to the minimum variation of the prediction error power. However, if the AR model order is infinite or the model is finite but contaminated by noise, the algorithm alone might not choose the correct order. For later notation need and implementation, the Levinson algorithm [31] is summarized as follows. The Levinson algorithm recursively computes the parameter sets {a1 [1], ρ1} , {a2 [1], a2 [2], ρ 2 } , …, {a p [1], a p [2],! , a p [ p ], ρ p }. The final set at order p is the desired solution of the Yule-Walker equa-
tions. If x[n] is an AR of order p , then a p [i ] = a [i ] for i = 1, 2,! , p and ρ p = σ 2 , where σ 2 is the excitation noise variance. The algorithm is first initialized by a1 [1] = −
rxx [1] rxx [0]
(4)
and
ρ1 = (1 − a12 [1]) rxx [0] ,
(5)
where rxx is the ACF of x[n] , and then the following parameter sets are computed by k −1
ak [k ] = −
rxx [k ] + ¦ ak −1 [l ]rxx [k − l ] l =1
ρ k −1
,
(6)
ak [i ] = ak −1 [i ] + ak [k ]ak −1 [k − i ] , i = 1, 2, ! , k − 1 ,
(7)
ρ k = (1 − ak2 [k ]) ρ k −1 ,
(8)
and
for k = 2, 3, !, p . 2.3. Model order selection Six common methods for order selection of the AR model include the final prediction error (FPE) [31–32, 37–38], Akaike information criterion (AIC) [31–32, 37–38], criterion autoregressive transfer function (CAT) [31, 37–38], minimum description length (MDL) criterion [32, 37–38], residual variance (RV) [37–38], as well as Hannan and Quinn (HQ) [37–38]. The FPE estimates the model order as the place minimizing the following equation: FPE (k ) =
N +k ρˆ k , N −k
(9)
where ρˆ k is the estimate of the white noise variance (prediction error power) [31], or called the average prediction-error power [32], for the kth order AR model. The second method is the AIC, which is defined as the following equation: AIC(k ) = N ln ρˆ k + 2k ,
(10)
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1045
As before, minimizing the AIC gives the order selected. The performance of the AIC is similar to the one of the FPE. For larger data records, the two criteria will produce almost identical model order results because they are mathematically related to each other. The third method is the MDL and is defined as MDL(k ) = N ln ρˆ k + k ln N .
(11)
The MDL is asymptotically consistent unlike the AIC, which is not a consistent estimator. Similarly, the model order estimate is obtained by minimizing the MDL. The fourth method, the CAT, is defined as CAT(k ) =
1 N
k
1
¦ ρ~ i =1
i
1 − ~ , ρ
(12)
k
where
ρ~i =
N ρˆ i . N −i
(13)
The model order estimate can be chosen by minimizing the CAT. The fifth method, called the RV, is defined as RV (k ) =
N −k ρˆ k . N − 2k − 1
(14)
Similarly, the model order estimate is obtained by minimizing the RV. The sixth method, called the HQ, is defined as HQ(k ) = N ln ρˆ k + 2k ln ln N .
(15)
The model order is estimated as the position of the minimum of the HQ. These methods work well on finite-order situations. However, according to Wold and Kolmogorov theorems, the AR model order for a regular process is infinite. Therefore, we would like to know whether or not they are applicable to an infinite order model. 3. Procedures of estimating the Hurst exponent and its application It has been shown that the DFGN is a regular process [19], and hence the transformation of the AR model into the MA model can be omitted. Experimental results have shown that the AR method is superior to the MA method in accuracy and efficiency. Since the AR model of the DFGN is an infinite order, the higher the model order is, the closer to a regular process it is. It can be seen from the true ACFs, where the true ACF rH [k ] behaves asymptotically as k 2 H −2 = k −α , α ∈ (0, 2 ) , never a zero for any lag [10]. Accordingly, the order of the AR model for estimated ACFs was chosen as the maximum possible order by intuition, i.e., the data length minus one [19]. However, for finite data, errors easily happen at higher orders because the AR model highly depends on inaccurate estimates of the ACFs. These inaccurate estimates further result in adverse effects on the estimation of the Hurst exponent. Therefore, a trade-off between theoretical demand and practical estimation will be essential. In this paper, the standard or referenced curve of the slopes of power spectral density (PSD) versus
1046
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
the Hurst exponent was obtained using the true ACFs, hence the higher the model order is, the closer to a true model it is. Unfortunately, an infinite order cannot be modeled in a finite time-space, and estimation is a very complex process. In order to simplify the complexity of comparison, a 1000 order of AR model was adopted in the paper and a fitting curve of 4 degrees was derived. The estimated slopes were obtained by the estimated ACFs based on the model order chosen from six model order selection methods. Finally, an estimated Hurst exponent was produced through an established standard curve ( H versus slope) beforehand, described as follows. First, every true Hurst exponent was used to produce an estimated slope of AR power spectrum estimates versus frequencies in a logarithmically scaled way sˆ = f (H ) , H ∈ {0.01, 0.02, ! , 0.99} ,
(16)
and then the standard curve was established using least squares polynomial curve-fitting [39] of 4 degrees as follows: Hˆ = c 4 sˆ 4 + c3 sˆ 3 + c2 sˆ 2 + c1 sˆ1 + c0 sˆ 0 ,
(17)
c0 = 0.4998 , c1 = -0.4327 , c2 = 0.0400 , c3 = 0.0192 , and c4 = 0.0040 .
(18)
where
The detailed steps are listed in the following procedure: Procedure 1: Production of one standard curve for N = 1000 1. 2. 3. 4.
Begin the Hurst exponent with H = 0.01 with increment 0.01. Calculate the true ACFs from (3). Estimate the AR parameters and the white noise variance (prediction error power) using these true ACFs via Levinson algorithm or from (4) – (8). Calculate Aˆ ( f ) = 1 + aˆ [1]exp {− j 2πf } + " + aˆ [ pˆ ]exp {− j 2πfpˆ } ,
5.
(19)
where pˆ is the estimated order of the AR model. For true ACFs, the order of the AR model is selected as N − 1 . Estimate the PSD of the AR model via the following equation 2 Sˆ ( f ) = ρˆ pˆ | Aˆ ( f ) | ,
(20)
where ρˆ pˆ is the estimate of the white noise variance (prediction error power) for order pˆ . 6.
7. 8.
{ (
)}
Estimate the slope, sˆ , of log Sˆ ( f i ) versus {log( f i )} over FM = { f1 , f 2 , ! , f M } , where f1 = 0.5 M , f i − f i −1 = 0.5 M for i = 2, 3, !, M with M being selected as 63 in this paper. If H = 0.99 , then go to next; or execute the next H . Use the least squares polynomial curve-fitting of 4 degrees to fit these 99 points of the Hurst exponent versus the estimated slope.
Once a standard curve established, a procedure for estimating the Hurst exponent and its performance evaluation is provided to determine whether or not the model order selection method can con-
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1047
tribute to the accuracy of estimation. For the purpose, six model order selection methods were chosen, including the FPE, AIC, MDL, CAT, RV, and HQ; six types of data length, including 128, 256, 512, 1024, 2048, and 4096; 13 Hurst exponents, including 0.01, 0.05, 0.1, 0.2, ! , 0.9 , 0.95, and 0.99. Then, 100 realizations of DFGN, x H [n ] , for different data length and Hurst exponent were simulated according to the algorithm of Lundahl et al. [5] The algorithm provides a fine correlation structure and longterm dependency. The detailed steps are listed in the following procedure: Procedure 2: Estimation of Hurst exponent and its performance evaluation
1. 2. 3. 4. 5. 6.
Begin the model order selection method with the FPE. Begin the data length with N = 128 . Begin the Hurst exponent with H = 0.01 . Input the first realization of x H [n ] for N and H . Modify each realization of x H [n ] into zero mean. Estimate ACFs for k = 0,1, ! , N − 1 using the following form rˆ[ k ] =
1 N
N
¦ x[n]x[n − k ] .
(21)
n = k +1
Estimate the AR parameter sets using these estimated ACFs via Levinson algorithm or from (4) – (8), and then memorize these parameters. 8. Use the model order selection method to choose the optimal order, pˆ . 9. Extract the corresponding AR parameter set with order pˆ from memory. 10. Estimate the PSDs of AR model using (19) and (20). 11. Calculate the slopes of log Sˆ ( f i ) versus {log( f i )} over FM = { f1 , f 2 , ! , f M } , where f1 = 0.5 M , f i − f i −1 = 0.5 M for i = 2, 3, ! , M with M being selected as 63. 12. Estimate the H through the established standard curve (17). 13. If the number of realization is 100, then go to next; or execute the next realization. 14. Compute the mean mean-squared error between the estimated H and true H over 100 realizations. 15. If H = 0.99 , then go to next; or execute the next H . 16. Compute the mean mean-squared error over these 13 Hurst exponents. 17. If N = 4096 , then go to next; or execute the next N . 18. If the model order selection method is equal to the HQ, then Stop; or execute the next method. 7.
{ (
)}
From Procedure 2, one relatively better model order selection method will be discovered. Users can use the recommended method of order selection to more effectively estimate the Hurst exponent in the signals of interest. For easy use, a complete procedure is described as follows: Procedure 3: Use of the proposed estimator to a practical signal
1. 2. 3.
Input a section of practical signal, x[n] . Modify x[n] into zero mean. Estimate ACFs for k = 0,1, ! , N − 1 using (21).
1048
4.
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
8.
Estimate the AR parameter sets using these estimated ACFs via Levinson algorithm or from (4)–(8), and then memorize these parameters. Use the model order selection method, for example, the RV, to choose the optimal order, pˆ . Extract the corresponding AR parameter set with order pˆ from memory. Estimate the PSDs of AR model using (19) and (20). Calculate the slopes of log Sˆ ( f ) versus {log( f )} over F = { f , f , ! , f } , where
9.
f1 = 0.5 M , f i − f i −1 = 0.5 M for i = 2, 3, ! , M with M being selected as 63. Estimate the H through the established standard curve (17).
5. 6. 7.
{ (
i
)}
i
M
1
2
M
4. Results and discussion
In order to verify the effectiveness of the estimator proposed in this paper, a Hurst exponent estimator based on AR model with order selection, two sets of data were simulated, each producing 100 realizations of DFGN. The purpose of another data set is to confirm that experimental results are not affected by different realizations. In addition, the original AR method without model order selection, simply called the AR method, was also implemented to compare. Other six model order selection methods considered into the AR model will be simply referred to as the FPE, AIC, MDL, CAT, RV, and HQ, not the AR-FPE, AR-AIC, etc. Their experimental results are listed in Tables 1–6. Table1 provides the mean mean-squared error for each data length over 13 Hurst exponents and 100 realizations on Set 1; Table 2, the mean values of orders selected for each data length over 13 Hurst exponents and 100 realizations on Set 1; Table 3, the standard deviation values of orders selected for each data length over 13 Hurst exponents and 100 realizations on Set 1. Tables 4–6 are the results of Set 2 corresponding to Tables 1–3. It is obvious from Tables 1 and 4 that the AR methods with considering order selection are more accurate than the original one without considering order selection. Moreover, the RV method is almost the best except for N = 4096 . The FPE, AIC, and CAT are almost equivalent to one another. This phenomenon is consistent with the description that the performance of the CAT has been found to be similar to that of the FPE and AIC [31]. The MDL is the worst among the six methods. As described by Haykin [32], the application of the FPE criterion has been found to produce excellent results in the case of a pure AR process. However, when the FPE criterion is applied to practical non-AR data, the criterion appears to be conservative in the sense that it often yields too low a value for the model order. As with the FPE criterion, the AIC has been found to select too low a value for the model order when it is applied to practical non-AR data. Experimental results also confirm this argument. The FPE and AIC indeed choose a lower order than the RV. Table 1 Comparison of methods for mean mean-squared error on Set 1 Methods
N = 128
N = 256
N = 512
N = 1024 N = 2048 N = 4096
AR FPE AIC MDL CAT RV HQ
3.37E-03 3.54E-03 3.54E-03 7.14E-03 3.55E-03 2.80E-03 4.89E-03
3.06E-03 1.38E-03 1.38E-03 3.27E-03 1.39E-03 1.21E-03 2.18E-03
3.41E-03 8.72E-04 8.72E-04 2.13E-03 8.74E-04 6.55E-04 1.33E-03
3.29E-03 4.64E-04 4.64E-04 1.09E-03 4.65E-04 3.71E-04 6.83E-04
3.01E-03 2.29E-04 2.29E-04 6.53E-04 2.30E-04 1.83E-04 3.52E-04
3.00E-03 1.36E-04 1.36E-04 4.15E-04 1.36E-04 1.48E-04 2.35E-04
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1049
Table 2 Comparison of methods for mean values of orders selected on Set 1 Methods FPE AIC MDL CAT RV HQ
N = 128 3.02 3.03 1.38 2.98 6.42 1.90
N = 256 4.51 4.52 1.87 4.43 10.74 2.70
N = 512 6.64 6.64 2.69 6.54 17.28 3.81
N = 1024 9.36 9.36 3.78 9.26 29.72 5.48
N = 2048 12.57 12.57 5.00 12.51 55.73 7.33
N = 4096 17.69 17.69 6.93 17.67 74.40 10.01
Table 3 Comparison of methods for standard deviation values of orders selected on Set 1 Methods FPE AIC MDL CAT RV HQ
N = 128 2.22 2.22 0.66 2.17 4.65 1.00
N = 256 2.60 2.61 0.78 2.53 7.13 1.34
N = 512 3.49 3.49 1.00 3.38 10.51 1.50
N = 1024 4.39 4.39 1.11 4.30 20.36 1.90
N = 2048 4.98 4.98 1.15 4.95 35.31 2.06
N = 4096 6.08 6.08 1.47 6.08 45.65 2.50
On the other hand, the AIC is statistically inconsistent in the sense that the probability of selecting the correct order does not approach unity as the data length N approaches infinity. In fact, there is a tendency for the AIC to overestimate the model order in the case of a pure AR process as N approaches infinity. Although the MDL was developed to overcome the asymptotically inconsistent behavior of the AIC, its purpose cannot work on this infinite order of AR model. Tables 2 and 5 show that the MDL excessively underestimates the model order in the case of an infinite AR process. Table 4 Comparison of methods for mean mean-squared error on Set 2 Methods
N = 128
N = 256
N = 512
AR FPE AIC MDL CAT RV HQ
4.61E-03 3.25E-03 3.12E-03 4.22E-03 1.88E-03 9.90E-04 4.21E-03 1.88E-03 9.90E-04 7.59E-03 4.09E-03 2.19E-03 4.20E-03 1.88E-03 9.90E-04 3.54E-03 1.65E-03 8.36E-04 5.32E-03 2.77E-03 1.35E-03
N = 1024 N = 2048 N = 4096 3.18E-03 3.05E-03 4.24E-04 2.39E-04 4.24E-04 2.38E-04 1.01E-03 6.47E-04 4.25E-04 2.40E-04 3.46E-04 1.71E-04 6.11E-04 3.79E-04
2.86E-03 1.38E-04 1.38E-04 4.15E-04 1.38E-04 1.59E-04 2.42E-04
Table 5 Comparison of methods for mean values of orders selected on Set 2 Methods FPE AIC MDL CAT RV HQ
N = 128
N = 256
N = 512
3.33 3.34 1.38 3.26 6.35 1.88
4.67 4.69 1.89 4.55 11.11 2.69
6.44 6.44 2.72 6.41 17.01 3.85
N = 1024 N = 2048 N = 4096 9.09 9.10 3.71 9.00 28.89 5.31
12.53 12.54 5.06 12.47 48.17 7.34
18.76 18.76 6.81 18.64 89.84 10.13
1050
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
Table 6 Comparison of methods for standard deviation values of orders selected on Set 2 Methods FPE AIC MDL CAT RV HQ
N = 128 2.98 2.98 0.64 2.83 4.55 1.01
N = 256 3.07 3.12 0.82 2.79 7.64 1.27
N = 512 3.17 3.17 1.05 3.14 10.85 1.60
N = 1024 N = 2048 N = 4096 4.34 5.12 7.12 4.35 5.13 7.12 1.08 1.20 1.49 4.25 5.06 7.05 19.56 32.60 55.81 1.84 2.05 2.58
5. Conclusions
In theory, an infinite order of autoregressive (AR) model can be used to model the discrete-time fractional Gaussian noise and its increment process, the discrete-time fractional Brownian motion (DFBM). Therefore, in a signal of DFGN or DFBM with finite data length, the order of an AR model should be chosen to be as large as possible by intuition. Unfortunately, the estimation of Hurst exponent using an AR model highly depends on the estimation of autocorrelation functions (ACFs), which easily leads to larger errors in larger lags. In order to resolve the dilemma between a theoretical requirement and an estimation difficulty, six common order selection methods were introduced to choose a suitable order for the AR model. Experimental results show that the estimation of the Hurst exponent with considering order selection into the AR model outperforms the original one without considering order selection. Moreover, the residual variance (RV) method is almost the best among these six methods except in the data length of 4096. 6. Acknowledgements
This work was supported by the National Science Council, Republic of China, under Grant NSC 101-2221-E-040 -010. References [1]
B.B. Mandelbrot and J.W.V. Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (1968), 422-437. [2] B.B. Mandelbrot, Fractals: Form, Chance and Dimension, W. H. Freeman, New York, 1977. [3] B.B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman, New York, 1983. [4] A.P. Pentland, Fractal-Based Description of Natural Scenes, IEEE Trans. Pattern Anal. Machine Intell. PAMI-6 (1984), 661-674. [5] T. Lundahl, J. Ohley, S.M. Kay, and R. Siffert, Fractional Brownian motion: A maximum likelihood estimator and its application to image texture, IEEE Trans. Med. Image, MI-5 (1986), 152-161. [6] P. Flandrin, On the Spectrum of Fractional Brownian Motions, IEEE Trans. Inform. Theory 35 (1989), 197-199. [7] K. Falconer, Fractal Geometry: Mathematical Foundations and Applications. John Wiley & Sons, New York, 1990. [8] N. Sarkar and B.B. Chaudhuri, An efficient approach to estimate fractal dimension of textural images, Pattern Recognit. 25 (1992), 1035-1041. [9] S.S. Chen, J.M. Keller, and R.M. Crownover, On the calculation of fractal features from images, IEEE Trans. Pattern Anal. Mach. Intell. 15 (1993), 1087-1090. [10] G. Samorodnitsky and M.S. Taqqu, Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance, Chapman & Hall, New York, 1994.
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1051
[11] N. Sarkar and B.B. Chaudhuri, An efficient differential box-counting approach to compute fractal dimension of image, IEEE Trans. Syst. Man Cybern. 24 (1994), 115-120. [12] X.C. Jin, S.H. Ong, and Jayasooriah, A practical method for estimating fractal dimension, Pattern Recognit. Lett. 16 (1995), 457-464. [13] M.S. Taqqu, V. Teverovsky, and W. Willinger, Estimators for long-range dependence: An empirical study, Fractals 3 (1995), 785-798. [14] S.-C. Liu and S. Chang, Dimension estimation of discrete-time fractional Brownian motion with applications to image texture classification, IEEE Trans. Image Process. 6 (1997), 1176-1184. [15] S. Chang, S.-T. Mao, S.-J. Hu, W.-C. Lin, and C.-L. Cheng, Studies of detrusor-sphincter synergia and dyssynergia during micturition in rats via fractional Brownian motion, IEEE Trans. Biomed. Eng. 47 (2000), 1066-1073. [16] R. Jennane, W.J. Ohley, S. Majumdar, and G. Lemineur, Fractal analysis of bone X-ray tomographic microscopy projections, IEEE Trans. Med. Imaging 20 (2001), 443-449. [17] A. Peréz, C.E. D’Attellis, M. Rapacioli, G.A. Hirchoren, and V. Flores, Analyzing blood cell concentration as a stochastic process, IEEE Eng. Med. Biol. Mag. 20 (2001), 170-175. [18] A.N. Esgiar, R.N.G. Naguib, B.S. Sharif, M.K. Bennett, and A. Murray, Fractal analysis in the detection of colonic cancer images, IEEE Trans. Inf. Technol. Biomed. 6 (2002), 54-58. [19] Y.-C. Chang and S. Chang, A fast estimation algorithm on the Hurst parameter of discrete-time fractional Brownian motion, IEEE Trans. Signal Process. 50 (2002), 554-559. [20] S. Chang, S.-J. Hu, and W.-C. Lin, Fractal dynamics and synchronization of rhythms in urodynamics of female Wistar rats, J. Neurosci. Meth. 139 (2004), 271-279. [21] J.-M. Lee, U. Yoon, J.-J. Kim, I.Y. Kim, D.S. Lee, J.S. Kwon, and S.I. Kim, Analysis of the hemispheric asymmetry using fractal dimension of a skeletonized cerebral surface, IEEE Trans. Biomed. Eng. 51 (2004), 1494-1498. [22] S. Chang, S.-J. Li, M.-J. Chiang, S.-J. Hu M.-C. Hsyu, Fractal dimension estimation via spectral distribution function and its application to physiological signals, IEEE Trans. Biomed. Eng. 54 (2007), 1895-1898. [23] S. Chang, M.-C. Hsyu, H.-Y. Cheng, and S.-H. Hsieh, Synergic co-activation of muscles in elbow flexion via fractional Brownian motion, Chin. J. Physiol. 51 (2008), 376-386. [24] S. Chang, M.-C. Hsyu, H.-Y. Cheng, S.-H. Hsieh, and C.-C. Lin, Synergic co-activation in forearm pronation, Ann. Biomed. Eng. 36 (2008), 2002-2018. [25] S. Chang, Physiological rhythms, dynamical diseases and acupuncture, Chin. J. Physiol. 53 (2010), 77-90. [26] S. Chang, Fractional Brownian motion in biomedical signal processing, physiology, and modern physics, The 5th ICBBE, May 10-12, 2011; Wuhan, China. [27] Y.-C. Chang, L.-H. Chen, L.-C. Lai, and C.-M. Chang, An efficient variance estimator for the Hurst exponent of discrete-time fractional Gaussian noise, IEICE Trans. Fundam. Electron. Commun. Comput. Sci. E95-A (2012), 15061511. [28] J.L. Doob, Stochastic Processes, John Wiley & Sons, New York, 1953. [29] M.B. Priestley, Spectral Analysis and Time Series, Volume 2: Multivariate Series, Prediction and Control, Academic Press, New York, 1981. [30] A.N. Shiryaev, Probability, 2nd ed., translated by R. P. Boas, Springer-Verlag, New York, 1996. [31] S.M. Kay, Modern Spectral Estimation: Theory & Application, Prentice-Hall, New Jersey, 1988. [32] S. Haykin, Modern Filters, Macmillan, New York, 1989. [33] A. Papoulis, Predictable processes and Wold’s decomposition: A review, IEEE Trans. Acoust. Speech, Signal Process. 33 (1985), 933-938. [34] A. Papoulis, Levinson's algorithm, Wold’s decomposition, and spectral estimation, SIAM Rev. 27 (1985), 405-441. [35] M. Deriche and A.H. Tewfik, Maximum likelihood estimation of the parameters of discrete fractionally differenced Gaussian noise process, IEEE Trans. Signal Process. 41 (1993), 2977-2989. [36] K.L. Chung, A Course in Probability Theory, 2nd ed., Academic, New York, 1973. [37] R. Palaniappan, Towards optimal model order selection for autoregressive spectral analysis of mental tasks using genetic algorithm, International Journal of Computer Science and Network Security 6 (2006), 153-162. [38] A.M. Aibinu, A.R. Najeeb, M.J.E. Salami, and A.A. Shafie, Optimal model order selection for transient error autoregressive moving average (TERA) MRI reconstruction method, World Academy of Science, Engineering and Technology 42 (2008), 161-165. [39] D. Hanselman and B. Littlefield, Mastering MATLAB 6: A Comprehensive Tutorial and Reference, Prentice-Hall, New Jersey, 2001. [40] C.W.J. Granger and R. Joyeux, An introduction to long-memory time series models and fractional differencing, J. Time Ser. Anal. 1 (1980), 15-29.
Copyright of Bio-Medical Materials & Engineering is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
1041
A Hurst exponent estimator based on autoregressive power spectrum estimation with order selection Yen-Ching Chang a,*, Li-Chun Lai b, Liang-Hwa Chenc, Chun-Ming Chang d and Chin-Chen Chuehe a
Department of Medical Informatics, Chung Shan Medical University and Department of Medical Imaging, Chung Shan Medical University Hospital, Taichung, 40201, Taiwan, ROC b Bachelor Program in Robotics, National Pingtung University of Education, Pingtung, 90003, Taiwan, ROC c Department of Computer Information and Network Engineering, Lunghwa University of Science and Technology, Taoyuan, 33306, Taiwan, ROC d Department of Applied Informatics and Multimedia, Asia University, Taichung, 41354, Taiwan, ROC e Department of Medical Informatics, Chung Shan Medical University, Taichung, 40201, Taiwan, ROC
Abstract. The discrete-time fractional Gaussian noise (DFGN) has been proven to be a regular process. According to Wold and Kolmogorov theorems, this process can be described as an autoregressive (AR) model of an infinite order. An estimator for the Hurst exponent based on autoregressive power spectrum estimation has been proposed, but without considering order selection. In this paper, six common order selection methods for the AR model were used to select appropriate orders of the AR model in order to raise the accuracy of estimating the Hurst exponent. Experimental results show that these six AR methods with considering order selection are more accurate than the original AR method without considering order selection. Keywords: Fractional Brownian motion, fractional Gaussian noise, Hurst exponent, autoregressive, order selection
1. Introduction In the fields of nature and biology, signals usually exhibit a strong long-term correlation and selfsimilarity [1–26]. Among them, Pentland addressed the problems of representing natural shapes and computing their description from image data [4]. Jennane et al. studied the fractal analysis of bone Xray tomographic microscopy projections [16]. Peréz et al. used the estimation of Hurst parameter and fractal dimension to analyze the hemogram components [17]. Esgiar et al. used the fractal dimension to separate the normal and cancerous images [18]. Chang et al. applied the fractal dimension to analyze the dynamics and synchronization of rhythms urodynamics of female Wistar rats [20]. Lee et al. used the fractal dimension of a skeletonized cerebral surface to analyze the hemispheric asymmetry *
Corresponding author. E-mail: [email protected].
0959-2989/14/$27.50 © 2014 – IOS Press and the authors. All rights reserved
1042
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
[21]. For the purpose of analysis, it is often necessary to model these signals. Among all models that are suitable for description, the fractional Brownian motion (FBM) [7, 10] is well known and has been invoked by engineers who engage in the field of statistical signal processing. The simplicity of FBM lies in the fact that a simple parameter, i.e., the Hurst exponent, with its values ranging from 0 to 1, can characterize a wide choice of biomedical and natural signals. Another model is called the fractional Gaussian noise (FGN), which is an increment process of FBM. On the contrary, FBM is an accumulative process of FGN. In practical applications, data are usually sampled from a continuous-time process. The sampled FBM is referred to as the discrete-time FBM (DFBM), and the sampled increment process of FBM as the discrete-time FGN (DFGN). There are many different methods estimating the Hurst exponent of DFBM and DFGN, including the variance method [1, 5, 27], box-counting method [8–9, 11–12], maximum likelihood estimator (MLE) [5], moving-average (MA) method [14], and autoregressive (AR) method [19]. The variance method needs to consider the order of linear regression, or its accuracy is weaker and even less reliable. Chang et al. has proved that the variance estimator can use the autocorrelation function (ACF) to estimate the Hurst exponent, and that 4-point linear regression is almost optimal in most cases [27]. Box-counting methods are rough and need to deal with the problems of scales or box sizes. Lundahl et al. [5] proposed an MLE method and applied it to the texture of the X-ray images of bones. This method possesses the optimal accuracy but requires substantial repeated matrix operations during estimation. Since its computational cost is very high, the MLE is not recommended to systems that need quick response. To improve the execution efficiency of the MLE, Liu and Chang [14] proposed a moving-average method and further applied it to the classification problem of image texture. The MA method exploits Wold decomposition theorem to decompose the DFGN into two uncorrelated processes: one is regular and the other is singular [28–30]. After the ACFs are calculated, Levinson algorithm [31–32] is used to estimate the parameters of the AR model, and then the algorithm [33–34] of Wold decomposition theorem is utilized to estimate the parameters of the MA model. This MA method is indeed faster than the MLE, but less accurate. If DFGN can be shown to be a regular process, then the transformation of the AR model into the MA one is no longer necessary, which will significantly reduce the computational cost. For the purpose, Chang and Chang employed a comparative method to show that the DFGN is indeed a regular process [19]. This comparative method is based on the fact that the ACFs between DFGN and fractionally differenced Gaussian noise (fdGn) [35] are similar and their approximate behavior on the ACFs is nearly the same. Thus, a more accurate and efficient estimator than the MA method, called the AR method, was proposed. The AR method was proved to be an “approximate” MLE [14, 31, 36]. In theory, an AR model requires an infinite order to describe a regular process. Unfortunately, the given data length is always finite and much larger errors will be induced by higher order estimates through the estimated ACFs. Therefore, it is essential for the AR method to take order selection into consideration. In this paper, six common order selection methods are considered to select the order of the AR model. The AR model with order selection can be used to derive a more accurate algorithm than the original AR model without considering order selection [19]. The rest of this paper is organized as follows. Section 2 briefly describes mathematical preliminaries. Section 3 provides the procedures of estimating the Hurst exponent and its application. Results are discussed in Section 4. Finally, the conclusions are drawn in Section 5.
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1043
2. Mathematical preliminaries In this section, we will introduce some related terminologies and notations for the rest of the paper. Besides, Levinson algorithm and common model order selection methods will be briefly introduced. The notation {x(t ), t ∈ R} is used to denote a continuous-time random process, and {x[n ], n ∈ Z } a discrete-time process. 2.1. FBM and its variants The FBM is a generalized Brownian motion, first officially proposed by Mandelbrot and Van Ness [29] and denoted as the random function BH (t , ω ) , where ω belongs to a sample space Ω , simplified as BH (t ) . It can be defined as follows: BH (0 ) = b0 , BH (t2 ) − BH (t1 ) =
1 Γ (H + 1 2 )
{³
t2
−∞
}
(t2 − s )H −1 2 dB(s ) − ³−∞ (t1 − s )H −1 2 dB(s ) , t1
(1)
where H is the Hurst exponent, the only parameter with a value lying between 0 and 1. It is easy to understand that FBM will become an ordinary Brownian motion at H = 0.5 . Unfortunately, FBM is not a stationary process, i.e., its spectrum is time-dependent. For obtaining a time-independent spectrum, its increment of FBM, denoted as BH′ (t ) and called FGN, is considered. FGN is a stationary process and its spectrum is expressed as follows: S BH′ (Ω ) =
1 Ω
2 H −1
.
(2)
Furthermore, the discrete version of FBM, called the DFBM, is denoted as BH [n] = BH (nTs ) , where Ts denotes the sampling time. The DFGN is denoted as x H [n] = BH [n] − BH [n − 1] and its corresponding ACF as rH [k ] =
(
)
VH 2H 2H 2H , k +1 − 2 k + k −1 2
(3)
where VH = Γ(1 − 2 H ) cos(Hπ ) Hπ [1, 6]. The ACF rH [k ] behaves asymptotically as k 2 H −2 = k −α , α ∈ (0, 2 ) [10]. Chang and Chang [19] have proved the DFGN is a regular process using comparisons between two types of ACFs of DFGN and fdGn. Because of the regular property, an AR model can be used to describe the behavior of DFGN. Furthermore, Levinson algorithm is utilized to efficiently obtain the coefficients of an AR model. Some of its characteristics have been summarized by Lundahl et al. [5] Moreover, Samorodnitsky and Taqqu [10] have a complete introduction to these topics. 2.2. Levinson algorithm The Levinson algorithm [31], or Levinson-Durbin recursion [32], is an efficient method for solving the Yule-Walker equation or the Wiener-Hopf equation. Its computational cost can further be reduced by using the Schur algorithm [32], in which inner products are no longer required. This algorithm can provide the AR parameters for all lower-order AR models fitting to the ACFs. This property is useful, especially when users do not know a prior the correct model order. If the AR model is finite and free
1044
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
of noise, one can use this algorithm to derive the correct model order and its corresponding parameters according to the minimum variation of the prediction error power. However, if the AR model order is infinite or the model is finite but contaminated by noise, the algorithm alone might not choose the correct order. For later notation need and implementation, the Levinson algorithm [31] is summarized as follows. The Levinson algorithm recursively computes the parameter sets {a1 [1], ρ1} , {a2 [1], a2 [2], ρ 2 } , …, {a p [1], a p [2],! , a p [ p ], ρ p }. The final set at order p is the desired solution of the Yule-Walker equa-
tions. If x[n] is an AR of order p , then a p [i ] = a [i ] for i = 1, 2,! , p and ρ p = σ 2 , where σ 2 is the excitation noise variance. The algorithm is first initialized by a1 [1] = −
rxx [1] rxx [0]
(4)
and
ρ1 = (1 − a12 [1]) rxx [0] ,
(5)
where rxx is the ACF of x[n] , and then the following parameter sets are computed by k −1
ak [k ] = −
rxx [k ] + ¦ ak −1 [l ]rxx [k − l ] l =1
ρ k −1
,
(6)
ak [i ] = ak −1 [i ] + ak [k ]ak −1 [k − i ] , i = 1, 2, ! , k − 1 ,
(7)
ρ k = (1 − ak2 [k ]) ρ k −1 ,
(8)
and
for k = 2, 3, !, p . 2.3. Model order selection Six common methods for order selection of the AR model include the final prediction error (FPE) [31–32, 37–38], Akaike information criterion (AIC) [31–32, 37–38], criterion autoregressive transfer function (CAT) [31, 37–38], minimum description length (MDL) criterion [32, 37–38], residual variance (RV) [37–38], as well as Hannan and Quinn (HQ) [37–38]. The FPE estimates the model order as the place minimizing the following equation: FPE (k ) =
N +k ρˆ k , N −k
(9)
where ρˆ k is the estimate of the white noise variance (prediction error power) [31], or called the average prediction-error power [32], for the kth order AR model. The second method is the AIC, which is defined as the following equation: AIC(k ) = N ln ρˆ k + 2k ,
(10)
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1045
As before, minimizing the AIC gives the order selected. The performance of the AIC is similar to the one of the FPE. For larger data records, the two criteria will produce almost identical model order results because they are mathematically related to each other. The third method is the MDL and is defined as MDL(k ) = N ln ρˆ k + k ln N .
(11)
The MDL is asymptotically consistent unlike the AIC, which is not a consistent estimator. Similarly, the model order estimate is obtained by minimizing the MDL. The fourth method, the CAT, is defined as CAT(k ) =
1 N
k
1
¦ ρ~ i =1
i
1 − ~ , ρ
(12)
k
where
ρ~i =
N ρˆ i . N −i
(13)
The model order estimate can be chosen by minimizing the CAT. The fifth method, called the RV, is defined as RV (k ) =
N −k ρˆ k . N − 2k − 1
(14)
Similarly, the model order estimate is obtained by minimizing the RV. The sixth method, called the HQ, is defined as HQ(k ) = N ln ρˆ k + 2k ln ln N .
(15)
The model order is estimated as the position of the minimum of the HQ. These methods work well on finite-order situations. However, according to Wold and Kolmogorov theorems, the AR model order for a regular process is infinite. Therefore, we would like to know whether or not they are applicable to an infinite order model. 3. Procedures of estimating the Hurst exponent and its application It has been shown that the DFGN is a regular process [19], and hence the transformation of the AR model into the MA model can be omitted. Experimental results have shown that the AR method is superior to the MA method in accuracy and efficiency. Since the AR model of the DFGN is an infinite order, the higher the model order is, the closer to a regular process it is. It can be seen from the true ACFs, where the true ACF rH [k ] behaves asymptotically as k 2 H −2 = k −α , α ∈ (0, 2 ) , never a zero for any lag [10]. Accordingly, the order of the AR model for estimated ACFs was chosen as the maximum possible order by intuition, i.e., the data length minus one [19]. However, for finite data, errors easily happen at higher orders because the AR model highly depends on inaccurate estimates of the ACFs. These inaccurate estimates further result in adverse effects on the estimation of the Hurst exponent. Therefore, a trade-off between theoretical demand and practical estimation will be essential. In this paper, the standard or referenced curve of the slopes of power spectral density (PSD) versus
1046
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
the Hurst exponent was obtained using the true ACFs, hence the higher the model order is, the closer to a true model it is. Unfortunately, an infinite order cannot be modeled in a finite time-space, and estimation is a very complex process. In order to simplify the complexity of comparison, a 1000 order of AR model was adopted in the paper and a fitting curve of 4 degrees was derived. The estimated slopes were obtained by the estimated ACFs based on the model order chosen from six model order selection methods. Finally, an estimated Hurst exponent was produced through an established standard curve ( H versus slope) beforehand, described as follows. First, every true Hurst exponent was used to produce an estimated slope of AR power spectrum estimates versus frequencies in a logarithmically scaled way sˆ = f (H ) , H ∈ {0.01, 0.02, ! , 0.99} ,
(16)
and then the standard curve was established using least squares polynomial curve-fitting [39] of 4 degrees as follows: Hˆ = c 4 sˆ 4 + c3 sˆ 3 + c2 sˆ 2 + c1 sˆ1 + c0 sˆ 0 ,
(17)
c0 = 0.4998 , c1 = -0.4327 , c2 = 0.0400 , c3 = 0.0192 , and c4 = 0.0040 .
(18)
where
The detailed steps are listed in the following procedure: Procedure 1: Production of one standard curve for N = 1000 1. 2. 3. 4.
Begin the Hurst exponent with H = 0.01 with increment 0.01. Calculate the true ACFs from (3). Estimate the AR parameters and the white noise variance (prediction error power) using these true ACFs via Levinson algorithm or from (4) – (8). Calculate Aˆ ( f ) = 1 + aˆ [1]exp {− j 2πf } + " + aˆ [ pˆ ]exp {− j 2πfpˆ } ,
5.
(19)
where pˆ is the estimated order of the AR model. For true ACFs, the order of the AR model is selected as N − 1 . Estimate the PSD of the AR model via the following equation 2 Sˆ ( f ) = ρˆ pˆ | Aˆ ( f ) | ,
(20)
where ρˆ pˆ is the estimate of the white noise variance (prediction error power) for order pˆ . 6.
7. 8.
{ (
)}
Estimate the slope, sˆ , of log Sˆ ( f i ) versus {log( f i )} over FM = { f1 , f 2 , ! , f M } , where f1 = 0.5 M , f i − f i −1 = 0.5 M for i = 2, 3, !, M with M being selected as 63 in this paper. If H = 0.99 , then go to next; or execute the next H . Use the least squares polynomial curve-fitting of 4 degrees to fit these 99 points of the Hurst exponent versus the estimated slope.
Once a standard curve established, a procedure for estimating the Hurst exponent and its performance evaluation is provided to determine whether or not the model order selection method can con-
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1047
tribute to the accuracy of estimation. For the purpose, six model order selection methods were chosen, including the FPE, AIC, MDL, CAT, RV, and HQ; six types of data length, including 128, 256, 512, 1024, 2048, and 4096; 13 Hurst exponents, including 0.01, 0.05, 0.1, 0.2, ! , 0.9 , 0.95, and 0.99. Then, 100 realizations of DFGN, x H [n ] , for different data length and Hurst exponent were simulated according to the algorithm of Lundahl et al. [5] The algorithm provides a fine correlation structure and longterm dependency. The detailed steps are listed in the following procedure: Procedure 2: Estimation of Hurst exponent and its performance evaluation
1. 2. 3. 4. 5. 6.
Begin the model order selection method with the FPE. Begin the data length with N = 128 . Begin the Hurst exponent with H = 0.01 . Input the first realization of x H [n ] for N and H . Modify each realization of x H [n ] into zero mean. Estimate ACFs for k = 0,1, ! , N − 1 using the following form rˆ[ k ] =
1 N
N
¦ x[n]x[n − k ] .
(21)
n = k +1
Estimate the AR parameter sets using these estimated ACFs via Levinson algorithm or from (4) – (8), and then memorize these parameters. 8. Use the model order selection method to choose the optimal order, pˆ . 9. Extract the corresponding AR parameter set with order pˆ from memory. 10. Estimate the PSDs of AR model using (19) and (20). 11. Calculate the slopes of log Sˆ ( f i ) versus {log( f i )} over FM = { f1 , f 2 , ! , f M } , where f1 = 0.5 M , f i − f i −1 = 0.5 M for i = 2, 3, ! , M with M being selected as 63. 12. Estimate the H through the established standard curve (17). 13. If the number of realization is 100, then go to next; or execute the next realization. 14. Compute the mean mean-squared error between the estimated H and true H over 100 realizations. 15. If H = 0.99 , then go to next; or execute the next H . 16. Compute the mean mean-squared error over these 13 Hurst exponents. 17. If N = 4096 , then go to next; or execute the next N . 18. If the model order selection method is equal to the HQ, then Stop; or execute the next method. 7.
{ (
)}
From Procedure 2, one relatively better model order selection method will be discovered. Users can use the recommended method of order selection to more effectively estimate the Hurst exponent in the signals of interest. For easy use, a complete procedure is described as follows: Procedure 3: Use of the proposed estimator to a practical signal
1. 2. 3.
Input a section of practical signal, x[n] . Modify x[n] into zero mean. Estimate ACFs for k = 0,1, ! , N − 1 using (21).
1048
4.
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
8.
Estimate the AR parameter sets using these estimated ACFs via Levinson algorithm or from (4)–(8), and then memorize these parameters. Use the model order selection method, for example, the RV, to choose the optimal order, pˆ . Extract the corresponding AR parameter set with order pˆ from memory. Estimate the PSDs of AR model using (19) and (20). Calculate the slopes of log Sˆ ( f ) versus {log( f )} over F = { f , f , ! , f } , where
9.
f1 = 0.5 M , f i − f i −1 = 0.5 M for i = 2, 3, ! , M with M being selected as 63. Estimate the H through the established standard curve (17).
5. 6. 7.
{ (
i
)}
i
M
1
2
M
4. Results and discussion
In order to verify the effectiveness of the estimator proposed in this paper, a Hurst exponent estimator based on AR model with order selection, two sets of data were simulated, each producing 100 realizations of DFGN. The purpose of another data set is to confirm that experimental results are not affected by different realizations. In addition, the original AR method without model order selection, simply called the AR method, was also implemented to compare. Other six model order selection methods considered into the AR model will be simply referred to as the FPE, AIC, MDL, CAT, RV, and HQ, not the AR-FPE, AR-AIC, etc. Their experimental results are listed in Tables 1–6. Table1 provides the mean mean-squared error for each data length over 13 Hurst exponents and 100 realizations on Set 1; Table 2, the mean values of orders selected for each data length over 13 Hurst exponents and 100 realizations on Set 1; Table 3, the standard deviation values of orders selected for each data length over 13 Hurst exponents and 100 realizations on Set 1. Tables 4–6 are the results of Set 2 corresponding to Tables 1–3. It is obvious from Tables 1 and 4 that the AR methods with considering order selection are more accurate than the original one without considering order selection. Moreover, the RV method is almost the best except for N = 4096 . The FPE, AIC, and CAT are almost equivalent to one another. This phenomenon is consistent with the description that the performance of the CAT has been found to be similar to that of the FPE and AIC [31]. The MDL is the worst among the six methods. As described by Haykin [32], the application of the FPE criterion has been found to produce excellent results in the case of a pure AR process. However, when the FPE criterion is applied to practical non-AR data, the criterion appears to be conservative in the sense that it often yields too low a value for the model order. As with the FPE criterion, the AIC has been found to select too low a value for the model order when it is applied to practical non-AR data. Experimental results also confirm this argument. The FPE and AIC indeed choose a lower order than the RV. Table 1 Comparison of methods for mean mean-squared error on Set 1 Methods
N = 128
N = 256
N = 512
N = 1024 N = 2048 N = 4096
AR FPE AIC MDL CAT RV HQ
3.37E-03 3.54E-03 3.54E-03 7.14E-03 3.55E-03 2.80E-03 4.89E-03
3.06E-03 1.38E-03 1.38E-03 3.27E-03 1.39E-03 1.21E-03 2.18E-03
3.41E-03 8.72E-04 8.72E-04 2.13E-03 8.74E-04 6.55E-04 1.33E-03
3.29E-03 4.64E-04 4.64E-04 1.09E-03 4.65E-04 3.71E-04 6.83E-04
3.01E-03 2.29E-04 2.29E-04 6.53E-04 2.30E-04 1.83E-04 3.52E-04
3.00E-03 1.36E-04 1.36E-04 4.15E-04 1.36E-04 1.48E-04 2.35E-04
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1049
Table 2 Comparison of methods for mean values of orders selected on Set 1 Methods FPE AIC MDL CAT RV HQ
N = 128 3.02 3.03 1.38 2.98 6.42 1.90
N = 256 4.51 4.52 1.87 4.43 10.74 2.70
N = 512 6.64 6.64 2.69 6.54 17.28 3.81
N = 1024 9.36 9.36 3.78 9.26 29.72 5.48
N = 2048 12.57 12.57 5.00 12.51 55.73 7.33
N = 4096 17.69 17.69 6.93 17.67 74.40 10.01
Table 3 Comparison of methods for standard deviation values of orders selected on Set 1 Methods FPE AIC MDL CAT RV HQ
N = 128 2.22 2.22 0.66 2.17 4.65 1.00
N = 256 2.60 2.61 0.78 2.53 7.13 1.34
N = 512 3.49 3.49 1.00 3.38 10.51 1.50
N = 1024 4.39 4.39 1.11 4.30 20.36 1.90
N = 2048 4.98 4.98 1.15 4.95 35.31 2.06
N = 4096 6.08 6.08 1.47 6.08 45.65 2.50
On the other hand, the AIC is statistically inconsistent in the sense that the probability of selecting the correct order does not approach unity as the data length N approaches infinity. In fact, there is a tendency for the AIC to overestimate the model order in the case of a pure AR process as N approaches infinity. Although the MDL was developed to overcome the asymptotically inconsistent behavior of the AIC, its purpose cannot work on this infinite order of AR model. Tables 2 and 5 show that the MDL excessively underestimates the model order in the case of an infinite AR process. Table 4 Comparison of methods for mean mean-squared error on Set 2 Methods
N = 128
N = 256
N = 512
AR FPE AIC MDL CAT RV HQ
4.61E-03 3.25E-03 3.12E-03 4.22E-03 1.88E-03 9.90E-04 4.21E-03 1.88E-03 9.90E-04 7.59E-03 4.09E-03 2.19E-03 4.20E-03 1.88E-03 9.90E-04 3.54E-03 1.65E-03 8.36E-04 5.32E-03 2.77E-03 1.35E-03
N = 1024 N = 2048 N = 4096 3.18E-03 3.05E-03 4.24E-04 2.39E-04 4.24E-04 2.38E-04 1.01E-03 6.47E-04 4.25E-04 2.40E-04 3.46E-04 1.71E-04 6.11E-04 3.79E-04
2.86E-03 1.38E-04 1.38E-04 4.15E-04 1.38E-04 1.59E-04 2.42E-04
Table 5 Comparison of methods for mean values of orders selected on Set 2 Methods FPE AIC MDL CAT RV HQ
N = 128
N = 256
N = 512
3.33 3.34 1.38 3.26 6.35 1.88
4.67 4.69 1.89 4.55 11.11 2.69
6.44 6.44 2.72 6.41 17.01 3.85
N = 1024 N = 2048 N = 4096 9.09 9.10 3.71 9.00 28.89 5.31
12.53 12.54 5.06 12.47 48.17 7.34
18.76 18.76 6.81 18.64 89.84 10.13
1050
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
Table 6 Comparison of methods for standard deviation values of orders selected on Set 2 Methods FPE AIC MDL CAT RV HQ
N = 128 2.98 2.98 0.64 2.83 4.55 1.01
N = 256 3.07 3.12 0.82 2.79 7.64 1.27
N = 512 3.17 3.17 1.05 3.14 10.85 1.60
N = 1024 N = 2048 N = 4096 4.34 5.12 7.12 4.35 5.13 7.12 1.08 1.20 1.49 4.25 5.06 7.05 19.56 32.60 55.81 1.84 2.05 2.58
5. Conclusions
In theory, an infinite order of autoregressive (AR) model can be used to model the discrete-time fractional Gaussian noise and its increment process, the discrete-time fractional Brownian motion (DFBM). Therefore, in a signal of DFGN or DFBM with finite data length, the order of an AR model should be chosen to be as large as possible by intuition. Unfortunately, the estimation of Hurst exponent using an AR model highly depends on the estimation of autocorrelation functions (ACFs), which easily leads to larger errors in larger lags. In order to resolve the dilemma between a theoretical requirement and an estimation difficulty, six common order selection methods were introduced to choose a suitable order for the AR model. Experimental results show that the estimation of the Hurst exponent with considering order selection into the AR model outperforms the original one without considering order selection. Moreover, the residual variance (RV) method is almost the best among these six methods except in the data length of 4096. 6. Acknowledgements
This work was supported by the National Science Council, Republic of China, under Grant NSC 101-2221-E-040 -010. References [1]
B.B. Mandelbrot and J.W.V. Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (1968), 422-437. [2] B.B. Mandelbrot, Fractals: Form, Chance and Dimension, W. H. Freeman, New York, 1977. [3] B.B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman, New York, 1983. [4] A.P. Pentland, Fractal-Based Description of Natural Scenes, IEEE Trans. Pattern Anal. Machine Intell. PAMI-6 (1984), 661-674. [5] T. Lundahl, J. Ohley, S.M. Kay, and R. Siffert, Fractional Brownian motion: A maximum likelihood estimator and its application to image texture, IEEE Trans. Med. Image, MI-5 (1986), 152-161. [6] P. Flandrin, On the Spectrum of Fractional Brownian Motions, IEEE Trans. Inform. Theory 35 (1989), 197-199. [7] K. Falconer, Fractal Geometry: Mathematical Foundations and Applications. John Wiley & Sons, New York, 1990. [8] N. Sarkar and B.B. Chaudhuri, An efficient approach to estimate fractal dimension of textural images, Pattern Recognit. 25 (1992), 1035-1041. [9] S.S. Chen, J.M. Keller, and R.M. Crownover, On the calculation of fractal features from images, IEEE Trans. Pattern Anal. Mach. Intell. 15 (1993), 1087-1090. [10] G. Samorodnitsky and M.S. Taqqu, Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance, Chapman & Hall, New York, 1994.
Y.-C. Chang et al. / A Hurst exponent estimator based on autoregressive power spectrum estimation
1051
[11] N. Sarkar and B.B. Chaudhuri, An efficient differential box-counting approach to compute fractal dimension of image, IEEE Trans. Syst. Man Cybern. 24 (1994), 115-120. [12] X.C. Jin, S.H. Ong, and Jayasooriah, A practical method for estimating fractal dimension, Pattern Recognit. Lett. 16 (1995), 457-464. [13] M.S. Taqqu, V. Teverovsky, and W. Willinger, Estimators for long-range dependence: An empirical study, Fractals 3 (1995), 785-798. [14] S.-C. Liu and S. Chang, Dimension estimation of discrete-time fractional Brownian motion with applications to image texture classification, IEEE Trans. Image Process. 6 (1997), 1176-1184. [15] S. Chang, S.-T. Mao, S.-J. Hu, W.-C. Lin, and C.-L. Cheng, Studies of detrusor-sphincter synergia and dyssynergia during micturition in rats via fractional Brownian motion, IEEE Trans. Biomed. Eng. 47 (2000), 1066-1073. [16] R. Jennane, W.J. Ohley, S. Majumdar, and G. Lemineur, Fractal analysis of bone X-ray tomographic microscopy projections, IEEE Trans. Med. Imaging 20 (2001), 443-449. [17] A. Peréz, C.E. D’Attellis, M. Rapacioli, G.A. Hirchoren, and V. Flores, Analyzing blood cell concentration as a stochastic process, IEEE Eng. Med. Biol. Mag. 20 (2001), 170-175. [18] A.N. Esgiar, R.N.G. Naguib, B.S. Sharif, M.K. Bennett, and A. Murray, Fractal analysis in the detection of colonic cancer images, IEEE Trans. Inf. Technol. Biomed. 6 (2002), 54-58. [19] Y.-C. Chang and S. Chang, A fast estimation algorithm on the Hurst parameter of discrete-time fractional Brownian motion, IEEE Trans. Signal Process. 50 (2002), 554-559. [20] S. Chang, S.-J. Hu, and W.-C. Lin, Fractal dynamics and synchronization of rhythms in urodynamics of female Wistar rats, J. Neurosci. Meth. 139 (2004), 271-279. [21] J.-M. Lee, U. Yoon, J.-J. Kim, I.Y. Kim, D.S. Lee, J.S. Kwon, and S.I. Kim, Analysis of the hemispheric asymmetry using fractal dimension of a skeletonized cerebral surface, IEEE Trans. Biomed. Eng. 51 (2004), 1494-1498. [22] S. Chang, S.-J. Li, M.-J. Chiang, S.-J. Hu M.-C. Hsyu, Fractal dimension estimation via spectral distribution function and its application to physiological signals, IEEE Trans. Biomed. Eng. 54 (2007), 1895-1898. [23] S. Chang, M.-C. Hsyu, H.-Y. Cheng, and S.-H. Hsieh, Synergic co-activation of muscles in elbow flexion via fractional Brownian motion, Chin. J. Physiol. 51 (2008), 376-386. [24] S. Chang, M.-C. Hsyu, H.-Y. Cheng, S.-H. Hsieh, and C.-C. Lin, Synergic co-activation in forearm pronation, Ann. Biomed. Eng. 36 (2008), 2002-2018. [25] S. Chang, Physiological rhythms, dynamical diseases and acupuncture, Chin. J. Physiol. 53 (2010), 77-90. [26] S. Chang, Fractional Brownian motion in biomedical signal processing, physiology, and modern physics, The 5th ICBBE, May 10-12, 2011; Wuhan, China. [27] Y.-C. Chang, L.-H. Chen, L.-C. Lai, and C.-M. Chang, An efficient variance estimator for the Hurst exponent of discrete-time fractional Gaussian noise, IEICE Trans. Fundam. Electron. Commun. Comput. Sci. E95-A (2012), 15061511. [28] J.L. Doob, Stochastic Processes, John Wiley & Sons, New York, 1953. [29] M.B. Priestley, Spectral Analysis and Time Series, Volume 2: Multivariate Series, Prediction and Control, Academic Press, New York, 1981. [30] A.N. Shiryaev, Probability, 2nd ed., translated by R. P. Boas, Springer-Verlag, New York, 1996. [31] S.M. Kay, Modern Spectral Estimation: Theory & Application, Prentice-Hall, New Jersey, 1988. [32] S. Haykin, Modern Filters, Macmillan, New York, 1989. [33] A. Papoulis, Predictable processes and Wold’s decomposition: A review, IEEE Trans. Acoust. Speech, Signal Process. 33 (1985), 933-938. [34] A. Papoulis, Levinson's algorithm, Wold’s decomposition, and spectral estimation, SIAM Rev. 27 (1985), 405-441. [35] M. Deriche and A.H. Tewfik, Maximum likelihood estimation of the parameters of discrete fractionally differenced Gaussian noise process, IEEE Trans. Signal Process. 41 (1993), 2977-2989. [36] K.L. Chung, A Course in Probability Theory, 2nd ed., Academic, New York, 1973. [37] R. Palaniappan, Towards optimal model order selection for autoregressive spectral analysis of mental tasks using genetic algorithm, International Journal of Computer Science and Network Security 6 (2006), 153-162. [38] A.M. Aibinu, A.R. Najeeb, M.J.E. Salami, and A.A. Shafie, Optimal model order selection for transient error autoregressive moving average (TERA) MRI reconstruction method, World Academy of Science, Engineering and Technology 42 (2008), 161-165. [39] D. Hanselman and B. Littlefield, Mastering MATLAB 6: A Comprehensive Tutorial and Reference, Prentice-Hall, New Jersey, 2001. [40] C.W.J. Granger and R. Joyeux, An introduction to long-memory time series models and fractional differencing, J. Time Ser. Anal. 1 (1980), 15-29.
Copyright of Bio-Medical Materials & Engineering is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
