Abstract:
A method providing an analytical technique introducing label information into an anomaly detection model. Effective utilization of label information is based on introducing the degree of similarity between samples. Assuming, for example, there is a degree of similarity between normally labeled samples and no similarity between normally labeled and abnormally labeled samples. Also each sensor value is generated by the linear sum of a latent variable and a coefficient vector specific to each sensor. However, the magnitude of observation noise is formulated to vary according to the label information for the sensor values, and set so that normal label≦unlabeled≦anomalously labeled. A graph Laplacian is created based on the degree of similarity between samples, and determines the optimal linear transformation matrix according to a gradient method. A optimal linear transformation matrix is used to calculate an anomaly score for each sensor in the test samples.

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application claims priority under 35 U.S.C. §119 from Japanese Patent Application No. 2012-134319 filed Jun. 14, 2012, the entire contents of which are incorporated herein by reference. 
     BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present invention relates to a technique for detecting anomalies in vehicles, industrial machinery, and the like. 
     2. Description of the Related Art 
     Because an accident occurring in industrial machinery at a railway or plant has significant social consequences, it is very important to detect any anomaly that can occur before an accident occurs. 
     In order to ensure safety, sensors have been installed in vehicles and industrial machinery at various locations to monitor operations, and the measurement data obtained from these sensors has been analyzed by computers to detect anomalies. 
     For example, the temperature at major locations in a vehicle can be used to detect anomalies. The temperature can be measured using a laser measuring device installed near the path of a vehicle. In this way, early anomaly detection can be performed based on the measured data. 
     Here, knowledge related to devices used to detect anomalies is incorporated into the computers performing the analysis. 
     However, anomaly detection in our knowledge base has not yet reached the point of being sufficiently reliable. At this point, the reliability of anomaly detection is being increased by using anomaly patterns detected in the past and reducing the possibility of overlooking cases similar to those in the past. The related technologies described in the following literature are known. 
     Laid-open Patent Publication No. 7-280603 describes the use of samples in an anomaly detecting method for machinery. 
     International Patent Publication No. WO2008/114863 describes the calculation of the degree of similarity between patterns of change in objects observed using diagnostic equipment. 
     Laid-open Patent Publication No. 2008-58191 describes the calculation of the degree of similarity between standard parameter values as confidence factors in a diagnostic method for rotating machinery. 
     Laid-open Patent Publication No. 2009-76056 describes the use of anomaly frequency measurements in a method used to identify anomalous values. 
     Laid-open Patent Publication No. 2010-78467 describes a method in which a correlation coefficient matrix is created with time-series data for testing purposes and normal time-series data for reference purposes, a sparse accuracy matrix is created in which each correlation coefficient matrix is an inverse matrix, and a localized probability distribution is created for the time-series data for testing purposes and the normal time-series data for reference purposes, preferably using the accuracy matrix in a multivariate Gaussian model. 
     X. Zhu, Z. Ghahramani, “Semi-Supervised Learning Using Gaussian Fields and Harmonic Functions” in Proceedings of the ICML, 2003 describes semi-supervised learning based on a Gaussian random field model, and discloses labeled data and unlabeled data represented as vertices in a weighted graph. 
     A. B. Goldberg, X. Zhu, and S. Wright, “Dissimilarity in Graph-Based Semi-Supervised Classification” in AISTATS, 2007 describes a semi-supervised classification algorithm in which learning occurs based on the degree of similarity and dissimilarity between labeled data and unlabeled data. 
     SUMMARY OF THE INVENTION 
     Various techniques applicable to anomaly detection have been described, including those with semi-supervised algorithms, but none have suggested the use of anomaly patterns detected in the past. In other words, the effective utilization of anomaly patterns detected in the past requires arbitrary preprocessing in prior art techniques, but this does not sufficiently increase the reliability of anomaly detection. 
     Therefore, it is an object of the present invention to provide an analytical technique introducing existing label information into an anomaly detection model. 
     It is another object of the present invention to provide an anomaly detection technique able to effectively utilize label information in data including a mix of both labeled samples and unlabeled samples. 
     The present invention is intended to solve these problems. The effective utilization of label information is based on the idea of introducing the degree of similarity between samples. This assumes, for example, that there is a degree of similarity between normally labeled samples and no similarity to abnormally labeled samples. 
     Also, it is assumed that an unlabeled sample has greater a degree of similarity to a normal sample than to an anomalous sample when it has been determined from past experience that a failure is unlikely to occur, and that an unlabeled sample has an equal degree of similarity to a normal sample and to an anomalous sample when there is no previous information. 
     Each normalized sample is expressed by a multi-dimensional vector in which each element is a sensor value. The present invention also assumes that each sensor value is generated by the linear sum of a latent variable and a coefficient vector specific to each sensor. However, the magnitude of observation noise is formulated to vary according to the label information for the sensor values. The observation noise is set so that normally labeled≦unlabeled≦anomalously labeled. 
     Next, a graph Laplacian is created based on the degree of similarity between samples, the graph Laplacian is used to determine the optimal linear transformation matrix according to the gradient method or the like. 
     When the optimal linear transformation matrix has been obtained, an anomaly score is calculated for each sensor in the test samples according to the technique described in the Patent Application No. 2011-206087 filed by the present applicant. 
     The present invention is able to reduce the arbitrariness of criteria for anomaly detection and increase the reliability of anomaly detection by incorporating samples of anomaly patterns and normal patterns detected in the past into an anomaly detection model. 
    
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         FIG. 1  is a block diagram showing an example of a computer hardware configuration used to embody the present invention. 
         FIG. 2  is a block diagram showing a function configuration used to embody the present invention. 
         FIG. 3  is a flowchart of the process for calculating model parameters for anomaly detection according to the present invention. 
         FIG. 4  is a flowchart of the process for calculating anomaly scores using model parameters and the like. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
     The following is an explanation of an example of the present invention with reference to the drawings. The same reference numbers are used to denote the same objects in all of the drawings except where otherwise indicated. The following explanation is a single embodiment of the present invention. The present invention is by no means intended to be limited to the content explained in the example. 
       FIG. 1  is a block diagram of computer hardware used to realize the system configuration and processing in an example of the present invention. In  FIG. 1 , a CPU  104 , main memory (RAM)  106 , a hard disk drive (HDD)  108 , a keyboard  110 , a mouse  112 , and a display  114  are connected to a system bus  102 . The CPU  104  is preferably based on 32-bit or 64-bit architecture, and can be a Pentium® 4, Core® 2 Duo or Xeon® from Intel Corporation, or an Athlon® from Advanced Micro Devices, Inc. The main memory  106  preferably has a capacity of 4 GB or more. The hard disk drive  108  preferably has a capacity of 500 GB or more in order to store a large amount of data. 
     While not shown in any of the drawings, the hard disk drive  108  includes a pre-installed operating system. The operating system can be any operating system compatible with the CPU  104 . Examples include Windows XP® or Windows® 7 from Microsoft Corporation, or MacOS® from Apple, Inc. 
     The hard disk drive  108  also contains, as explained below with reference to  FIG. 2 , a main program  202 , labeled data  204 , unlabeled data  206 , a parameter group  208 , a graph Laplacian calculation routine  210 , a parameter optimization routine  212 , and an anomaly detection routine  214 . The main program  202 , graph Laplacian calculation routine  210 , parameter optimization routine  212 , and anomaly detection routine  214  can be written in any existing programming language, including Java®, C, C ++ , or C#. 
     The keyboard  110  and mouse  112  operate on the operating system or main program  202  loaded from hard disk drive  108  into main memory  106  and displayed on display  114 , and are used to enter characters. 
     The display  114  is preferably a liquid crystal display. Any resolution can be used, including XGA (resolution: 1024×768) or UXGA (resolution: 1600×1200). While not shown in the drawings, display  114  is used to display operating windows for entering parameters and starting programs, and to display parameter calculation results and the like. 
     The following is an example of a functional configuration of the processing in the present invention with reference to the block diagram in  FIG. 2 . 
     In  FIG. 2 , the main program  202  is a program with functions integrating all of the processing. This is used by the operator to set a parameter group  208 , start the graph Laplacian calculation routine  210 , parameter optimization routine  212 , and anomaly detection routine  214 , execute calculations, and display results on display  114 . 
     Labeled data  204  includes data detected in the past that has been found to be anomalous or normal. An anomaly label is applied to data found to be anomalous, and a normal label is applied to data found to be normal. 
     Unlabeled data  206  includes unlabeled data that has not been found to be either anomalous or normal. Depending on the situation, it is treated as either labeled data  204  or unlabeled data  206 . 
     A single unit of data (called a sample) is a D-dimensional real vector consisting of type-D sensor values. A set of sensor data can be expressed by the equation X=[X1, . . . , XN] T εR N×D , where N is the number of samples. 
     Sensor data set X is preferably data normalized based on the original sensor data set X′=[X′1, . . . , X′N] T εR N×D . The normalization is performed based on the following equation. Here, Xn,d is the d th  element of vector Xn. The same is true of X′n,d. 
     
       
         
           
             
               
                 
                   
                     
                       μ 
                       d 
                     
                     = 
                     
                       
                         1 
                         N 
                       
                       ⁢ 
                       
                         
                           ∑ 
                           
                             n 
                             = 
                             1 
                           
                           N 
                         
                         ⁢ 
                         
                           X 
                           
                             n 
                             , 
                             d 
                           
                           ′ 
                         
                       
                     
                   
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     
                       σ 
                       d 
                     
                     = 
                     
                       
                         
                           1 
                           N 
                         
                         ⁢ 
                         
                           
                             ∑ 
                             
                               n 
                               = 
                               1 
                             
                             N 
                           
                           ⁢ 
                           
                             
                               ( 
                               
                                 
                                   X 
                                   
                                     n 
                                     , 
                                     d 
                                   
                                   ′ 
                                 
                                 - 
                                 
                                   μ 
                                   d 
                                 
                               
                               ) 
                             
                             2 
                           
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       X 
                       
                         n 
                         , 
                         d 
                       
                     
                     = 
                     
                       
                         
                           X 
                           
                             n 
                             , 
                             d 
                           
                           ′ 
                         
                         - 
                         
                           μ 
                           d 
                         
                       
                       
                         σ 
                         d 
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   1 
                 
               
             
           
         
       
     
     Also, label information Y=[Y1, . . . , YN] T εR N×D  is provided for each sensordata set X=[X1, . . . , XN] T εR N×D . While not shown in the drawings, this is stored along with the labeled data  204  and the unlabeled data  206  in the hard disk drive  108 . 
     The label information Y is defined as follows. 
     
       
         
           
             
               
                 
                   
                     Y 
                     
                       n 
                       , 
                       d 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           0 
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               normal 
                             
                             , 
                           
                         
                       
                       
                         
                           1 
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               anomaly 
                             
                             , 
                           
                         
                       
                       
                         
                           NaN 
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               unlabel 
                             
                             , 
                           
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   2 
                 
               
             
           
         
       
     
     Here, NaN is any real number other than 0 or 1. 
     In the present invention, it is assumed that each sensor value Xn,d in each normalized sample Xn is expressed as follows with a latent variable ZnεR D′  (D′≦D), coefficients for the magnitude of noise for each label, snormal, sanomaly, sunlabel, and Gaussian noise c with a mean of 0 and a variance of 1. Here, snormal corresponds to normal, sanomaly corresponds to anomalous, and sunlabel corresponds to unlabeled. Also, D′ is usually equal to D, but D′ is set to about 100 when D is very large and the number of data units N is small. 
     
       
         
           
             
               
                 
                   
                     X 
                     
                       n 
                       , 
                       d 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               
                                 
                                   W 
                                   d 
                                   T 
                                 
                                 ⁢ 
                                 
                                   Z 
                                   n 
                                 
                               
                               + 
                               
                                 
                                   s 
                                   normal 
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 ɛ 
                               
                             
                             , 
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 Y 
                                 
                                   n 
                                   , 
                                   d 
                                 
                               
                             
                             = 
                             0 
                           
                         
                       
                       
                         
                           
                             
                               
                                 
                                   W 
                                   d 
                                   T 
                                 
                                 ⁢ 
                                 
                                   Z 
                                   n 
                                 
                               
                               + 
                               
                                 
                                   s 
                                   anomaly 
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 ɛ 
                               
                             
                             , 
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 Y 
                                 
                                   n 
                                   , 
                                   d 
                                 
                               
                             
                             = 
                             1 
                           
                         
                       
                       
                         
                           
                             
                               
                                 
                                   W 
                                   d 
                                   T 
                                 
                                 ⁢ 
                                 
                                   Z 
                                   n 
                                 
                               
                               + 
                               
                                 
                                   s 
                                   unlabel 
                                 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 ɛ 
                               
                             
                             , 
                           
                         
                         
                           otherwise 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   3 
                 
               
             
           
         
       
     
     Here, the setting is snormal sunlabel sanomaly. Specific examples include snormal=1, sunlabel=3, sanomaly=5 if nothing is found; snormal=1, sunlabel=2, sanomaly=5 if the unlabeled data is found to be mostly normal; and snormal=1, sunlabel=4, sanomaly=5 if the unlabeled data is found to be mostly anomalous. 
     The parameter group  208  includes parameters such as noise magnitudes snormal, sanomaly, sunlabel, a scale parameter A, and the numbers of dimensions D, D′. These are stored in the hard disk drive  108 , and can be set by the user. 
     The parameter group  208  also includes values used to determine a similarity matrix R. The similarity matrix is a N×N square matrix, where N is the number of samples, each row and each column correspond to samples (for example, row i/column j corresponds to the degree of similarity between the i th  and j th  samples), an element corresponding to a normal (labeled) sample and a normal sample is positive number a, an element corresponding to a normal sample and an anomalous sample is non-positive number b, an element corresponding to a normal sample and an unlabeled sample is c, an element corresponding to an anomalous sample and an anomalous sample is d, an element corresponding to an anomalous sample and an unlabeled sample is e, and an element corresponding to an unlabeled sample and an unlabeled sample is f. Here, a, b, c, d, e, and f satisfy the relationships b≦c≦a and e≦d≦f. 
     Preferably, a, b and d above are set as a=1, b=0, d=0.2. 
     As in the case of sunlabel, the settings for c, e and f depend on what the algorithm user has discovered regarding the unlabeled data in the application data. 
     Namely: 
     c=0.5, e=0.1, f=0.5, for example, if nothing is found; 
     c=0.8, e=0, f=0.8, for example, if the unlabeled data is found to be mostly normal; and 
     c=0, e=0.1, f=0.2, for example, if the unlabeled data is found to be mostly anomalous. 
     The graph Laplacian calculation routine  210  creates a similarity matrix R based on the values a, b, c, d, e, f set in the parameter group  208 , and then calculates a graph Laplacian L from the resulting similarity matrix R in the following way. 
     
       
         
           
             
               
                 
                   
                     
                       K 
                       
                         i 
                         , 
                         i 
                       
                     
                     = 
                     
                       
                         ∑ 
                         
                           d 
                           = 
                           1 
                         
                         N 
                       
                       ⁢ 
                       
                         R 
                         
                           i 
                           , 
                           d 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     L 
                     = 
                     
                       K 
                       - 
                       R 
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   4 
                 
               
             
           
         
       
     
     The latent variable Z≡[Z1, . . . , ZN] T εR N×D  is realized by means of the graph Laplacian L as follows. 
     
       
         
           
             
               
                 
                   
                     Pr 
                     ⁡ 
                     
                       ( 
                       Z 
                       ) 
                     
                   
                   ∝ 
                   
                     exp 
                     ⁢ 
                     
                       { 
                       
                         
                           - 
                           
                             λ 
                             2 
                           
                         
                         ⁢ 
                         
                           tr 
                           ⁡ 
                           
                             ( 
                             
                               
                                 Z 
                                 T 
                               
                               ⁢ 
                               LZ 
                             
                             ) 
                           
                         
                       
                       } 
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   5 
                 
               
             
           
         
       
     
     Because the probability Pr(X|W, Z, s) of X≡[X1, . . . XN] T εR N×D  can be regarded as a likelihood function of parameter W≡[W1, . . . WD] T εR D×D′  and Z, parameter optimization routine  212  seeks (W*, Z*) using, for example, the gradient method so that the posterior probability is optimized. This process will be explained in greater detail below with reference to the flowchart in  FIG. 3 . 
     The anomaly detection routine  214  calculates the anomaly score for each variable based on (W*, Z*) obtained in this manner. The anomaly detection routine  214  preferably uses the technique described in Patent Application No. 2011-206087 filed by the present applicant. The processing in the anomaly detection routine  214  will be explained in greater detail below with reference to the flowchart in  FIG. 4 . 
     The following is an explanation of the processing used to determine the model parameters (optimal linear transformation matrix) W* and the like with reference to the flowchart in  FIG. 3 . 
     In Step  302  of  FIG. 3 , main program  202  inputs training data {X′εRN×D,y} by retrieving labeled data  204  and unlabeled data  206  from hard disk drive  108 , normalizes the data in the manner described above, and stores the mean μd and standard deviation ad of each column d calculated when each column was normalized. 
     In Step  304 , main program  202  retrieves parameter D′, scale parameter λ, snormal, sunlabel, sanomaly, a, b, c, d, e, and f from the parameter group  208  or enters them into a setting screen (not shown) using keyboard  110  and mouse  112 . The scale parameter λ can be set, for example, to 0.1, and the noise magnitude and the like are determined as indicators using the cross-validation method. 
     In Step  306 , main program  202  calls up the graph Laplacian calculation routine  210 , and a graph Laplacian L is calculated using label information Y and a, b, c, d, e, and f. Because the graph Laplacian L calculation has already been explained with reference to  FIG. 2 , further explanation has been omitted here. 
     In Step  308 , the main program  202  initiates WεR D×D′  and ZεR N×D′ . Any method can be used to perform the initialization. However, W and Z are initialized with a standard normal distribution, that is, a value of each element of W or Z is set to a realized value of a normal distribution in which the mean is 0 and the standard deviation is 1. 
     In Step  310 , the main program  202  sets the time variable t to 1. 
     In Step  312 , the main program  202  updates W in accordance with the following equation.
 
 W:=W−α[{S ·( X−ZW   T )} T   Z+N ( WW   T ) −1   W]   Equation 6
 
     Here, S is described as follows. 
     
       
         
           
             
               
                 
                   
                     
                       { 
                       S 
                       } 
                     
                     
                       n 
                       , 
                       d 
                     
                   
                   = 
                   
                     { 
                     
                       
                         
                           
                             
                               1 
                               
                                 s 
                                 normal 
                               
                             
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 Y 
                                 
                                   n 
                                   , 
                                   d 
                                 
                               
                             
                             = 
                             0 
                           
                         
                       
                       
                         
                           
                             
                               1 
                               
                                 s 
                                 anomaly 
                               
                             
                           
                         
                         
                           
                             
                               if 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 Y 
                                 
                                   n 
                                   , 
                                   d 
                                 
                               
                             
                             = 
                             1 
                           
                         
                       
                       
                         
                           
                             
                               1 
                               
                                 s 
                                 unlabel 
                               
                             
                           
                         
                         
                           otherwise 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   7 
                 
               
             
           
         
       
     
     The operation S·(X−ZW T ) means elements n, d of matrix (X−ZW T ) are multiplied by elements n, d of S. Also, a is the learning rate and is set, for example, to 0.1. The value of a needs not be constant. It can be reduced with each iteration. 
     In Step  314 , the main program  202  updates Z in accordance with the following equation.
 
 Z:=Z−α[{S ·( X−ZW   T )} W+λLZ]   Equation 8
 
     This equation is used to perform calculations so that the parameters are updated in accordance with an update equation with a term that reduces the penalty based on the degree of similarity. This includes a term that reduces the penalty based on the degree of similarity to the latent variable of each observation. More specifically, it has been formulated so that the penalty based on the degree of similarity is the Mahalanobis distance based on the similarity matrix (or graph Laplacian). It is then calculated to converge in accordance with the gradient method. 
     Step  312  and Step S 314  do not have to be calculated in this order. The order can be switched. 
     After Step  314 , main program  202 , in Step  316 , determines the termination conditions. Here, the Frobenius norm is calculated for the matrix W′ calculated in the previous loop and the matrix W calculated in the current loop, and the termination conditions are satisfied when this is within, for example, 0.001 of a predetermined threshold value. 
     
       
         
           
             
               
                 
                   
                     
                        
                       
                         
                           W 
                           ′ 
                         
                         - 
                         W 
                       
                        
                     
                     F 
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           1 
                         
                         D 
                       
                       ⁢ 
                       
                         
                           ∑ 
                           
                             j 
                             = 
                             1 
                           
                           
                             D 
                             ′ 
                           
                         
                         ⁢ 
                         
                           
                             ( 
                             
                               
                                 W 
                                 
                                   i 
                                   , 
                                   j 
                                 
                                 ′ 
                               
                               - 
                               
                                 W 
                                 
                                   i 
                                   , 
                                   j 
                                 
                               
                             
                             ) 
                           
                           2 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   9 
                 
               
             
           
         
       
     
     In Step  318 , main routine  202  increases t by “1”, and returns to Step  312  when the termination conditions have not been satisfied. 
     In Step  320 , main program  202  outputs the model parameters W, snormal, μ=[μ1, . . . , μD], and σ=[σ1, . . . , σD] when the termination conditions have been satisfied. 
     The following is an explanation of the anomaly score calculation processing performed in anomaly detection routine  214  with reference to the flowchart in  FIG. 4 . 
     In Step  402  of  FIG. 4 , main program  202  calls up anomaly detection routine  214 , and provides model parameters W, snormal, μ=[μ1, . . . , μD], σ=[σ1, . . . , σD]. 
     In Step  404 , the anomaly detection routine  214  inputs test data {X′εR N×D , y} by retrieving labeled data  204  and unlabeled data  206  from hard disk drive  108 , X′ in each column is normalized according the equation described above using μ and σ, and X is obtained. 
     In Step  406 , the anomaly detection routine  214  calculates the correlation anomaly score vector snεR D  using the following equation.
 
 s   n   ≡s   0 +½diag( AX   n   X   n   T   AB   −1 )  Equation 10
 
     Provided, 
     
       
         
           
             
               
                 
                   
                     A 
                     = 
                     
                       
                         ( 
                         
                           
                             
                               W 
                               T 
                             
                             ⁢ 
                             W 
                           
                           + 
                           
                             
                               s 
                               normal 
                               2 
                             
                             ⁢ 
                             I 
                           
                         
                         ) 
                       
                       
                         - 
                         1 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     B 
                     ≡ 
                     
                       
                         diag 
                         2 
                       
                       ⁡ 
                       
                         ( 
                         A 
                         ) 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       
                         ( 
                         
                           s 
                           0 
                         
                         ) 
                       
                       i 
                     
                     ≡ 
                     
                       
                         1 
                         2 
                       
                       ⁢ 
                       ln 
                       ⁢ 
                       
                         
                           2 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           π 
                         
                         
                           A 
                           
                             i 
                             , 
                             i 
                           
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   11 
                 
               
             
           
         
       
     
     Here, I is a unit matrix. The algorithm used to calculate the correlation anomaly score vector based on the optimal linear transformation matrix W is described in Patent Application No. 2011-206087 filed by the present applicant. It is not described in detail here. 
     In Step  408 , anomaly detection routine  214  outputs anomaly score vectors s1, . . . , sN based on these calculations. Each element of s1, . . . , or sN is an anomaly score for each sensor of the first, second, . . . , or N th  test sample, that is, each dimension of s=each variable. A higher value indicates an anomaly. 
     The anomaly detection for industrial machinery at a railway or plant in the present invention was explained with reference to an example. However, the present invention is not limited to this. It can be applied to any example in which anomaly detection is performed based on a plurality of measurement parameters.