System for Filtering Events Occurred Within a Predetermined Period of Time

The present invention relates to a system for filtering events which occur within a predetermined period of time.

There exists no system for filtering a large number of events in an event log database, which occur due to a single cause, with considering time sequential causality.

It is therefore an object of the present invention to provide a system for filtering events occurring within a predetermined period, whereby simultaneously generated mutual causalities with a large number of parameters can be effectively filtered.

According to the present invention, a system for filtering events occurring within a predetermined period of time is provided. The system includes a first generation unit for generating association rules between the events (parameters) which simultaneously occur within the predetermined period of time, a second generation unit for generating groups of events suited for the VAR (Vector Auto Regressive) model analysis based upon the association rules generated in the first generation unit, an analysis unit for analyzing the groups of events generated in the second generation unit, the analysis of the event groups being executed by using the VAR model analysis, and a filtering unit for verifying reliability of the association rules generated in the first generation unit, the verification being executed by using the VAR model.

Thus, any cross causalities of a large number of events (parameters), which causalities have not been predicted beforehand and simultaneously occurring can be effectively discovered.

It is preferred that the first generation unit includes a unit for generating association rules between the events based upon a basket analysis method using support considering Poisson error and confidence considering error.

It is also preferred that the second generation unit includes a unit for dividing the events into groups of events, each of the groups having the number of events equal to or less than a predetermined number.

Preferably, the analysis unit includes a unit for forcibly resetting coefficients used in the VAR model when there is no association rule between the events corresponding to the coefficient.

It is further preferred that the filtering unit includes a unit for verifying reliability of the association rules by using the STP analysis method.

Preferably, the filtering unit includes a unit for forcibly resetting coefficients used in the VAR model when there is no association rule between the events corresponding to the coefficient.

Further objects and advantages of the present invention will be apparent from the following description of the preferred embodiments of the invention as illustrated in the accompanying drawings.

FIG. 1 is a block diagram schematically showing a constitution of an event filtering system according to the present invention.

In FIG. 1, reference numeral 11 denotes an association rule generation unit, 12 a parameter group generation unit, 13 an analysis unit, 14 a filtering unit, and 15 a display unit.

First, the operation of the association rule generation unit 11 is described in detail.

This generation unit 11 detects association rules of events simultaneously occurring within a predetermined period of time (window) to analyze mutual causality. More particularly, the unit 11 recognizes events occurring within a predetermined period of time before and after each sampling time as a transaction and analyzes it.

The unit 11 adopts a modified basket analysis method. The basic basket analysis method extracts only important relationships by defining confidence which depends upon the number of occurred transactions (sets of items), when relationships between items such as purchase of bread and purchase of milk are extracted as an association rule of purchase of bread.fwdarw.purchase of milk.

In the modified basket method, definition of itemset, definition of support in itemset, definition of confidence, production method of important itemset and production method of association rules are modified.

The itemset is defined with considering time sequence. Events occurred within a predetermined period of time t2 from a time point t is considered as a transaction. Thus, the items are events a.sub.t 2,i ={s.sub.t,i : t=[t,t+t2], i=1, . . . , L, with value of 1} occurred within a time period [t,t+t2]. The itemset has to contain the event occurred at a reference time point t. For example, a transaction x may be x=<a.sub.0,1,a.sub.1,5,a.sub.2,6,a.sub.2,15 >where t2=0 means occurrence at the reference time point.

L transactions can be obtained by executing the above prosecutions for the all data at L points with changing the reference time point t. A set of transactions is indicated by D={x.sub.j } where j=1, . . . , L.

The association rule is defined as follows. If partial sets B and H (B.LAMBDA.H={.phi.}) are defined in the itemset x of the transaction x, it can be considered that "H is satisfied when B is satisfied". This relationship is described as "an association rule : B.fwdarw.H". B is called as a body and H is called as a head.

In the basic basket analysis method, although the support is defined by (the number of interested events)/(the number of all events), information with respect to the absolute occurrence number of these interested events is truncated. In the modified basket method, the support of the itemset is defined as follows on the assumption that the events are pursuant to Poisson process. ##EQU1## If the number of events in x is expressed by nx, .sqroot.the number of events in x is expressed by enx, the number of all transactions is expressed by L, and .sqroot.the number of all transactions is expressed by eL, support(x) is expressed as ##EQU2## This is based upon the fact that true average of Poisson distribution with a frequency m is m and that its standard deviation is .sqroot.m.

Definition of confidence in an association rule is substantially the same as that of the basic method except that this system adopt confidence considering its error.

If support is expressed as ##EQU3##

In order to extract useful itemset, the minimum threshold value of confidence (minconf) and the minimum threshold value of certainty (mincert) which are references for adoption/rejection of the itemset are used.

Certainty of itemset Cert(x) is expressed as ##EQU4## wherein .lambda. is the number of occurred itemsets x. The certainty represents cumulative density in the standard deviation (.sqroot..lambda.) in Poisson distribution defined by the average .lambda.. When the number of occurred events increases, this certainty increases and thus ambiguity is so decreased.

The itemsets and the association rules are evaluated as follows.

Itemsets with large uncertainty are rejected. The certainty of events based upon Poisson distribution is defined only by the number of its occurrence. Thus, the uncertainty of each itemset is evaluated in accordance with certainty so as to reject the itemset with a low certainty less than the mincert.

Association rules with a low confidence less than the minconf is also rejected.

Negative association rules are rejected. If confidence of the rule is lower than support of its head (H), this rule is rejected. If first index of change of confidence lift1(B.fwdarw.H)=confidence(B.fwdarw.H)/support(H)(including error) is lower than 1, then the rule is rejected. If confidence of the rule is lowered by adding a condition (D) to body (C), this rule is rejected. If second index of change of confidence lift2(C+D.fwdarw.H)=confidence(C+D.fwdarw.H)/confidence(C.fwdarw.H)(includ ing error) is lower than 1, then the rule is rejected.

Association rule with lower importance item added to its body is also rejected. In order to evaluate independency between the head and the body, Akaike's Information Criteria (AIC) may be used.

Algorithm of the association rule generation unit 11 is as follows.

Step 0: set of important association rules: R={.phi.}

Step 1: set of important items: I.sub.1 ={.phi.}

Step 2: q=1, . . . , N

Step 3-1: y.sub.q ={a.sub.q }, the number of occurrence (=n(y.sub.q)) is counted

Step 3-2: itemset containing no item occurred at the reference time point (a.sub.0,*) is rejected

Step 4: if Cert(y.sub.q)<mincert, then go to Step 2 (itemset with large uncertainty is rejected)

Step 5: y.sub.q is included into I.sub.1, go to Step 2

Step 6: p=2, . . . , N

Step 7: p items set is generated by composition of sets in

I.sub.p-1 : Y.sub.p

Y.sub.p ={y.sub.q :q=1, . . . M.sub.p }, M.sub.p : the number of sets contained

in the set Y.sub.p

Step 8: set of important items: I.sub.p ={.phi.}

Step 9: q=1, . . . , M.sub.p

Step 10-1: y.sub.q ={a.sub.q }, the number of occurrence (=n(y.sub.q)) is counted

Step 10-2: itemset containing no item occurred at the reference time point (a.sub.0,*) is rejected

Step 11: if Cert(y.sub.q)<mincert, then go to Step 9 (itemset with large uncertainty is rejected)

Step 12: flag=0

Step 13: r=1, . . . , (p-1)

Step 14: partial set of r items set Y.sub.q : Z.sub.r Z.sub.r ={z.sub.s :s=1, . . . , 0.sub.r }, 0.sub.r : the number of sets contained in the set Z.sub.r

Step 15: s=1, . . . , 0.sub.r

Step 16-1: association rule: B.fwdarw.H, body: B={y.sub.p -z.sub.s } and head: H=z.sub.s are defined

Step 16-2: association rule B={y.sub.p -z.sub.s } containing no item occurred at the reference time point (a.sub.0,*) is rejected, go to Step 15

Step 17: if confidence(B.fwdarw.H)<minconf, then go to Step 15

(association rule with low confidence is rejected)

Step 18: lift1(B.fwdarw.H)=confidence(B.fwdarw.H)/support(H) if lift1(B.fwdarw.H)<1 then go to Step 15 (negative association rule is rejected)

Step 19: if AIC(DM)>AIC(IM) then go to Step 15 (association rule with no dependency between head and body is rejected)

Step 20: t=1, . . . , (r-1)

Step 21: partial set of t items set B: W.sub.t

W.sub.t ={w.sub.u :u=1, . . . , P.sub.t }, P.sub.t : the number of sets contained in the set W.sub.t

Step 22: u=1, . . . , P.sub.t

Step 23: C={B-w.sub.u } and D=w.sub.u are defined

Step 24: if set of rules R contains no association rule of C.fwdarw.H, then go to Step 22

Step 25: lift2(C+D.fwdarw.H)=confidence(C+D.fwdarw.H)/ confidence(C.fwdarw.H)

if lift2(C+D.fwdarw.H)<1 then go to Step 15 (negative association rule is rejected)

Step 26: if AIC(DM)>AIC(IM) then go to Step 15

(association rule with added lower importance item is rejected)

Step 27: go to Step 22

Step 28: go to Step 20

Step 29 association rule of B.fwdarw.H is included in the set of rules R, flag=1

Step 30 go to Step 15

Step 31: go to Step 13

Step 32: if flag=1, then itemset y.sub.q is included in the important itemset I.sub.p

Step 33: go to Step 9

Step 34: go to Step 6

The operation of the parameter group generation unit 12 shown in FIG. 1 is now described in detail.

This generation unit 12 divides parameters (events) into groups, each including a predetermined or less number of parameters used for VAR (Vector or Multivariate Auto-Regressive model), in accordance with the association rules between the parameters (events), derived by the unit 11. This division is necessary because the VAR method has the upper limit of the number of parameters to be analyzed.

In the unit 12, the graphic separation method according to the graph theory is basically used. Basic matrix of cut set of the objected graph, which is used for a reference at the graphic separation is generated, and then groups of the parameters with the upper limit of the parameter number are generated by the division.

First, definitions with respect to graph are described.

graph: dots (vertexes) and edges connecting the dots

vertex set: V={v.sub.1, . . . , v.sub.i, . . . , v.sub.I } (v.sub.i : each vertex, I: total number of vertexes in the graph)

edge: e.sub.j ={v.sub.p, v.sub.q } (e.sub.j : line across v.sub.p and v.sub.q)

edge set: E={e.sub. 1, . . . , e.sub.j, . . . , e.sub.J } (e.sub.j : each edge, J: total number of edges in the graph)

formal representation of the graph: G(V,E) (this graph consists of vertexes V and edges E)

tree: non-directed graph wherein all vertex pairs are just connected by one column of edges

generation tree: tree including all vertexes in the graph

chord if T represents the tree of the graph, edge which is not included in T is called as chord with respect to T

closed loop: column of edges starting from a point and terminated to the same point without including the same edge twice

cut set: edge set E.sub.1 if this edge set is removed from the connected graph G(V,E), the graph becomes a non-connected graph.

Next, definitions of matrixes accompanying with the graph are described.

connection matrix (B) (matrix B is a connection matrix of graph G(V,E)): each row of B corresponds to vortex of G, each column of B corresponds to edge of G, b.sub.ij =1 when j-th edge is connected with i-th vortex and b.sub.ij =0 when otherwise, and all columns corresponding to loops are 0

closed loop matrix (C) (matrix C is a closed loop matrix of graph G(V,E)): each row of C corresponds to closed loop of G, each column of C corresponds to edge of G, and c.sub.ij =l when j-th edge is included in i-th closed loop and c.sub.ij =0 when otherwise

cut set matrix (K) (matrix K is a cut set matrix of graph G(V,E)) each row of K corresponds to cut set of G, each column of K corresponds to edge of G, and k.sub.ij =1 when j-th edge is included in i-th cut set and k.sub.ij =0 when otherwise.

The maximum number of parameters in each group finally obtained is defined as maximum parameter number (MP). Since the number of variants allowed for the VAR method is up to 5, it is desired to determine as MP=5 in general.

(1) In order to extract groups of the parameters, first, composite rules are transformed to simple rules. The composite rule is an association rule in which its head or its body is constituted by a plurality of parameters, for example [Tires].fwdarw.[Automotive Services]. Whereas the simple rule is an association rule in which each of its head and body is constituted by a single parameter, for example [Tires] and [Auto Accessories].fwdarw.[Automotive Services]. Since the VAR method cannot evaluate composite rule, it is necessary to transform such composite rule to simple rules if the association rule contains composite rule. This transformation is executed by searching composite association rule in which its head or its body is constituted by a plurality of parameters, and by generating all combinations of simple rules which are defined by the respective parameters in the body and by the respective parameters in the head. No generation of simple rule is necessary if this simple rule has been existed already. Support and confidence of the newly generated simple rules depend on those of the composite rule.

(2) Then, directed graphs are changed into non-directed graphs. The association rules are expressed as the directed graphs in which directions are defined for indicating relationships between parameters. Since such directional information is unnecessary to this end, graph is changed to non-directed graph. This change is executed by canceling directional information in the association rule (relationships between parameters are expressed as pairs of two parameters), by adopting one rule when there two bidirectional association rules between two parameters and adopting larger support and confidence, by translating the parameters into vortexes in the graph in accordance with the graph theory, and finally by expressing all relationships on the graphs by means of connection matrix.

(3) Next, connected graphs are extracted from all the graphs. Group of connected parameters is focussed from all parameters. This extraction is executed by checking all the edges and extracting groups of the connected vortexes to express the connected graphs constituted by the respective groups of the connected vortexes as G.sub.n (n=1, . . . , N, N is the total number of the connected graphs), by generating a connection matrix B.sub.n : n=1, . . . , N for each of the connected graphs G.sub.n, and by storing partial graphs contained in each graph G.sub.n and provided with the number of the vortexes equal to or less than MP into "a set of final partial graphs R".

(4) Then, basic matrix of cut sets is generated from the graphs. Namely, for each of the graphs G.sub.n, its cut set K.sub.n is extracted. To this end, a generation tree is generated from the connection matrix of the graph, then a basic matrix of closed loops is generated, and finally the basic matrix of cut sets is generated.

(4-1) The generation tree (T) is obtained from the connection matrix (B) as follows. The number of the vortexes is indicated as I.

Step 1: edge set T to be included in the generation tree is determined. Initially, T=0. Labels of "0" are attached to all the edges.

Step 2 edge e.sub.1 is included in T. Both vortexes connected to e.sub.1 are checked. Row of e.sub.1 in the connection matrix B is canceled.

Step 3: 1 is added to labels of all edges e.sub.j which are elements b.sub.ij =1 of B, for already checked vortexes v.sub.i.

Step 4: all edges with label of 2 are canceled. If there is no edge with label of 1, go to Step 5. Otherwise, an optional edge e.sub.j with label of 1 within the already checked raw is selected and unchecked vortex which is connected with the selected edge is checked. Then, e.sub.j is included in T and the raw including e.sub.j is canceled from the matrix.

Step 5 when (I-1) edges are entered in T, this T provides the generation tree. If all the labels are 0, there exists no generation tree. Otherwise, go to Step 3.

(4-2) The basic matrix of closed loop (C) is obtained from the connection matrix (B) by utilizing the generation tree (T) as follows.

Step 1: p is defined as index representing the number of rows of closed loop basic matrix. Initially, p=1. A set of chords of the generation tree is determined as F=E-T.

Step 2: p-th chord in F is selected to determine it as an edge e.sub.k. Rows in B corresponding to all chords other than e.sub.k are canceled.

Step 3: for i=1, . . . , (the number of vortexes)

Step 4: for j=1, . . . , (the number of edges)

Step 5 column i and row j are canceled when elements b.sub.ij and b.sub.im (m.noteq.j) of B are b.sub.ij =1 and b.sub.im =0. Go to Step 4.

Step 6: Go to Step 3.

Step 7: c.sub.pq =1 for all rows q which were not canceled.

Step 8: p=p+1. x is determined to a value of the number of edges minus the number of vortexes. If p is equal to x+2, the process is ended. C at that time is the closed loop basic matrix to be obtained. If p is smaller than x+2, all the canceled rows and columns of B are recovered. Then, go to Step 2.

(4-3) The basic matrix of cut set (K) is obtained from the connection matrix (B) by utilizing the closed loop basic matrix (C) as follows.

Step 1: column which forms a unit matrix is discovered in the closed loop basic matrix (C). By canceling this column, C.sub.0 is obtained.

Step 2: a matrix of ##EQU5## is generated by adding the unit matrix to C.sub.0 as a row. The cut set matrix is obtained as the transposed matrix of this matrix.

(5) Then, partial graphs are extracted by using the basic matrix of the cut set. The cut set represents how to cancel edges for dividing the graph. Thus, by utilizing the cut set, the partial graph including vortexes of the number equal to or less than a determined value is obtained as follows.

Step 1: basic cut set contained in the cut set basic matrix (K) is canceled if the number of its vortexes which are contact with the edge set is more than MP.

Step 2: For each basic cut set, the maximum values of support and confidence between the contacting vortexes are confirmed.

Step 3: For each basic cut set, priority is determined in decreasing order of the number of vortexes, the maximum support and the maximum confidence.

Step 4: from the edge set E, edges contained in the basic cut set are canceled in decreasing order of the priority. However, if any one of the edge to be canceled has been already canceled, this basic cut set is rejected.

Step 5: the vortexes contacting with the basic cut set which can be canceled is confirmed, and then the partial graph constituted only by these vortexes is accommodated into "a set of final partial graphs R".

Step 6: When E={.phi.} or K={.phi.}, the process is ended. In the latter case, the remainder of the set E is stored into R. However, when the number of the remaining edge is equal to 1, the reminder is rejected without storing into R.

(6) Finally, groups of parameters are outputted. For each graph contained in the set of final partial graphs R, groups of contacting vortexes are provided. These provided vortexes are groups of parameters to be obtained.

The operation of the analysis unit 13 shown in FIG. 1 is now described in detail.

This analysis unit 13 executes the VAR analysis with respect to the focused groups of parameters. In the analysis unit 13, the VAR model is used as the prediction model. By using the VAR model, an appropriate model can be effectively obtained from observed data of the parameters. Also, in the unit 13, the method of least squares is utilized for estimating the VAR model. As a result, an approximate value of the maximum likelihood estimator with respect to the auto regressive model can be obtained.

The characteristic variants of the observed data are defined as .DELTA.t: measuring time width (sec), N: the number of observation points, and T: total observation time (=N.times..DELTA.t sec).

In the VAR model, the multivariate time sequence is defined as y.sub.n =[y.sub.n (1),y.sub.n (2), . . . , y.sub.n (L)].sup.t, where t is a transposition matrix, y.sub.n (i) is the value of the parameter #i at time n, and L is the number of the focused variates. Thus, the VAR model which represents the current value of the multivariate time sequence by means of the past values (y.sub.n-1, . . . , y.sub.n-m) and white noise is defined as follows. ##EQU6## where A.sub.i.sup.m ={a.sub.i (p,q)} is an AR coefficient matrix (L.times.L matrix: p=1, . . . , L, q=1, . . . , L), v.sub.n is a L-th dimension white noise, and m is an AR degree which designates what number of the past values should be used.

In the AR method, generally the model is not applied to all of the observed data, and thus a window width is determined. Namely, the whole is divided into several segments and estimator in each of the segments is averaged so as to lower estimation error. The maximum frequency which can be measured is given as f.sub.max =1/2.DELTA.t (Hz) based upon the Nyquist theorem, and the minimum frequency which can be measured is given as f.sub.min =1/T=1/N.DELTA.t (Hz). Thus, the estimation error is in inverse proportion to a square root of the number of the segments. Namely, the greater number of the segments, th lower estimation error. Also, the minimum frequency which can be measured is decreased resulting the observation frequency range to narrow.

A prediction model is constituted by using the VAR model as follows.

Following processes are independently executed for each group V.sub.j of the parameters. Therefore, although the description relates to the evaluation of only one group, all the groups (V.sub.j : j=1, . . . , M) are similarly evaluated.

(a) In order to constitute estimation model of each parameter, the AR model is estimated in accordance with the observed data. The least square method is adopted for the estimation algorithm. Thus, it is possible to use different AR degrees for the respective variates resulting that more flexible modeling than the generally utilized Yule-Walker method can be provided. The optimum AR degree m* is determined by the Akaike's information criteria (AIC).

In the estimation of the VAR model using the least square method, a modified VAR model is used as follows. ##EQU7## where the matrix B.sub.0.sup.m is a lower triangular matrix in which all elements on the diagonal and in its lower area are "0" as follows. ##EQU8##

Each of the coefficient matrixes B.sub.i.sup.m is expressed as; ##EQU9##

The VAR model can be extracted from A.sub.i.sup.m =(I-B.sub.0.sup.m).sup.-1 B.sub.i.sup.m.

Herein after, algorithm of the VAR model estimation is described in detail.

(1) The maximum degree of the AR degree to be evaluated is determined (m).

(2) The following (N-m).times.(Lm+L) matrix X is defined depending upon the observed data. ##EQU10## (3) The following processes are executed to j=1, . . . , L.

(3-1) The matrix X is transformed by means of the Householder transformation which transforms using the orthogonal transformation an optional matrix into an upper triangular matrix to obtain the following form.

If j=1 ##EQU11##

(3-2) For p=1, . . . , m, a residual variance and AIC of the p degree model are calculated as follows:

residual variance of the p degree model: ##EQU12## AIC of the p degree model:

(3-3) The degree p in which the minimum AIC can be obtained is determined to an optimum degree for j-th variate (m*(j)=p).

(3-4) The AR coefficients can be obtained as a solution of the following linear equation. According to the equation, the solution can be easily obtained by backwardly substituting. If j=1 ##EQU13## where Lm* degree vector:

else ##EQU14## where (Lm*+j-1) degree vector:

(3-5) The AR coefficients expressed by C are stored in the corresponding matrixes B.sub.0,B.sub.1, . . . , B.sub.m*, respectively.

(4) The AR coefficients in the general VAR model are calculated.

where i=1, . . . , M, M is the maximum value of the AR coefficients of all the variates.

(a') It is desired that the calculated AR coefficients are reset depending upon the association rules between the parameters provided from the generation unit 11. In fact, the AR coefficient between the parameters with no association rule is forcefully reset to "0". As a result, possible error or apparent association which may be often occurred in the multivariate analysis can be prevented from generation.

The processes of this operation is as follows.

(1) With respect to the focused parameters, combinations each consisting of a pair of the parameters are determined.

(2) For the pair (y(i),y(j)), it is judged whether "association rule of that y(i) is contained in the body and y(j) is contained in the head" exists. If the rule exists, go to (4).

(3) If the association rule does not exist, the AR coefficient which contributes relationship from y(i) to y(j) is reset to "0". In other words, "element on column corresponding to y(i) and on raw corresponding to y(j)" is reset to "0".

(4) For the same pair (y(i),y(j)), it is judged whether "association rule of that y(j) is contained in the body and y(i) is contained in the head" exists. If the rule exists, go to (6).

(5) If the association rule does not exist, the AR coefficient which contributes relationship from y(j) to y(i) is reset to "0". In other words, "element on column corresponding to y(j) and on raw corresponding to y(i)" is reset to "0".

(6) Jump to (2) to consider another pairs. (b) Then, one point forward prediction of each variate is executed by using the VAR model. Namely, the value of the focused variate at current time n is predicted in accordance with the past observed values at times n-1,n-2, . . . , n-m.

The processes of this prediction operation is as follows.

(1) One point forward prediction model is defined by canceling noise source in th VAR model.

The one point forward prediction model is expressed as; ##EQU15## where A.sub.i is an AR coefficient matrix (L.times.L matrix: i=1, . . . , M), and M is the maximum value of AR degrees for all variates.

(2) The one point forward prediction of the respective variates is executed by applying the model to the observed data.

The operation of the filtering unit 14 shown in FIG. 1 is now described in detail. This unit 14 extracts the association rules for filtering. More concretely, this unit verifies reliability of the association rules between the events which are obtained in the unit 11 by means of the Signal Transmission Path (STP) analysis method. By adopting the STP analysis method, if the relationship between the parameters has linearity, the appeared relationship between the parameters can be quantitatively analyzed.

The flow of processes in the filtering unit 14 is (a) estimation of a VAR model, (b) calculation of noise contribution ratio (NPC) based upon the VAR model, (c) calculation of partial noise contribution ratio (PNPC) based upon the VAR model, (d) calculation of partial coherence (PCH) based upon the VAR model, (e) determination of focussing frequency, and (f) evaluation of the association rules by comparing NPC with PNPC-PCH.

The definition of the characteristic variants of the observed data, the multivariate time sequence and the VAR model are the same as those in the VAR analysis.

The sample cross covariance matrix is defined as follows where [y.sub.1 (i)y.sub.2 (i), . . . , y.sub.N (i)],i=1, . . . , L is observed. ##EQU16##

In filtering process of the association rule by means of the STP analysis method, a window width and threshold values for evaluation based upon NPC, PNPC and PCH are determined. The window width is determined as the same as that in the aforementioned VAR analysis.

The threshold values for evaluation based upon NPC, PNPC and PCH are determined by selecting one of the following two definitions.

absolute definition of the threshold value: although relationships between the parameters are evaluated by the STP analysis, the threshold value used is not clear. Thus, in the comparison of the estimators, it is necessary to have an acceptable range of the threshold value. In general, the threshold will be determined to about 10 % of the estimator, but it is desired to determine it by cut-and-try methods.

relative definition of the threshold value: with respect to relationships between the parameters, a predetermined number of threshold values with greater estimators are adopted.

Filtering of the association rules is executed as follows.

Following processes are independently executed for each group V.sub.j of the parameters. Therefore, although the description relates to the evaluation of only one group, all the groups (V.sub.j : j=1, . . . , M) are similarly evaluated.

(a) In order to execute the STP analysis, the AR model is estimated in accordance with the observed data. AR coefficients to be obtained are used as information sources for the STP analysis. The Yule-Walker method is adopted to the estimation algorithm. The optimum AR degree m* is determined by the Akaike's information criteria (AIC).

Herein after, algorithm of the VAR model estimation is described in detail.

(1) The maximum degree of the AR degree to be evaluated is determined (m).

(2) The sample cross covariance matrix for k=0,1, . . . , m is calculated (C.sub.0,C.sub.1, . . . , C.sub.m).

(3) V.sub.0 =U.sub.0 =C.sub.0

AIC.sub.0 =N(L log2.pi.+log.vertline.V.sub.0 .vertline.+L)+L(L+1)

(4) For p=1, . . . , m, ##EQU17## (4.5) AIC.sub.p =N(L log2.pi.+log.vertline.V.sub.p .vertline.+L)+L(L+1)+2L.sup.2 p

where A.sub.i.sup.0 is a L dimensional unit matrix, and D.sup.-1 is an inverse matrix of the matrix D.

(5) m*=p: AIC.sub.p is minimum. At this time, the AR coefficients to be obtained are expressed as A.sub.i.sup.m * (i=1, . . . , m*).

(a') It is desired that the calculated AR coefficients are reset depending upon the association rules between the parameters provided from the generation unit 11. In fact, the AR coefficient between the parameters with no association rule is forcefully reset to "0". As a result, possible error or apparent association which may be often occurred in the multivariate analysis can be prevented from generation.

The processes of this operation is as follows.

(1) With respect to the focused parameters, combinations each consisting of a pair of the parameters are determined.

(2) For the pair (y(i),y(j)), it is judged whether "association rule of that y(i) is contained in the body and y(j) is contained in the head" exists. If the rule exists, go to (4).

(3) If the association rule does not exist, the AR coefficient which contributes relationship from y(i) to y(j) is reset to "0". In other words, "element on column corresponding to y(i) and on raw corresponding to y(j)" is reset to "0".

(4) For the same pair (y(i),y(j)), it is judged whether "association rule of that y(j) is contained in the body and y(i) is contained in the head" exists. If the rule exists, go to (6).

(5) If the association rule does not exist, the AR coefficient which contributes relationship from y(j) to y(i) is reset to "0". In other words, "element on column corresponding to y(j) and on raw corresponding to y(i)" is reset to "0".

(6) Jump to (2) to consider another pairs.

(b) The noise contribution ratio (NPC) is an estimator indicating contribution of a certain noise source causing fluctuation in the focused variate. The NPC with a value closer to 1 indicates that it is the main source of the fluctuation. The NPC is calculated as follows.

(1) A Power Spectral Density (PSD) function is calculated by P(f)=2.DELTA.t{A(f)}.sup.-1 W{A(f)*}.sup.-1, where D* is a complex conjugate of the complex matrix D, A(f) is the AR coefficient expressed in a function of frequency ##EQU18## (I: L.times.L unit matrix, .DELTA.t: observation time width), and W is a covariance matrix of noise. ##EQU19## P(f) is a L.times.L matrix. An (i,i) element P.sub.i,i (f) is called as Auto Power Spectral Density (APSD) of the variate y(i), and (i,j) element p.sub.ij (f) is called as Cross Power Spectral Density (CPSD) of the variates y(i) and y(j).

(2) The NPC is calculated as a noise contribution rate from a noise source of the variate y(j) to the variate y(i) as follows. ##EQU20## where .vertline.D.vertline. is a determinant of the matrix D. (c) The NPC evaluates the influences of the noise sources transmitted through all paths. Thus, the calculated NPC may be sometimes high when there exists no direct path but there exists an indirect path (apparent direct path) over the variate with that noise source. In order to avoid the problems, the partial noise contribution ratio (PNPC) which considers only the direct path between two variates is used. The PNPC is calculated as a partial noise contribution rate from a noise source of the variate y(J) to the variate y(i) as follows. ##EQU21## (d) As similar to the conventional NPC, the general coherence (PCH) which represents degree of association of frequency components in the fluctuation between two variates is also influenced from the indirect paths. Thus, the partial coherence (PCH) which considers only the direct path between two variates is used. The PCH is calculated as a partial coherence from a noise source of the variate y(j) to the variate y(i) as follows. ##EQU22## (e) Since the causality estimators are defined over the entire measuring frequency range, the evaluation is very difficult. Thus, estimation is executed in only focussed important frequency range between the two variates as follows.

(1) The frequency range is selected from one of all frequency range and a particular frequency band.

(2) The important frequency is estimated and selected from one of:

a frequency at which the maximum APSD of the variate for evaluating the noise contribution from another variate is obtained, or several relative maximum values, and

a frequency at which the maximum coherence between the two variates for evaluating the noise contribution is obtained, or several relative maximum values.

The coherence between the variate y(i) and the variate y(j) with respect to (i,j) element p.sub.ij (f) of the PSD function P(f) is defined as:

(f) Screening of the association rules is executed based upon the three causality estimators NPC, PNPC and PCH. First, in accordance with the NPC, it is confirmed whether there are direct or indirect paths. Then, it is confirmed whether there is a direct path between the two parameters in accordance with PNPC-PCH comparison. The evaluation of the association rules are executed for all combinations of two variates y(i) and y(j) of the focused variates or parameters as follows.

(1) If the NPC value at the frequency selected at (e) is within a tolerance error range and NPC.sub.ij (f)=0, it is judged that there is no association rule between these two variates.

(2) If the PNPC and PCH values at the frequency selected at (e) are within respective tolerance error ranges and:

(2-1) PCH.sub.ij =PNPC.sub.ij =PNPC.sub.ij =0, it is judged that there is no association rule between these two variates,

(2-2) PCH.sub.ij =PNPC.sub.ij and PNPC.sub.ji =0, it is judged that there is an association rule of the variate y(i).fwdarw.variate y(j), and

(2-3) PCH.sub.ij =PNPC.sub.ji and PNPC.sub.ij =0, it is judged that there is an association rule of the variate y(j).fwdarw.variate y(j).

(3) If PNPC values at all frequency satisfy, within a tolerance error range, the next two conditions, there is a bidirectional association rule (feedback) between the two variates:

0<PNPC.sub.ij (f)<1 and 0<PNPC.sub.ji (f)<1, and

patterns of PNPC differ from patterns of PCH in all frequency ranges (these can evaluated from the distances between both patterns and inner product values).

(4) The association rule with no contradiction with respect to (1)-(3) is adopted. When combinations of all parameters are evaluated, the process is completed with providing the finally adopted association rule.

The display unit 15 shown in FIG. 1 displays thus extracted association rule between the parameters in graphs which will provide clear understanding. This unit 15 utilizes basically the STT method which is one of automatic plotting and graphing methods. The algorithm of the STT method consists of layering, normalization, determination of arranging order of vortexes, determination of arranging coordinates of vortexes and shaping of graph. Since this display unit 15 can be constituted any another configuration, detailed explanation of the operation in this unit 15 is omitted in the specification.

Many widely different embodiments of the present invention may be constructed without departing from the spirit and scope of the present invention. It should be understood that the present invention is not limited to the specific embodiments described in the specification, except as defined in the appended claims.