Chinese Science Bulletin © 2007

SCIENCE IN CHINA PRESS

Springer

Information flow among neural networks with Bayesian estimation LI Yan, LI XiaoLi†, OUYANG GaoXiang & GUAN XinPing Centre for Networking Control and Bioinformatics, Institute of Electrical Engineering, Yanshan University, Qinhuangdao 066004, China

Estimating the interaction among neural networks is an interesting issue in neuroscience. Some methods have been proposed to estimate the coupling strength among neural networks; however, few estimations of the coupling direction (information flow) among neural networks have been attempted. It is known that Bayesian estimator is based on a priori knowledge and a probability of event occurrence. In this paper, a new method is proposed to estimate coupling directions among neural networks with conditional mutual information that is estimated by Bayesian estimation. First, this method is applied to analyze the simulated EEG series generated by a nonlinear lumped-parameter model. In comparison with the conditional mutual information with Shannon entropy, it is found that this method is more successful in estimating the coupling direction, and is insensitive to the length of EEG series. Therefore, this method is suitable to analyze a short time series in practice. Second, we demonstrate how this method can be applied to the analysis of human intracranial epileptic electroencephalogram (EEG) recordings, and to indicate the coupling directions among neural networks. Therefore, this method helps to elucidate the epileptic focus localization. phase synchronization, coupling direction, conditional mutual information, Bayesian estimation, epileptic EEG

The synchronization analysis in the dynamical system has been a focus of attention in various disciplines[1]. Synchronization phenomena have been found in physical systems and biological systems. A typical example is the cardio respiratory interaction[2]. The synchronization between brain areas is associated with human motion[3], sleep[4], learning[5] and consciousness[6]. However, excessive synchronization is likely to result in epileptic seizures[7]. Therefore, it is necessary to develop a method to detect the synchronization occurrences, in particular the coupling direction between brain areas. Rosenblum and his colleagues[1,2] put forward two methods to identify the coupling direction: evolution map approach (EMA) and instantaneous period approach (IPA). The two methods detect a driving or casual relationship through the phase dynamics between two channel signals. EMA and IPA have been successfully utilized to detect the coupling direction between www.scichina.com

www.springerlink.com

the heart and lungs. As the heart rate and respiration signals are regular (periodical) and less noisy, the EMA and IPA can obtain clear results. It is noted that since the EMA is sensitive to noises, it is not suitable for analyzing noisy and nonstationary EEG recordings. In refs. [8,9], state-space and phase-dynamics approaches have been proposed to detect the weak coupling direction. The same drawbacks as those of EPA and IPA exist. A directionality index based on conditional mutual information is proposed and applied to the instantaneous phases of weakly coupled oscillators[10], with which the coupling direction between oscillators can be identified. This method is abbreviated to IM in this paper. The advantage of IM is that it can reveal and quantify the posReceived November 12, 2006; accepted March 2, 2007 doi: 10.1007/s11434-007-0272-3 † Corresponding author (email: [email protected]) Supported by the National Natural Science Foundation of China (Grant No. 60575012)

Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

1 The methods 1.1 Estimation of signal phase In order to analyze the phase synchronization between two series, the first step is to estimate the signal phase. Three typical methods are addressed in ref. [12]. In this paper, the signal phase is estimated by Hilbert transform:

ς (t ) = s (t ) + isH (t ) = A(t )e

iφ ( t )

,

(1)

where sH(t) is a Hilbert transform of the signal s(t); φ (t) and A(t) are the phase and amplitude of the signal s(t), respectively. Hilbert transform does not need to set any parameters, so it is widely applied to estimate the phase of sequence series. The limitation of Hilbert transform is that it is only suitable for the analysis of a narrow frequency-band signal. To analyze a wide frequency-band signal, a decomposition method may be employed to gain several narrow frequency-band signals. The narrow frequency-band signal phases are then estimated by a Hilbert transform. In this paper, our first step is to use a Bior wavelet transform to decompose EEG data[13] into several narrow frequency-band signals, and then extract 8―16 Hz frequency-band EEG signals[14]. 1.2 Directionality index Before estimating the directionality index, we need to calculate the phase synchronization strength. The reason is that once two signals possess the phase synchronization, the directionality index may possibly exist. The phase synchronization is calculated by

Synab =

1 N

N −1

∑ exp(i(φa ( jΔt ) − φb ( jΔt ))) , j =1

(2)

− H (X1 , X 2 |X 3 ),

ARTICLES

where N is the length of the random discrete variable, 1/Δt is sampling frequency, and the value of Syn ranges from 0 to 1. The value of Syn is equal to 0, meaning that there is no phase synchronization between two signals, and the larger the value of Syn, the stronger the phase synchronization. Palus et al.[15,16] proposed a method for estimating the directionality index with conditional mutual information based on Shannon entropy. The conditional mutual information is defined as I (X 1 ∩ X 2 |X 3 ) = H (X1|X 3 ) + H (X 2 |X 3 )

(3)

where X1, X2, X3 are random variables, Xj = {xj(ti)}, ti = iΔt, i = 1, … NX, and Δt is the sampling time. Considering the two processes X and Y, their phases φx(t) and φy(t) are estimated by Hilbert transform. The conditional mutual information can be calculated through the joint probability distribution function (PDF), joint entropy and conditional entropy: (4) I xy = I (φx (t ) ∩ Δτ φ y (t + τ ) | φ y (t )), where Ixy is the conditional mutual information, and τ is the delay. The directionality index is calculated by I xy − I yx (5) IM ( x, y ) = . I xy + I yx The value of IM(x, y) ranges from −1 to 1. IM(x, y) > 0 means that the progress X(t) drives Y(t), IM(x, y) < 0 means that the progress Y(t) drives X(t), and IM(x, y) = 0 means that the interactions between X(t) and Y(t) are symmetrical. 1.3 Directionality index with Bayesian estimation

In the above method, the estimation of mutual information requires a long series. If the series is too short or includes a great deal of noises, the directionality index is inaccurate. In eq. (3), the joint conditional mutual information is calculated through two one-dimension conditional entropies and a two-dimension joint entropy; the two-dimension conditional entropy is calculated by a two-dimension joint entropy and a three-dimension joint entropy. In ref. [11], Bayesian entropy was adopted to replace Shannon entropy. In this paper, a two-dimension conditional entropy is decomposed into a two-dimension joint entropy and a three-dimension joint entropy, and the two-dimension and three-dimension joint entropy with Bayesian is calculated. (For the details, see Appendix in electronic version online). The directionality in-

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

2007

INFORMATION SCIENCE AND SYSTEM SCIENCE

sible asymmetry in the couplings. Unfortunately, the IM estimator is sensitive to the length of series, i.e. the bias and variance of IM decreases with increasing series length. The simulation results show that the IM can reliably detect the coupling direction between two chaotic oscillators once the length of data is larger than 104. In this paper, we propose another method with Bayesian estimation. The Bayesian estimation is based on the least square error[11]. The test results of the simulated and real EEG data show that Bayesian estimation can gain accurate entropy from the noisy and nonstationary data; therefore a more reliable directionality index can be estimated to describe the information flow among neural networks.

dex with Bayesian conditional mutual information is given by Iˆxy − Iˆyx . IM bys ( x, y ) = (6) Iˆxy + Iˆyx A drawback of Bayesian entropy is that it requires a long computation time. In order to overcome this drawback, we applied a Taylor series method to estimate Bayesian entropy. Therefore the directionality index estimation with Bayesian is much faster.

2 Simulation 2.1 LPCC model

A nonlinear lumped-parameter cerebral cortex (LPCC) ― model[17 20] is used to test the performance of this new method. In this model there are three important parameters to be considered: excitatory neuron parameter A, inhibitory neuron parameter B and connective strength K. The parameters A and B modulate the balance of excitation and inhibition (he(t) = u(t)Acte−at and h(t) = u(t)Bbte−bt), and the parameter K manages the degree and coupling direction between neural networks (NNs)[21]. When neurons are fired up, a static nonlinear sigmoid function[18] transforms the average membrane potential into the average pulse density of potentials, and then transfers this to other neurons along the axon of the firing neuron. In this study, three coupled NNs (NN1, NN2 and NN3) are built to generate three surrogated EEG signals. The details of the LPCC can be seen in ref. [20]. Three parameters in the NN1 model are set at A = 3.25 and B = 22. The model does not generate the epileptic signals; after 80 s, the parameter A of the NN1 increases to 3.585, and the NN1 generates an epileptic EEG series. After 150 s, parameter A of the NN1 reduces to 3.25, and the seizures disappear. Finally parameter A of the NN1 increases to 3.585 again and the NN1 generates an epileptic EEG series again, as shown in Figure 1(a). Different values of parameter K12 are set, the connection strength between NN1 and NN2, and the K23 between NN2 and NN3 is set as a Gaussian distributed random coefficient (mean = 130, variance = 50). The connection strength K12 = 200 ensures to transfer the NN1 excitation to the NN2. The EEG signal of the NN2 is plotted in Figure 1(b). As K23 is a value generated by a Gaussian distributed function, the NN2 drives the NN3 randomly. The generated EEG data in the NN3 is shown in Figure 1(c). 2008

Figure 1 The simulation of the LPCC models. (a) The EEG signals of the NN1; (b) the EEG signals of the NN2, driven by the NN1; (c) the EEG signals of the NN3, which is randomly driven by the NN2.

2.2 Comparison of IM index and IMbys index between simulated EEG

Before estimating the directionality index, a delay time τ should be set up. In this study, τ is set at 100 point. A moving-window method is used to estimate the phase synchronization among NNs over time. The effect of the data length on the directionality index is also analyzed. The simulated EEG in Figure 1 is analyzed using a moving-window of the length of 39 s (10k) with an overlap of 50%. As can be seen in Figure 2(a), the IM method can detect the coupling direction between the

Figure 2 The directionality index and phase synchronization with a moving window of 39 s (5k). (a) The directionality index between NNs with IM method; (b) the directionality index between NNs with Bayesian; (c) phase synchronization index. The solid lines represent the directionality index and phase synchronization between the NN1 and NN2; the dash line represents the directionality index and phase synchronization between the NN1 and NN3; the dash-dot line represents the directionality index and phase synchronization between the NN2 and NN3.

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

Figure 3 The directionality index and phase synchronization with a moving window of 19 s (5k). (a) The directionality index between NNs with the IM method; (b) the directionality index with Bayesian; (c) phase synchronization index.

Further shorting the moving-window as the length of 4 s (1k), as shown in Figure 4(a), the IM method failed to detect the coupling direction among NNs. However, Figure 4(b) shows that the IMbys method is still able to detect the coupling direction among NNs.

3 The information flow of brain areas during epileptic seizures 3.1 The EEG recordings

The Epilepsy Centre at Freiburg University in Germany provided the EEG data for analysis in this study. The EEG data was recorded in the intracranial depth elec-

ARTICLES Figure 4 The directionality index and phase synchronization with a moving window of 4 s (1k). (a) The directionality index between NNs with the IM method; (b) the directionality index with Bayesian; (c) phase synchronization index.

trodes implanted in a patient and inputted into a computer using a 16-bit analog-to-digital converter and sampled at 256 Hz. Three channels situated in the left hippocampus were analyzed. The EEG signals with an epileptic seizure are plotted in Figure 5.

Figure 5 The EEG recordings with intracranial depth electrodes implanted in a patient situated in the left hippocampus. The shadow represents ictal states.

3.2 Comparison of IM index and IMbys index between EEG recordings

A discrete wavelet transform is used to decompose the EEG recordings; the signals at the Alpha band (8―16 Hz) are extracted for further analysis. The phase of signals at the Alpha band is estimated with Hilbert transform. Then, the phase synchronization interrelationships among EEG series are described by IM index and IMbys index. The 390 s EEG recordings, including interictal and ictal states, are analyzed (Figure 5). A moving window of 10 s with an overlap of 50% is applied in this

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

2009

INFORMATION SCIENCE AND SYSTEM SCIENCE

NN1 and the NN2. However, it is unable to detect the coupling direction between the NN2 and the NN3; actually there is no interaction at the duration from 150 to 250 s, but the directionality index is set at zero. Figure 2(b) shows that the IMbys can successfully detect the directionality index. In particular, in the interictal state, the IMbys is close to zero. As can be seen in Figure 2(c), the phase synchronization index decreases to 0.8 during the seizures; this finding is the same as that in refs. [22― 25]. The moving-window is shorted to 19 s (5k), as can be seen in Figure 3(a). The directionality index IM is insensitive, and the directionality index between the NN1 and the NN2 (strong coupling) can still be detected. The IM method cannot detect the directionality index of weakly coupled NNs. Figure 3(b) shows that the IMbys method can indicate the coupling direction. The entire directionality index is greater than zero in the ictal state and close to zero in the interictal state.

study. Figure 6(c) shows that the phase synchronization strength index between the first and second EEG series significantly decreases in the ictal state (shadow). In terms of the phase synchronization theory, this process is called desynchronization. The phase synchronization index increases to a normal level after seizures. As can be seen in Figures 6―8(b), in the interictal state the interaction among EEG series is small. In Figures 6―8(a), the IM index oscillates around zero, unlike IMbys index at the interictal state, which can approximate to zero. In the ictal state, IM index is unable to indicate the coupling direction among the EEG series; but the IMbys index can significantly distinguish the interaction difference among EEG series at the different states. Figures 6(b) and 8(b) show that the interaction between the first and second EEG series is stronger than that between the first and third EEG series in the ictal state. The interaction between the second and third EEG series is the

Figure 6 The phase synchronization interrelationship between the first and second EEG series in Figure 5. (a) IM index; (b) IMbys index; (c) phase synchronization index.

Figure 7 The phase synchronization interrelationship between the second and third EEG series in Figure 5. (a) IM index; (b) IMbys index; (c) phase synchronization index. 2010

Figure 8 The phase synchronization interrelationship between the first and third EEG series in Figure 5. (a) IM index; (b) IMbys index; (c) phase synchronization index.

strongest in the ictal state (Figure 7). Figures 6―8 show that the IM method fails to distinguish the interaction differences among NNs at different states, but IMbys can. Several papers analyzed the topological structure of multi-channel signals. In ref. [2], an EMA method is proposed to analyze the topological structure of three unidirectional coupling Van Der Pol chaotic oscillators. In ref. [26], a conditional Granger causality method (based on the frequency decomposition) is proposed to analyze the multi-channel of a three-dimension ARMA model. In ref. [27], partial directed coherence (PDC) is applied to analyze the topological structure of a 5-dimension ARMA model, four coupling Van Der Pol chaotic oscillators and essential tremor. This paper concentrates on the phase information of multi-channel series to explore the coupling relationship among NNs. In terms of the findings in Figures 6―8, a topological structure among EEG series in different states is plotted in Figure 9. The EEG series of 100 s is adopted to con-

Figure 9 The topological structure of three EEG series in the interictal and ictal states. (a) The directionality index and topological structure in the interictal state; (b) the directionality index and topological structure in the ictal state.

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

4 Conclusion This paper proposed a new method to estimate the directionality index among NNs with a Bayesian estimator, avoiding the effect of data length on the directionality index. The directionality index estimated by the method proposed by Palus et al.[10] needs long data. The simulation studies show that when a long moving-window is 1 2 3

4 5

6

7

8

9

10

11

12

13 14

Rosenblum G M, Pikovsky A. Detecting direction of coupling in interacting oscillators. Phys Rev E, 2001, 64: 045202 Rosenblum G M. Identification of coupling direction: Application to cardio respiratory interaction. Phys Rev E, 2002, 65: 041909 Roelfsema P R, Engel A K, Konig P, et al. Visuomotor integration is associated with zero time-lag synchronization among cortical areas. Nature, 1997, 385: 157―161 Steriade M, McCormick D A. Thalamocortical oscillations in the sleeping and aroused brain. Science, 1993, 262: 679―685 Miltner W H, Braun C, Arnold M, et al. Coherence of gamma-band EEG activity as a basis for associative learning. Nature, 1999, 397: 434―436 Rodriguez E, George N, Lachaux J P, et al. Perception’s shadow: Long-distance synchronization of human brain activity. Nature, 1999, 397: 430―433 Niedermeyer E, Fernando L S. Electroencephalography: Basic Principles, Clinical Applications and Related Fields. Baltimore: Winlliams & Wilkins, 1993. 1097 Smirnov D A, Andrzejak R G. Detection of weak directional coupling: Phase-dynamics approach versus state-space approach. Phys Rev E, 2005, 71: 036207 Smirnov D A, Bezruchko B P. Estimation of interaction strength and direction from short and noisy time series. Phys Rev E, 2003, 68: 046209 Palus M, Stefanovska A. Direction of coupling from phases of interacting oscillators: An information-theoretic approach. Phys Rev E, 2003, 67: 055201 Yokota Y. An approximate method for Bayesian entropy estimation for a discrete random variable. In: Proceedings of the 26th Annual International Conference of the IEEE EMBS, San Francisco, USA, 2004. 99―102 Rosenblum G M, Pikovsky A. Phase Synchronization: From Theory to Data Analysis. Handbook of Biological Physics. Vol 4, Neuro-informatics. Elsevier Science, 2001. 279―321 Li X L, Yao X, John R G, et al. Interaction dynamics of neuronal oscillations with wavelet. J Neurosci Meth, 2007, 160(1): 178―185 Li X L. Temporal structure of neuronal population oscillations with empirical model decomposition. Phys Lett A, 2006, 356(3):

ARTICLES

applied, the two methods can both successfully detect the coupling direction among NNs. However, once the moving-window length is shortened, the IM method fails to indicate the coupling direction among NNs. The IM method is sensitive to the data length; however, the improved method is insensitive to the data length during the estimation of directionality index. The hippocampus EEG recordings of a patient are also used to compare the performance of the two methods. The results show that IM method can hardly gain a reliable estimation in the coupling direction because of the effect of data length. However, the improved method can successfully detect the coupling direction from the short EEG series. The topological structure based on the directionality index could help to understand the information flow in brain areas. This new method can enable us to check the information flow among NNs during cognitive processes as well. We thank the Epilepsy Centre at the Freiburg University in Germany for providing the EEG data.

15

16 17 18

19

20

21

22

23 24

25

26

27

237―241 Palus M, Komarek V, Hrncir Z, et al. Synchronization as adjustment of information rates: Detection from bivariate time series. Phys Rev E, 2001, 63: 046211 Cover T M, Thomas J A. Elements of Information Theory. New York: John Wiley & Sons, 1991 Lopes da Silva F H, Hoek A, Smith H, et al. Model of brain rythmic activity. Kybernetic, 1974, 15: 27―37 Jansen B, Rit V G. Electroencephalogram and visual evoked potential generation in a mathermatical model of coupled cortical columns. Biol Cybern, 1995, 73(4): 357―366 Jansen B H, Zouridakis G, Brandt M E. A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol Cybern, 1993, 68(3): 275―283 Wendling F, Bellanger J J, Bartolomei F, et al. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biol Cybern, 2000, 83(4): 367―378 Li Y, Li X L, Ouyang G X, et al. Strength and direction of phase synchronization of neural networks. In: International Symposium on Neural Networks, Chongqing, China, 2005. 314―319 Mormann F, Andrzejak R G, Kreuz T. Automated detection of a preseizure state based on a decrease in synchronization in intracranial electroencephalogram recordings from epilepsy patients. Phys Rev E, 2003, 67: 021912 Li X L, Ouyang G X. Nonlinear similarity analysis for epileptic seizures prediction. Nonlinear Anal-Theor, 2006, 64(8): 1666―1678 Li X L, Kapiris P G, Polygiannakis J, et al. Fractal spectral analysis of pre-epileptic seizures phase: In terms of criticality. J Neural Eng, 2005, 2: 11―15 Li X L, Ouyang G X, Yao X, et al. Dynamical characteristics of pre-epileptic seizures in rats with recurrence quantification analysis. Phys Lett A, 2004, 333: 164―171 Chen Y H, Bressler S L, Ding M Z. Frequency decomposition of conditional Granger causality and application to multivariate neural field potential data. J Neurosci Meth, 2006, 150: 228―237 Schelter B, Winterhalder M, Eichler M, et al. Testing for directed influences among neural signals using partial directed coherence. J Neurosci Meth, 2005, 152: 210―219

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

2011

INFORMATION SCIENCE AND SYSTEM SCIENCE

struct the topological structure in the interictal state. As shown in Figure 9(a), no significant coupling relationship is found among the EEG series. However, a significant coupling relation is found from the three EEG recordings of 280 s in the ictal state. The first EEG series drives the second and third EEG series, and therefore the site of the first EEG series is a potential focus of the seizures. The epileptic seizure either starts from the first site, spreading to the second before progressing to the third; or travels directly from the first site to the third site.

Appendix In this Appendix, two-dimension Bayesian entropy is described. Three-dimension entropy can be obtained by the extension of the two-dimension conditional entropy. A conditional entropy of a binary random variable can be defined by H ( X | Y ) = ent (r | s ) = ent (rx, y ) − ent ( s y ), (A1) where X, Y are two random variables, r is two dimensions joint PDF, and s is one dimension marginal PDF. The entropy estimation can be obtained through the Bayesian estimation of the occurrence probabilities r and s. Bayesian conditional entropy is defined below: H bys ( X | Y ) = E (ent (r | s ))

The Dirichlet distribution is derived from a multinomial distribution, so we may write (n1,1 !) (n p ,q !) n n , ∫ r1,11,1 rp,pq,q dr = p q ( pq − 1 + ∑∑ ni , j )! i =1 j =1

∫

n sq q ds

s1n1

(A2)

where Er |n [•] means the expectation of • with res|n

spect to the posterior probability density distribution p(r|n) and p(s|n). The posterior probability density distribution p(r|n) and p(s|n) can be calculated by p( n | r ) p (r ) p ( n | r ) p (r ) p ( r | n) = = , p ( n) p ( | ) p ( )d n r r r ∫ (A3) p (n | s ) p ( s ) p(n | s ) p( s) p ( s | n) = = . p ( n) ∫ p(n | s) p( s)ds If a priori knowledge about occurrence probabilities r and s is unknown, the prior probability density distribution p(r) of occurrence probabilities r and s may be considered as a constant[11]: ⎧⎪ p q ⎫⎪ p (r ) = V , ∀r ∈ Rr ≡ ⎨r | ∑∑ r (i, j ) = 1, r (i, j ) ≥ 0 ⎬ , ⎪⎩ i =1 j =1 ⎪⎭ ⎧⎪ q ⎫⎪ p ( s ) = C , ∀s ∈ Rs ≡ ⎨ s | ∑ s ( j ) = 1, s ( j ) ≥ 0 ⎬ . ⎩⎪ j =1 ⎭⎪ Combining eq. (A3) to (A2), it is p(n | r ) p ( r | n) = , p ( s | n) = ∫ p(n | r )dr

p(n | s )

∫ p(n | s)ds

,

(A4)

(A5)

and p (n | r ) = PN ,r (n)

=

N! n r1,11,1 (n1,1 !) (n p ,q !)

p (n | s) = PN , s (n) = 2

q

(q − 1 + ∑ n j )!

(A7)

.

In terms of eqs. (A6), (A7) and

∑i=1 ni = N , ∑i=1 ni = N , p

q

the denominator of eq. (A3) is expressed by N! n n ∫ p(n | r )dr = (n1,1 !) (n p,q !) ∫ r1,11,1 rp,pq,q dr N! , ( pq − 1 + N )!

=

N! n ∫ p(n | s)ds = ∫ (n1 !) (nq !) r1 1

=

(A8) n rq q ds

N! . (q − 1 + N )!

Therefore, applying eqs. (A5) and (A7) to eq. (A3), the posterior PDF p(r|n) and marginal PDF p(s|n) are represented by ( pq − 1 + N )! p ( r | n) = p(n | r ) N! ( pq − 1 + N )! n1,1 n = r1,1 rp ,pq,q , (n1,1 !) (n p ,q !) p ( s | n) =

(q − 1 + N )! n1 r1 (n1 !) (nq !)

(A9)

n

rq q .

Substituting eq. (A8) to (A2), it yields ⎡ ( pq − 1 + N )! n1,1 Hˆ bys ( X | Y ) = ∫ ∫ ⎢ent (r ) r1,1 (n1,1 !) (n p ,q !) ⎢⎣

n

rp ,pq,q

n ⎤ (A10) rq q ⎥ drds. ⎥⎦ As for large p and q, it is difficult in the calculation for the multiple integrals in eq. (A9). It is noted that eq. (A4) is constant, so a multi-Taylor series is applied to estimate the entropy function ent(r) and ent(s), the details are below:

−ent ( s )

p

(q − 1 + N )! n1 r1 (n1 !) (nq !)

q

ent (r ) = −∑∑ ri , j log(ri , j )

n

rp ,pq,q ,

N! r1n1 (n1 !) (nq !)

(nq !) j =1

r |n s|n

= ∫ ∫ [ent (r ) p(r | n) − ent ( s ) p( s | n)]drds,

=

(n1 !)

(A6) n

rq q

i =1 j =1

p q ⎡ = ent (r 0 ) − ∑∑ ⎢(1 + log(ri0, j )(ri , j − ri0, j )) i =1 j =1 ⎣

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

q

p

(A11)

ent ( s ) = ∑ s j log( s j ) = ent ( s 0 )

⎤ (−1)k ( s 0j ) − k +1 ( s 0j − s j ) k ) ⎥ , k = 2 k ( k − 1) ⎦ ∞

+∑

p

q 1 ∑∑ (1 + log(ri0,j ))(ni, j + 1) N + pq i =1 j =1

+

1 q (1 + log( s 0j ))(n j + 1) ∑ N + q j =1

k ⎛ k ⎞ ( N + q − 1)! 1 (−1)l ∑ ⎜ ⎟ l k k − N + q + l − ( 1) ( 1)! k =2 l =0 ⎝ ⎠ K

xy

q

(n j + j )!

j =1

ni , j !

k ⎛ k ⎞ ( N + pq − 1)! 1 (−1)l ∑ ⎜ ⎟ l k k N pq l − + + − ( 1) ( 1)! k =2 l =0 ⎝ ⎠

(ni , j + j )!

k ⎛ k ⎞ ( N + q − 1)! 1 +∑ (−1)l ∑ ⎜ ⎟ l k ( k − 1) ( N + q + l − 1)! k =2 l =0 ⎝ ⎠ K

j =1

ni , j !

,

(A13)

bys

x

τ y

y

− Hˆ bys (φx (t ) ∩ Δτ φ y (t + τ ) | φ y (t )),

(A14)

Iˆyx = Iˆbys (φ y (t ) ∩ Δτ φx (t + τ ) | φ x (t )) = Hˆ bys (φ y (t ) | φx (t )) + Hˆ bys (Δτ φ x (t + τ ) | φ x (t ))

ni , j !

(n j + j )!

,

= Hˆ bys (φx (t ) | φ y (t )) + Hˆ bys (Δτ φ y (t + τ ) | φ y (t ))

K

q

ni , j ,λ !

Applying eqs. (A12) and (A11) to eq. (3), we get the Bayesian conditional mutual information below: Iˆ = Iˆ (φ (t ) ∩ Δ φ (t + τ ) | φ (t ))

−∑

l =1 j =1

(ni , j ,λ + j )!

+∑

∑ (si0, j )− j +1

−

∑ (si0, j )− j +1

υ

l =1 j =1 λ =1

j =1

∑∑ (ri0, j )−i − j +1

q

∑∑∑ (ti0, j ,λ )−i − j −υ +1

q

q

1 q ∑ (1 + log(s0j ))(n j + 1) N + q j =1

s 0j = 1 / 2, j = 1 q , which

p

p

+

k ⎛ k ⎞ ( N + pqυ − 1)! 1 (−1)l ∑ ⎜ ⎟ l k = 2 k (k − 1) l = 0 ⎝ ⎠ ( N + pqυ + l − 1)!

Hˆ bys ( X | Y ) = ∑∑ ri0, j − ∑ s 0j i =1 j =1

υ 1 (1 + log(ti0, j ,λ ))(ni , j ,λ + 1) ∑∑∑ N + pqυ i =1 j =1 λ =1 q

K

will guarantee the entropy function is convergence. The proven details are in ref. [11]. Applying eq. (A10) to (A9), the entropy can be estimated by q

j =1

−

p

−∑

where r0 and s0 should satisfy r0 ≥ 1/2 and s0 ≥ 1/2,

p

q

i =1 j =1 λ =1

q ⎡ −∑ ⎢(1 + log( s 0 )( s j − s 0j ) j =1 ⎣

p, j = 1 q,

υ

q

Hˆ bys ( X ∩ Y | Z ) = ∑∑∑ ti0, j ,λ − ∑ s 0j

j =1

ri0,j = 1 / 2, i = 1

ARTICLES

ditional joint entropy is also calculated by

(A12)

where K=3. It is sufficient to meet the approximation precision and stabilization. In this paper, Bayesian con-

− Hˆ bys (φ y (t ) ∩ Δτ φ x (t + τ ) | φ x (t )).

Applying eq. (A13) to (4), the directionality index with the Bayesian conditional mutual information is calculated by Iˆxy − Iˆyx (A15) IM bys ( x, y ) = . Iˆ + Iˆ

LI Yan et al. Chinese Science Bulletin | July 2007 | vol. 52 | no. 14 | 2006-2011

xy

yx

3

INFORMATION SCIENCE AND SYSTEM SCIENCE

⎤ (−1) k (ri0, j )− k +1 (ri , j − ri0, j ) k ⎥ , k = 2 k ( k − 1) ⎦ ∞

+∑