Autoregressive model in the Lp norm space for EEG analysis.

G Model

ARTICLE IN PRESS

NSM 7078 1–9

Journal of Neuroscience Methods xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Journal of Neuroscience Methods journal homepage: www.elsevier.com/locate/jneumeth

1

Computational Neuroscience

2

Autoregressive model in the Lp norm space for EEG analysis

3 4 5 6 7 8 9

Q1

Peiyang Li a,1 , Xurui Wang a,1 , Fali Li a , Rui Zhang a , Teng Ma a , Yueheng Peng b , Xu Lei c , Yin Tian d , Daqing Guo a , Tiejun Liu a , Dezhong Yao a , Peng Xu a,∗ a Key Laboratory for NeuroInformation of the Ministry of Education, School of Life Science and Technology, University of Electronic Science and Technology of China, Chengdu, China b School of Microelectronics and Solid-State Electronics, University of Electronic Science and Technology of China, Chengdu, China c School of Physiology, Southwest University, Chongqing, China d School of Bioinformatics, Chongqing university of Posts and Telecommunications, Chongqing, China

10

11 12 13 14 15 16

h i g h l i g h t s • • • •

Designed an Lp (p ≤ 1) norm-based residual model to estimate autoregressive (AR) parameters. The Lp (p ≤ 1) norm AR model estimates parameters more robustly than an L2 norm-based AR mode for time series with outliers. The Lp (p ≤ 1) norm AR holds a lower relative error of AR parameters and higher prediction accuracy than L2 norm-based methods. A resting EEG power spectrum estimated by the Lp (p ≤ 1) norm AR model is less influenced by ocular artifacts compared with L2 norm-based AR.

17

18 31

a r t i c l e

i n f o

a b s t r a c t

19 20 21 22 23 24

Article history: Received 6 September 2014 Received in revised form 8 November 2014 Accepted 10 November 2014 Available online xxx

25

30

Keywords: Autoregressive model Lp norm EEG Power spectrum

32

1. Introduction

26 27 28 29

33 34 35 36 37 38 39 40

The autoregressive (AR) model is widely used in electroencephalogram (EEG) analyses such as waveform fitting, spectrum estimation, and system identification. In real applications, EEGs are inevitably contaminated with unexpected outlier artifacts, and this must be overcome. However, most of the current AR models are based on the L2 norm structure, which exaggerates the outlier effect due to the square property of the L2 norm. In this paper, a novel AR object function is constructed in the Lp (p ≤ 1) norm space with the aim to compress the outlier effects on EEG analysis, and a fast iteration procedure is developed to solve this new AR model. The quantitative evaluation using simulated EEGs with outliers proves that the proposed Lp (p ≤ 1) AR can estimate the AR parameters more robustly than the Yule–Walker, Burg and LS methods, under various simulated outlier conditions. The actual application to the resting EEG recording with ocular artifacts also demonstrates that Lp (p ≤ 1) AR can effectively address the outliers and recover a resting EEG power spectrum that is more consistent with its physiological basis. © 2014 Published by Elsevier B.V.

Power spectral density (PSD) and variable states are the two important measures for characterizing the physiological information underlying EEGs. There are two main approaches to estimating these measurements, nonparametric and parametric methods. Nonparametric methods (Xu et al., 2011), such as Fourier transforms and periodograms, use the observed data to directly perform the estimation. However, this type of approach is usually problematic because of leakage and frequency resolution in PSD estimates,

∗ Corresponding author. Tel.: +86 028 83206978. E-mail address: [email protected] (P. Xu). 1 Both the authors equally contribute to this work.

which require a large number of samples. In the parametric approach, a random signal is characterized by the parameters estimated from the finite record of the data, so there is no need to make assumption about how the data were generated. Currently, three types of parametric models, the autoregressive (AR) model, the moving average (MA) model and the autoregressive moving average (ARMA) model, are used for related estimations (Xu et al., 2011), among which the AR model is by far the most widely used due to the following merits. First, with a suitable order, it can approximate any stationary random process. Second, the AR model is suitable for representing spectra with narrow peaks. Third, the AR model has a set of very simple linear equations for parameter estimation so that many efficient algorithms are available (Jain and Dandapat, 2005). Compared with AR, both the MA and ARMA models, as a general rule, require more coefficients to represent the

http://dx.doi.org/10.1016/j.jneumeth.2014.11.007 0165-0270/© 2014 Published by Elsevier B.V.

Please cite this article in press as: Li P, et al. Autoregressive model in the Lp norm space for EEG analysis. J Neurosci Methods (2014), http://dx.doi.org/10.1016/j.jneumeth.2014.11.007

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

G Model

ARTICLE IN PRESS

NSM 7078 1–9

P. Li et al. / Journal of Neuroscience Methods xxx (2014) xxx–xxx

2

signal spectrum information. In various studies, AR has been proven to be able to build the more meaningful spectrum information of EEG than the FFT-based analysis (Güler et al., 2001). 58 In recent years, there have been a large number of applications 59 using AR-related models in various EEG-related studies such as EEG 60 system identification, EEG power spectrum estimation, EEG net61 work analysis, and brain computer interfaces (Chen et al., 2013; Qiu, 62 2011; Wang et al., 2014b). Wang et al. (2014b) proposed a corrected 63 approach to remove the ocular artifact using AR-based system iden64 tification, resulting in substantial performance improvement. In 65 Chen et al. (2013), AR is used to extract the phase and frequency 66 information from intracranial EEGs. Dahal et al. (2014) used a time 67 varying AR model to extract features to delineate the attention and 68 distraction in the motor driving experiment. Wang et al. (2010) 69 developed an improved feature extraction method based on the 70 multivariate adaptive autoregressive (MVAAR) model for the clas71 sification of motor imagery. 72 In real-world applications, EEGs will usually be contaminated 73 with outliers due to eye blinks or head movement, which will 74 greatly influence the AR estimation. Most of the previous work 75 in this area paid less attention to the outlier effect, even though 76 the adopted schemes such as moving averaging and sparse con77 straint may actually compress the outlier effect to some degree 78 (Songsiri, 2013). Theoretically, the previously used AR model and 79 its variants are based on the L2 norm structure, and the L2 norm 80 will exaggerate the outlier effect due to the square property of the 81 L2 norm (Blankertz et al., 2007; Songsiri, 2013). Compared with the 82 L2 norm, the Lp (p ≤ 1) norm has been proven to be robust to out83 liers, and it has been widely used in a diversity of signal processing 84 applications such as denoising, EEG inverse problem, MRI/CT recon85 struction, feature extraction (Chartrand, 2009; Li et al., 2013; Lustig 86 87Q2 et al., 2007; Xu et al., 2007). In consideration of the merits of the Lp (p ≤ 1) norm to compress the outlier effect, we will restructure 88 the residual equation for AR parameter estimation in the Lp (p ≤ 1) 89 norm space and establish a fast iteration procedure to solve this 90 new AR model. 91 56 57

In essence, the AR model minimizes the residual errors for all the observed samples:

2 q N N 2 2 E[un ] = un = wk x(i − k) x(i) − i=1

i=1

(4)

x(q)

⎢ x(q + 1) ⎢ A=⎢ .. ⎣ .

x(q − 1)

···

x(1)

x(q)

···

x(2)

.. .

x(N − 1)

93

94

2. Materials and methods 2.1. Autoregressive model The AR model is usually expressed as: q

95

x(n) = −

wk x(n − k) + u(n)

(1)

k=1

96 97 98 99 100

where u(n) is the input sequence to the system and is usually considered to be zero-mean white Gaussian noise with a variance 2 . x(n) is the observed data, representing the output sequence. w {wk , 1 ≤ k ≤ q} is the corresponding AR parameters with q being the order of the AR model. The system transfer function is given by

.. .

x(N − 2)

···

⎥ ⎥ ⎥ ⎦

(5)

H(z) =

1 B(z) = q C(z) 1+ w z −k k=1 k

103 104

105

2 Pxx (f ) = w |H(f )|2 =

(6)

By taking the derivative of (6) with respect to W under the condition df/dw = 0, we can obtain the follow formulation: T

T

2A AW − 2A Y = 0

(7)

and the objective parameters W can be estimated as W = (AT A)

−1 T

A Y

2 w

112

114

(8)

∗

W = arg minf (W ) = arg min||Y W N−q

p − AW ||p

|xq+i − A(i, :)W |p

117

119

120 121 122 123 124 125 126 127

129 130 131 132 133 134 135 136 137 138 139 140 141 142

143

(9)

144

i=1

where || • ||p denotes the Lp (p ≤ 1) norm of a vector; we refer to this model as the Lp-AR estimation. The gradient for this function is

116

128

In real-world applications, outliers will create an unexpected effect on related analyses such as spectrum estimation, signal prediction. To improve the robustness of the AR parameters estimation, some schemes such as sparse constraint with Lp (p ≤ 1) norm terms are proposed in various AR variant versions to alleviate the noise effect (Ping-bo and Zhi-ming, 2006). However, most of them mainly focus on the imposing restrictions on the parameters, leaving the main structure of the objective function in the L2 norm space. Unfortunately, the L2 norm object function will inevitably exaggerate the outlier effect no matter how the AR parameters are emphasized. We will define the AR object function in the Lp (p ≤ 1) norm space, aiming to improve the AR robustness to the outlier effect. The AR object function is defined in Lp (p ≤ 1) norm space as ∗

115

118

In addition to the least square algorithm, other approaches such as the Yule–Walker (Y-W) equations and the Burg method (Xu et al., 2011) are also used to estimate the AR parameters that can minimize the sum of the residual errors in (4). No matter what scheme is used to estimate the AR parameters, the inherent L2 norm structure in (4) and (6) indicates that the influence of outliers will be exaggerated due to the square property of the L2 norm, resulting in the biased AR parameters.

= arg min

q 2 1 + w e−2kf k=1 k

111

W

(2)

where the C(z) and B(z) represent the poles and zeros of the system response, respectively. Based on the transfer function, the AR based spectrum estimation at frequency f has the form

110

113

arg minf (W ) = ||Y − AW ||22

W 102

109

x(N − q)

Eq. (4) can be formatted as

W 101

108

⎤

2.2. Lp (p ≤ 1) norm based autoregressive model 92

107

k=1

Let W = [w1 , w2 , . . ., wq ]T ; Y = [x(q + 1), x(q + 2), . . ., x(N)]T , with N being the length of signal; || • ||2 denotes the L2 norm of a matrix or a vector; and A ∈ R(N−q)×q be the delay array:

⎡

106

145 146 147

n−q

(3)

g=p

|xq+i − A(i, :)W |p−1 sgn(i)(−AT (i, :))

(10)

i=1


148

G Model

ARTICLE IN PRESS

NSM 7078 1–9

P. Li et al. / Journal of Neuroscience Methods xxx (2014) xxx–xxx 149

where the polarity function sgn(i) is set as

150

151 152 153 154

sgn(i) =

xq+i − A(i, :)W > 0

1

−1 xq+i − A(i, :)W ≤ 0

⎧H ⎨ k ⎩H

k

155

, if Sk ≤ 0

+

ˇk Wk WkT

gk = gk+1 − gk , ˇk = 1 +

−

Hk gk WkT

Sk Wk = Wk+1 − Wk

−

Wk gkT Hk

, if Sk > 0

(12)

gkT Hk gk Sk

Sk = WkT gk 156 157 158

159 160 161 162 163 164 165 166 167

Here, gk indicates the change of the gradient and Wk is the displacement of W in the kth iteration. Using the gradient in equation (10) and the Hessen matrix in (12), the Lp-AR is solved by

(1) Setting the iteration number to k = 0, the termination error to ε ∈ [0, 1], initializing W as a random non-zero q-dimension vector, and setting the initial pseudo Hesse matrix H0 with unit matrix Iq ; (2) Computing the gradient gk following Eq. (10). If ||gk || < ε, W* = Wk and end the iteration; (3) Solving the linear coupled equations Hk d = −gk and obtaining the search direction dk ; (4) Finding the optimal learning velocity ˛k : ˛k = arg minf (Wk + ˛

168 169

170 171 172

173

174 175 176 177 178 179 180

˛dk ) and updating the AR parameter W as Wk+1 = Wk + ˛k dk ; (5) k = k + 1; using Eq. (12) to update Hk , go to step 2. With the above iteration, the optimal W can be estimated. In the current work, the Lp-AR with various p (p = 1.0, 0.8, 0.6, 0.4, 0.2) is adopted and the termination error ε is 1.0e−6.

3. Evaluation index In our current work, we mainly focus on the performance difference between the conventional AR parameter-estimation algorithms (e.g. Yule–Walker, Burg, LS) and Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR when the dataset is contaminated by outliers. One concern is the effect of outliers on AR parameters, and another concern is the effect on the fitting performance. The Euclidean distance is utilized to denote the performance as

181

1 ∧ d = || s − s||2 N

182

where s is the standard reference for AR parameters or time series, ∧

183 184 185 186 187 188 189 190

4. Results

191

4.1. Simulation studies

192

(11)

Based on above gradient, we adopted BFGS (Battiti and Masulli, 1990; Head and Zerner, 1985; Li and Fukushima, 2001) to find the optimal W for Lp-AR. In BFGS, the pseudo Hessen matrix is updated as follows:

Hk+1 =

3

(13)

s is the estimated counterpart of s, and N is the number of samples or the relevant AR parameters. For each experiment, the Akaike Information Criterion (AIC) was utilized to select the optimal order for AR methods within the range of 8–60 (Hardin, 2005). The AR variants, including Y-W, Burg, and least square (LS)-based AR, and our developed Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR were then used for the corresponding parameter estimation.

The simulations consist of two aspects: (1) evaluating the effect of outliers on the AR parameters and (2) evaluating the effect on the power spectrum estimation when outliers are introduced.

4.1.1. AR parameter estimation The P300 at electrode P3 of one subject is used for this evaluation. The EEGs are recorded using the EEG systems developed in our group. The P300 is evoked by the standard oddball paradigm. Each trial begins with the fixation point with 250 ms duration; 500 ms after the fixation point disappears, a stimulus, i.e., an upper (standard stimulus) or down (target stimulus) triangle will appear randomly on the screen center with a duration of 200 ms. The occurrence frequency ratio of the standard stimulus versus target stimulus is 4:1. After the stimulus disappears, the intertrial intervals (ITIs) range randomly between 1000 and 1200 ms. Subjects were instructed to press key “1” for the standard stimulus and key “2” for the target stimulus. Sixteen Ag/AgCl electrodes covering the sensorimotor area were used to record the EEG, and the signals were sampled at 1000 Hz. Each run consists of 250 trials, with 200 trials for the standard stimulus and 50 trials for the target stimulus. All experiments were performed in accordance with the Ethical Committee of University of Electronic Science and Technology of China (UESTC). The EEG data samples for our study were extracted from the interval between 0 ms and 800 ms (800 points) after the stimulus. The 50 trials for the target stimulus are averaged, resulting in 0.8 s long (800 points) P300 waveforms. The 0.8 s long P300 waveform is divided into two segments, with each being 0.4 s long. The first 0.4 s long segment is randomly corrupted with an outlier and the AR parameters are estimated from this segment. Based on the estimated AR parameters, the second segment is used to evaluate the fitting performance of the AR model. Using this P300 waveform, we quantitatively evaluated the possible effects of the outlier occurrence rate and outlier strengths on the model parameter estimation.

4.1.1.1. Effect of the outlier occurrence rate. To generate the datasets contaminated with different numbers of outliers, we added the outliers to the datasets while varying the number of outliers from 4 to 20 into the first 0.4-s-long P300 segment, where the time positions for the outliers’ corruption are randomly determined. The outlier is generated from the Gaussian distribution with fixed variance 2 = 1 and mean u = 52 ␮v. Taking the AR parameters estimated from the original P300 waveform as the reference, two bias errors, one for the AR parameters and another for the fitted waveform, are estimated. For each outlier number, the same procedure is repeated 100 times and the averaged bias error is estimated from the 100 runs. Fig. 1 gives an example illustration for the simulation in one of the 100 runs, where the first 0.4-s-long segment is corrupted with outliers and the second 0.4 s segment is used for curve fitting based on the AR parameters estimated from the first segment. The corresponding mean bias errors for the AR parameters and fitted waveforms are given in Tables 1 and 2, respectively. The paired t-test is used to investigate the difference for the eight AR variants, and the test reveals that the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR has statistically smaller bias errors than the other L2 AR methods for both the parameter estimation and waveform fitting. The bold values indicate the smallest bias errors among the eight AR variants, and “*, *, *, *, * ” indicate a significant difference (ps < 0.05) between the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR and the corresponding AR counterparts.


193 194 195

196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225

226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250

G Model

ARTICLE IN PRESS

NSM 7078 1–9


4

Original P300 waveform 100

Amplitude/uv

80 60 40 20 0 -20 -40

A 0

100

200

300

400

500

600

700

500

600

700

800

Time/ms P300 waveform with outliers 100

Amplitude/uv

80 60 40 20 0 -20 -40

B 0

100

200

300

400

800

Time/ms Fig. 1. P300 waveforms used for evaluation. Time series on the left side of the red dotted line are used for the AR parameters training; the data on the other side indicate the fitting result of AR model. (A) The fitting performance of the Lp (p ≤ 1) AR model without outlier influence. (B) The fitting performance of the Lp (p ≤ 1) AR model, when the training data mingle with the outlier. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1 Parameter estimation bias when the outlier number varied. Method

Y-W AR Burg AR LS AR L1 AR L0.8 AR L0.6 AR L0.4 AR L0.2 AR

251 252 253 254 255 256

257 258

Outlier number 4

8

12

16

20

1.22 ± 0.08***** 0.92 ± 0.06***** 0.92 ± 0.06*****

1.34 ± 0.05***** 0.96 ± 0.03***** 0.96 ± 0.03*****

1.38 ± 0.08***** 0.97 ± 0.05***** 0.97 ± 0.05*****

1.40 ± 0.08***** 0.97 ± 0.05***** 0.97 ± 0.05*****

1.42 ± 0.09*****

0.64 0.56 0.62 0.62 0.62

± ± ± ± ±

0.07 0.21 0.24 0.24 0.24

0.68 0.60 0.69 0.70 0.70

± ± ± ± ±

0.06 0.21 0.24 0.24 0.24

Specifically, the waveform fitting performance for the second 0.4-s-long waveform is shown in Fig. 2 when 8 outliers are introduced into the first segment of the P300 waveforms. Consistent with the fitting errors in Table 2, the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR model accurately reconstructs the waveforms based on the AR parameters learned in the first segment. 4.1.1.2. Effect of outlier strength. This simulation investigates the possible effect of outlier strengths on the AR variants. Similar to the

0.70 0.61 0.70 0.71 0.74

± ± ± ± ±

0.09 0.19 0.21 0.21 0.26

0.74 0.67 0.81 0.81 0.86

± ± ± ± ±

0.14 0.19 0.24 0.25 0.42

0.98 ± 0.05***** 0.98 0.75 0.67 0.81 0.82 0.96

± ± ± ± ± ±

0.05**** 0.13 0.19 0.23 0.24 0.64

above simulation, the outliers with different strengths are added to the first 0.4-s-long segment. The outlier number is 16, and the outlier is generated from the Gaussian distribution with fixed variance 2 = 1 and varied mean from 17 to 87. For each of the noise amplitudes, the procedure is repeated 100 times. The bias errors for the AR parameters and waveform fittings are listed in Tables 3 and 4. The paired t-test also demonstrates that the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR model consistently shows smaller bias errors for both the AR parameters and waveform fittings.

Table 2 Fitting errors when the outlier occurrence varied. Method

Outlier number 4

8

12

16

20

Y-W AR

1.22 ± 0.06*****

1.46 ± 0.08*****

1.59 ± 0.12*****

1.69 ± 0.13*****

1.79 ± 0.14*****

Burg AR

1.31 ± 0.08*****

1.55 ± 0.10*****

1.66 ± 0.16*****

1.74 ± 0.17*****

1.84 ± 0.17*****

LS AR L1 AR L0.8 AR L0.6 AR L0.4 AR L0.2 AR

1.31 ± 0.08*****

1.53 ± 0.10*****

1.64 ± 0.17*****

1.71 ± 0.18*****

1.80 0.75 0.98 1.21 1.21 1.32

0.25 0.60 0.70 0.71 0.71

± ± ± ± ±

0.07 0.20 0.24 0.25 0.27

0.40 0.72 0.92 0.91 0.93

± ± ± ± ±

0.16 0.21 0.24 0.24 0.35

0.55 0.84 1.05 1.06 1.10

± ± ± ± ±

0.11 0.21 0.27 0.28 0.34

0.65 0.94 1.17 1.18 1.21

± ± ± ± ±

0.19 0.24 0.32 0.32 0.46

± ± ± ± ± ±

0.18***** 0.49 0.24 0.36 0.37 0.73


259 260 261 262 263 264 265 266 267

G Model

ARTICLE IN PRESS

NSM 7078 1–9


Y-W AR

Burg AR

-10 -20

10 0 -10 -20

A 0

100

200

300

-30

400

0

100

L0.8 AR

300

-30

400

10 0 -10 -20

predicted value actual value

10 0 -10

200

300

400

-30

0 -10 -20

C 0

100

200

300

-30

400

100

Time/ms

200

300

400

200

300

400

L0.2 AR 30


20 10 0 -10

-30

100

Time/ms

-20

F 0

D 0

L0.4 AR

-20

E

10

30

20


20

Time/ms

Amplitude/mv


100

-10

L0.6 AR

Amplitude/mv

Amplitude/mv

200

30

0

0

Time/ms

30

-30

10

-20

B

Time/ms

20

20

Amplitude/mv

0

20

30 predicted value actual value

Amplitude/mv

10

L1 AR


Amplitude/mv

20

-30

LS AR


Amplitude/mv

Amplitude/mv

30

5


20 10 0 -10 -20

G 0

100

Time/ms

200

300

-30

400

H 0

Time/ms

100

200

300

400

Time/ms

Fig. 2. The predicted waveforms of eight AR model methods. The green lines depict the observed signal, corresponding to the right side of the red dotted line in Fig. 1. The red lines illustrate the predicted result of the eight AR models. (A) The predicted performance of the Y-W AR model when outliers are mixed in the original data. (B) The predicted performance of the Burg AR model. (C) The predicted performance of the LS AR model. (D)–(H) The predicted performance of the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3 Parameter estimation bias when the outlier strength varied. Method

17.39

34.78

52.17

69.56

86.95

Y-W AR

1.05 ± 0.07fx1

Burg AR

1.40 ± 0.08***** 0.98 ± 0.04**** 0.97 ± 0.04****

1.44 ± 0.08***** 0.99 ± 0.04**** 0.99 ± 0.04****

1.44 ± 0.09*****

0.97 ± 0.05***** 0.97 ± 0.05*****

1.31 ± 0.07***** 0.94 ± 0.05***** 0.94 ± 0.04*****


268 269 270 271 272 273 274 275

Noise amplitude

0.54 0.61 0.65 0.66 0.66

± ± ± ± ±

0.05 0.30 0.33 0.33 0.35

0.67 0.64 0.76 0.78 0.80

± ± ± ± ±

0.03 0.19 0.26 0.26 0.28

4.1.2. Power spectrum estimation One important application of the AR model is the power spectrum estimation. Here, we simulated a sinusoidal signal y = A1 cos (2w1 t) + A2 cos(2w2 t) + s(t), where A1 = 4.1, w1 = 8, A2 = 3.8, w1 = 20, and s(t) is a random signal with = 0, 2 = 0.5. Then, 8 outliers derived from the Gaussian distribution with 2 = 1 and u = 24 are added to investigate the influence of an outlier in the power spectrum estimation for various AR variants. The outlier

0.72 0.73 0.89 0.90 0.98

± ± ± ± ±

0.09 0.19 0.26 0.27 0.56

0.77 0.74 0.91 0.92 0.97

± ± ± ± ±

0.14 0.21 0.25 0.26 0.30

0.99 ± 0.04***** 0.99 0.84 0.77 0.88 0.90 0.94

± ± ± ± ± ±

0.04***** 0.17 0.21 0.22 0.23 0. 28

positions are randomly determined, and the procedure is repeated 200 times. For each run, the SNR at concerned frequency points (8 Hz and 20 Hz) is used to delineate the performance. SNR at frequency f is defined as

SNR(f ) =

20 log P(f )

f +5

i=f −5

(14)

P(i)

Table 4 Fitting errors when the outlier strength varied. Method

Noise amplitude 17.39

34.78

52.17

Y-W AR

0.97 ± 0.04*****

1.41 ± 0.09*****

1.72 ± 0.13*****

2.00 ± 0.12*****

2.30 ± 0.12*****

Burg AR

1.06 ± 0.06*****

1.51 ± 0.12*****

1.78 ± 0.16*****

2.02 ± 0.16*****

2.29 ± 0.14*****


1.06 0.42 0.85 0.96 0.95 0.96

± ± ± ± ± ±

0.06***** 0.09 0.19 0.32 0.32 0.35

1.50 0.52 0.91 1.24 1.26 1.29

± ± ± ± ± ±

0.12***** 0.05 0.25 0.36 0.38 0.41

1.75 0.61 1.08 1.48 1.49 1.50

± ± ± ± ± ±

69.56

0.17***** 0.16 0.33 0.45 0.45 0.51

1.97 0.82 1.15 1.51 1.52 1.54

± ± ± ± ± ±

86.95

0.17***** 0.81 0.36 0.38 0.39 0.47

2.21 1.79 1.25 1.58 1.60 1.62

± ± ± ± ± ±

0.15***** 2.33 0.37 0.47 0.49 0.57


276 277 278 279

280

G Model

ARTICLE IN PRESS

NSM 7078 1–9


6 Table 5

Q6 The SNR for the two frequencies. Methods

No outliers

With outliers

Burg AR

−7.39 ± 0.16***** −4.19 ± 0.14*****


−3.83 ± 0.07 −3.84 ± 0.36 −3.85 ± 0.50 −3.84 ± 0.56 −3.83 ± 0.60 −3.91 ± 0.80

−11.26 ± 0.92***** −7.15 ± 2.58***** −6.18 ± 0.58*****

SNR for 20 Hz (dB) Y-W AR

−4.85 ± 0.92 −4.82 ± 0.94 −5.19 ± 1.16 −5.77 ± 1.13 −5.98 ± 1.24

SNR for 8 Hz (dB)

281 282 283 284 285 286 287 288 289 290 291 292 293 294

Y-W AR

−1.37 ± 0.06*****

−3.67 ± 0.47*****

Burg AR

−0.69 ± 0.03*

−1.80 ± 0.37*****


−0.62 ± 0.02 −0.64 ± 0.18 −0.67 ± 0.26 −0.69 ± 0.29 −0.72 ± 0.35 −0.76 ± 0.47

−1.72 ± 0.29**** −0.93 ± 0.40 −0.90 ± 0.44 −1.18 ± 0.66 −1.48 ± 0.81 −1.64 ± 0.86

where P(f) is the amplitude of spectrum at frequency f. Obviously, a large SNR delineates a high ability for frequency identification. Fig. 3A shows the original data and Fig. 3B shows the outliercontaminated data in one simulation. The eight AR variants are then applied to the two datasets to estimate the power spectrum following Eq. (3). Fig. 4A–C are the corresponding power spectra of the L2-norm-based AR models (Yule–Walker, Burg, LS), and the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR-based power spectra are plotted in Fig. 4D–H. The averaged SNRs for the eight ARs are listed in Table 5 for the two concerned frequencies (8 and 20 Hz), where the paired ttest is used to investigate whether the difference between the L2 norm-based AR models and the Lp (p ≤ 1) AR model is significant (ps < 0.05). The bold values indicate the best SNRs among the eight

297

AR variants, and “*, *, *, *, * ” indicate the significant difference (ps < 0.05) between the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR model and the corresponding L2-AR counterparts.

298

4.2. Application to actual EEG recordings

295 296

299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321

The above simulation study reveals that the outliers have an obvious influence on the AR parameter estimation, resulting in a biased power spectrum estimation. In EEG recordings, ocular artifacts are inevitable, and they will unexpectedly distort the signals. We recorded a 2-min-long EEG signal using a BP system (Brain Products GmbH, Munich, Germany) for one male subject. During recording, we asked the subject to open his eyes and remain as relaxed as possible. Then, 64-channel EEGs are recorded, and we select the EEG at electrode F1 for further analysis. For the 2-min-long recording, the subject will occasionally blink his eye and ocular artifacts are introduced. We visually selected two 3-s-long segments from the whole recording, with one segment contaminated by ocular outliers and another segment free from ocular outliers. Based on the similar results for various Lp ARs revealed in the previous simulation study and the relative high efficiency of L1 AR, in this application of real data, we mainly focus on the performance difference between the conventional AR parameter-estimation algorithms (e.g., Yule–Walker, Burg, LS) and L1 AR. The four AR variants are applied to the two segments to estimate the corresponding power spectrum. Fig. 5 gives the two 3-s-long EEG segments used for the power spectrum, and Fig. 6 plots the corresponding PSD for the two segments when different AR variants are used.

As for the EEGs from the same subject, when a subject is kept relatively relaxed, the power spectrum will be kept relatively stable (Klimesch, 1999; Mantini et al., 2007; Qin et al., 2010). The plots in Fig. 6 reveal that the L2 norm-based AR variants are obviously influenced by the ocular outliers, resulting in different shapes of power spectra for the two segments. L1 AR has a relatively consistent power spectrum estimation due to the L1 norm merits of compressing outlier effect.

5. Discussion Outliers are inevitable in scalp EEGs. The simulation based on the P300 EEG dataset quantitatively evaluates the possible effects of outliers on AR parameter estimations, examining outlier occurrence rate and strength. The comparison reveals that all four AR models are influenced by outliers. As shown in Tables 1 and 2, when the occurrence number of outliers is increased from 4 to 20, the bias errors for both the AR coefficients and the waveform fitting are increased. Tables 3 and 4 further reveal that the stronger outliers will estimate more biased AR parameters. However, for all the tested outlier conditions, the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR model consistently shows a better performance than the other models with significantly smaller errors (ps < 0.05) for both the AR parameters and waveform fittings. When outliers are introduced into the dataset, outliers will be one of the main contributions to the residual errors in the AR object function, with their large amplitude. Obviously, the L2 norm will undesirably exaggerate the outliers because of its square properties. Therefore, L2 norm-based AR will provide a heavy emphasis on the outliers while the detailed information for the relevant signals are neglected. In essence, though the strategies used for AR parameters are different, Y-W, Burg and LS ARs are originally derived from the L2 norm space, which accounts for the poor performance when outliers are introduced into the signals. The Lp (p ≤ 1) AR model constructs the AR object function in Lp (p ≤ 1) norm space, and the Lp (p ≤ 1) norm will compress the outlier effect as shown in previous studies (Kwak, 2008). Following Eq. (3), the estimation of the power spectrum is based on AR parameters. Accordingly, the outlier effect on AR parameter estimation will be transferred to the power spectrum. The simulation study in Figs. 3 and 4 demonstrates the obvious effect of outliers on power spectrum estimation. As for the power spectrum, one important aspect is the discrimination of frequency components, i.e., the amplitude of the desired frequency component. The SNRs in Table 5 quantitatively reveal that without outliers, all eight ARs can recover the frequency information well with relatively higher SNRs and that no significant difference is observed for the eight ARs. When outliers are introduced, the SNRs for all eight ARs are dramatically lowered. However, compared with L2 norm-based approaches, Lp (p ≤ 1) AR can still hold better results. This consistent result is further revealed in Fig. 4, in which baseline drifts and compressed frequency peaks can be observed under the outlier conditions. The important aspect of Fig. 4 is that the Lp (p ≤ 1) AR model visually demonstrates superior ability to compress the outlier effect, where the very stable frequency peak information can still be recovered compared with the three other L2 norm-based ARs. Another aspect revealed in the simulation study is that the p values in Lp-AR actually influence the corresponding performance. In theory, a smaller p will facilitate compressing the outlier effect. However, when a p less than one is adopted, the gradient of the Lp object function will involve the rooting operation, which may dramatically increase the algorithm complexity. Moreover, when p is close to zero, the gradient may be prone to being disturbed by noise, resulting in the failure to estimate the optimal AR parameters. The different mathematical characteristics of the Lp norm account for


322 323 324 325 326 327 328 329

330

331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384

G Model

ARTICLE IN PRESS

NSM 7078 1–9


7

Original Sin waveform 40

Amplitude/mV

30 20 10 0

-10 -20

A 0

200

400

600

800

1000

1200

1400

1600

1800

1200

1400

1600

1800

2000

Time/ms Sin waveform with outliers 40

Amplitude/mV

30 20 10 0

-10 -20

B 0

200

400

600

800

1000

2000

Time/ms Fig. 3. Simulated sin signals for power spectrum estimation. (A) The original sinusoidal signal simulated by the equation in Section 4.1.2. (B) The waveform change when the original sinusoidal signal in (A) mingled with the outliers, where the spikes indicate the outlier influence.

Y-W AR

Burg AR

Orignal PSD PSD with Outliers


-10 -20 -30 -40

0 -10 -20 -30 -40

A 10

20

-50

30

0

L0.8 AR Orignal PSD PSD with Outliers

-20 -30 -40

30


-20 -30

0

10

20

20

Frequency/Hz

-30

-50

30

D 0

30

20

30

L0.2 AR 20



10

0 -10 -20 -30

-50

10

Frequency/Hz

-40

F 10

-20

-40

C

10

0

0

-10

L0.4 AR

-10

-50

0

Frequency/Hz

-40

E Frequency/Hz

-30

-50

30

Power/dB

0

20

-20

20

10

-10

10

-10

L0.6 AR

Power/dB

Power/dB

20

20

0

0

Frequency/Hz

20

-50

10


10

-40

B

Frequency/Hz

10


Power/dB

0

20

10

Power/dB

0

L1 AR

20

10

Power/dB

Power/dB

10

-50

LS AR

20

Power/dB

20

0 -10 -20 -30 -40

G 0

10

20

Frequency/Hz

30

-50

H 0

10

20

30

Frequency/Hz

Fig. 4. The power spectrum estimations of eight AR model methods. The blue solid line curve represents the original PSD estimated from the artifact-free data, and the red one indicates the PSD estimated from the outlier-contaminated data. (A) The power spectrum estimation performance of the Y-W AR model when the outliers are mixed in the original data. (B) The power spectrum estimation performance of the Burg AR model. (C) The power spectrum estimation performance of the LS AR model. (D)–(H) The power spectrum estimation performance of the Lp (p = 1.0, 0.8, 0.6, 0.4, 0.2) AR model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)


G Model

ARTICLE IN PRESS

NSM 7078 1–9


8

Resting EEG waveform 20 10

Amplitude/uv

0 -10 -20 -30 -40 -50

A 0

500

1000

1500

2000

2500

2000

2500

3000

Time/ms EEG waveform with ocular artifacts 20 10

Amplitude/uv

0 -10 -20 -30 -40 -50

B 0

500

1000

1500

3000

Time/ms Fig. 5. Two EEG segments used for power spectrum estimation. (A) The EEG signal without ocular outliers. (B) The EEG signal mingled with 3 ocular outliers.

385 386 387 388 389 390 391 392 393 394 395

the different Lp ARs performances revealed in the simulation study. Though a difference could be observed for various Lp ARs with p less than 1.0, the performance of L1 AR is very close to those Lp ARs (p = 0.8, 0.6, 0.4, 0.2), and no distinct difference exists among them. Therefore, in the practical applications, L1 version is usually considered due to its high efficiency and accepted performance. As for EEGs, ocular outliers will greatly influence the subsequent analysis (Lei et al., 2009; Wang et al., 2014b). However, during the EEG experiment, it is impossible to totally avoid ocular movement or blinks that will introduce very obvious outliers. The two segments in Fig. 5 are from the 2-min-long resting EEG recordings;

three ocular outliers are observed in the second segment. Fig. 6 demonstrates that the ocular outliers actually influence the power spectrum estimation as revealed in the previous simulation study. The resting EEG is usually characterized with the obvious spectrum peak in the alpha band (Tian and Yao, 2013). For the 3-s-long segment without an outlier, expect for the Y-W AR, the Burg, LS and L1 ARs revealed the alpha peak at approximately 10 Hz. When the EEG is contaminated with ocular outliers, the power spectra estimated with the Burg and LS AR models are very different from those estimated under ocular artifact-free conditions. More importantly, the alpha peaks are missed due to the ocular outlier effects when the

Y-W AR 50

PSD for pure EEG PSD for EEG with artifact

Power/dB

Power/dB

50

Burg AR

0

A

-50 0

5

10

15

20

25


0

B

-50 30

0

5

Frequency/Hz LS AR 50


0

C

-50 0

5

10

15

Frequency/Hz

15

20

25

30

Frequency/Hz L1 AR

Power/dB

Power/dB

50

10

20

25


0

D

-50 30

0

5

10

15

20

25

30

Frequency/Hz

Q5 Fig. 6. Different PSDs for the two EEG segments. The blue solid line curve represents the original PSD estimated from the artifact-free data, and the red one indicates the PSD estimated from the outlier-contaminated data. (A) The power spectrum estimation performance of the Y-W AR model when ocular outliers are mixed in the original data. (B) The power spectrum estimation performance of the Burg AR model. (C) The power spectrum estimation performance of the LS AR model. (D) The power spectrum estimation performance of the L1 AR model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)


396 397 398 399 400 401 402 403 404 405 406

G Model NSM 7078 1–9

ARTICLE IN PRESS P. Li et al. / Journal of Neuroscience Methods xxx (2014) xxx–xxx

407 408 409 410 411 412 413 414 415 416 417 418 419 420 421

L2 norm-based AR model is adopted. Similar to the results of the simulations, the L1 AR model shows the reliable ability to address the outliers in that the estimated spectrum from the contaminated condition is very close to that estimated from the uncontaminated EEG. Another merit in the spectrum of L1-AR is that the alpha peak is relatively clear for both two segments compared to the spectra estimated from the other L2-based ARs. The EEG is contaminated with various types of noise; even if the EEG does not have obvious ocular artifacts, other noise is undesirably introduced. Therefore, the L1 AR model can compress the noise effect and recover the strong alpha peak information as shown in Fig. 5. Various studies have proposed the L1 norm regularizationbased AR model to improve AR performance. Essentially, the L1 norm regularization-based AR model imposes the L1 norm constraints on projection w with the general form of arg min||Y − w

422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444

445

446 447 448 449 450 451 452 453 454 455 456

457

458 Q3 459 460 461

Aw||2 + ||w||1 (Hsu et al., 2008; Nardi and Rinaldo, 2011). It mainly focuses on imposing restrictions on the parameter w while leaving the main structure of the objective function defined in L2 norm space. In essence, the L1 norm regularization-based AR model aims to find the sparse solution of w, and the L2 norm-based residual object function will inevitably exaggerate the outlier effect though AR parameters that are emphasized with sparse constraints. Moreover, the determination of the regularization parameter is still an unsolved issue and limits its practical application in clinical situations. Compared with L1 norm regularization-based AR model, the proposed Lp-AR in this paper measures the residual errors in the Lp space, which will not exaggerate the outlier effect like the L2 norm and is a new structure for AR. In the future, we will develop the regularization AR in the Lp (p ≤ 1) residual space, which may provide more competitive ARs. Furthermore, in current neuroscience studies, brain network analysis using coherence as edge strengths attracts wide attention (Blankertz et al., 2008; Wang et al., 2014a; Xu et al., 2013, 2014; Zhang et al., 2013a,b). When coherence is estimated by the cross-spectrum, the effect of outliers on the power spectrum will be delivered to coherence. Therefore, the Lp (p ≤ 1) AR can be alternatively used for this issue to construct the more reliable EEG brain networks. 6. Conclusion In real-world applications with an AR model for EEG analysis, it is necessary to tackle the outlier effect. In essence, the current AR models are derived from the L2 norm space, which will be influenced by outlier noise. Compared to the conventional ARs, the Lp (p ≤ 1) norm AR we developed estimates the models in Lp (p ≤ 1) norm space, which can provide a more powerful ability to compress the effects of artifacts. The reported results in the current work illustrate that the Lp (p ≤ 1) norm can provide more meaningful merit for artifact compressing as proven in other reported studies (Blankertz et al., 2007; Kwak, 2008; Li et al., 2013; Oh and Kwak, 2013; Wang et al., 2013). Acknowledgements We acknowledge the support of the 973 Program 2011CB707803, the National Nature Science Foundation of China (#61175117, #81401484, #31200857 and #31100745), the Program for New Century Excellent Talents in University (#NCET-12-0089) and the 863 Project 2012AA011601.

9

References Battiti R, Masulli F. BFGS optimization for faster and automated supervised learning. In: International neural network conference; 1990. p. 757–60. Blankertz B, Dornhege G, Krauledat M, Müller K-R, Curio G. The non-invasive Berlin brain–computer interface: fast acquisition of effective performance in untrained subjects. NeuroImage 2007;37:539–50. Blankertz B, Tomioka R, Lemm S, Kawanabe M, Muller K-R. Optimizing spatial filters for robust EEG single-trial analysis. Signal Process Mag IEEE 2008;25: 41–56. Chen LL, Madhavan R, Rapoport BI, Anderson WS. Real-time brain oscillation detection and phase-locked stimulation using autoregressive spectral estimation and time-series forward prediction. Biomed Eng IEEE Trans 2013;60:753–62. Dahal N, Nandagopal DN, Cocks B, Vijayalakshmi R, Dasari N, Gaertner P. TVAR modeling of EEG to detect audio distraction during simulated driving. J Neural Eng 2014;11:036012. Güler I˙ , Kiymik MK, Akin M, Alkan A. AR spectral analysis of EEG signals by using maximum likelihood estimation. Comput Biol Med 2001;31:441–50. Hardin JW. Generalized estimating equations (GEE). Wiley Online Library; 2005. Head JD, Zerner MC. A Broyden–Fletcher–Goldfarb–Shanno optimization procedure for molecular geometries. Chem Phys Lett 1985;122:264–70. Hsu N-J, Hung H-L, Chang Y-M. Subset selection for vector autoregressive processes using lasso. Comput Stat Data Anal 2008;52:3645–57. Jain N, Dandapat S. Constrained autoregressive (CAR) model. In: INDICON, 2005 annual IEEE; 2005. p. 255–7. Klimesch W. EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res Rev 1999;29:169–95. Kwak N. Principal component analysis based on L1-norm maximization. Pattern Anal Mach Intell IEEE Trans 2008;30:1672–80. Lei X, Yang P, Yao D. An empirical Bayesian framework for brain–computer interfaces. Neural Syst Rehabil Eng IEEE Trans 2009;17:521–9. Li D-H, Fukushima M. A modified BFGS method and its global convergence in nonconvex minimization. J Comput Appl Math 2001;129:15–35. Li P, Xu P, Zhang R, Guo L, Yao D. L1 norm based common spatial patterns decomposition for scalp EEG BCI. Biomed Eng Online 2013:12. Mantini D, Perrucci MG, Del Gratta C, Romani GL, Corbetta M. Electrophysiological signatures of resting state networks in the human brain. Proc Natl Acad Sci U S A 2007;104:13170–5. Nardi Y, Rinaldo A. Autoregressive process modeling via the lasso procedure. J Multivar Anal 2011;102:528–49. Oh JH, Kwak N. Generalization of linear discriminant analysis using Lp-norm. Pattern Recognit Lett 2013;34:679–85. Ping-bo W, Zhi-ming C. LS-EM algorithm of parameters estimation for Gaussian mixture autoregressive model. In: IMACS multiconference on computational engineering in systems applications; 2006. p. 1–5. Qin Y, Xu P, Yao D. A comparative study of different references for EEG default mode network: the use of the infinity reference. Clin Neurophysiol 2010;121: 1981–91. Qiu H. A research to the autoregressive model-based efficient interpolation algorithm. In: 2011 Fourth international symposium on knowledge acquisition and modeling (KAM); 2011. p. 165–8. Songsiri J. Sparse autoregressive model estimation for learning Granger causality in time series. In: 2013 IEEE international conference on acoustics, speech and signal processing (ICASSP); 2013. p. 3198–202. Tian Y, Yao D. Why do we need to use a zero reference? Reference influences on the ERPs of audiovisual effects. Psychophysiology 2013;50:1282–90. Wang J, Xu G, Wang L, Zhang H. Feature extraction of brain–computer interface based on improved multivariate adaptive autoregressive models. In: 2010 3rd international conference on biomedical Engineering and Informatics (BMEI); 2010. p. 895–8. Wang H, Lu X, Hu Z, Zheng W. Fisher discriminant analysis with L1-norm. IEEE Trans Q4 Cybern 2013. Wang H, Li Y, Long J, Yu T, Gu Z. . An asynchronous wheelchair control by hybrid EEG–EOG brain–computer interface. Cogn Neurodyn 2014a:1–11. Wang Z, Xu P, Liu T, Tian Y, Lei X, Yao D. Robust removal of ocular artifacts by combining independent component analysis and system identification. Biomed Signal Process Control 2014b;10:250–9. Xu P, Yang P, Lei X, Yao D. An enhanced probabilistic LDA for multi-class brain computer interface. PLoS ONE 2011;6:e14634. Xu P, Tian C, Zhang Y, Jing W, Wang Z, Liu T, et al. Cortical network properties revealed by SSVEP in anesthetized rats. Sci Rep 2013:3. Xu P, Xiong X, Xue Q, Li P, Zhang R, Wang Z, et al. Differentiating between psychogenic nonepileptic seizures and epilepsy based on common spatial pattern of weighted EEG resting networks. IEEE Trans Biomed Eng 2014;61:1747–55. Zhang Y, Xu P, Guo D, Yao D. Prediction of SSVEP-based BCI performance by the resting-state EEG network. J Neural Eng 2013a;10:066017. Zhang Y, Xu P, Huang Y, Cheng K, Yao D. SSVEP response is related to functional brain network topology entrained by the flickering stimulus. PLOS ONE 2013b;8:e72654.


462

463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539

Generalized 2-D Principal Component Analysis by Lp-Norm for Image Analysis.

Lp-norm regularization in volumetric imaging of cardiac current sources.

A Hybridized GA with LP-LP Model for the Management of Confined Groundwater.

No-reference image sharpness assessment in autoregressive parameter space.

Adaptive Autoregressive Model for Reduction of Noise in SPECT.

The Role of Familiarity for Representations in Norm-Based Face Space.

EEG Source Reconstruction.

Is First-Order Vector Autoregressive Model Optimal for fMRI Data?

Bounds for the norm of lower triangular matrices on the Cesàro weighted sequence space.

Functional connectivity analysis in EEG source space: The choice of method.

An Lp (0 ≤ p ≤ 1)-norm regularized image reconstruction scheme for breast DOT with non-negative-constraint.

Pre-Surgical fMRI Data Analysis Using a Spatially Adaptive Conditionally Autoregressive Model.

Calibration of optical tweezers based on an autoregressive model.

Mixture of autoregressive modeling orders and its implication on single trial EEG classification.

Autoregressive model for high-resolution wavenumber estimation in a shallow water environment using a broadband source.

Long-range dependencies make the difference-Comment on "A stochastic model for EEG microstate sequence analysis".

Application of an autoregressive integrated moving average model for predicting injury mortality in Xiamen, China.

Learning in the model space for cognitive fault diagnosis.

A model to predict the thermal reaction norm for the embryo growth rate from field data.

Neutral space analysis for a Boolean network model of the fission yeast cell cycle network.

Model space diabatization for quantum photochemistry.

Kinematic analysis following implantation of the PRESTIGE LP.

Hölder norm estimate for the fractal Hilbert transform in Douglis analysis.

Individual level covariate adjusted conditional autoregressive (indiCAR) model for disease mapping.