WO2023223461A1 - 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム - Google Patents

因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム Download PDF

Info

Publication number
WO2023223461A1
WO2023223461A1 PCT/JP2022/020680 JP2022020680W WO2023223461A1 WO 2023223461 A1 WO2023223461 A1 WO 2023223461A1 JP 2022020680 W JP2022020680 W JP 2022020680W WO 2023223461 A1 WO2023223461 A1 WO 2023223461A1
Authority
WO
WIPO (PCT)
Prior art keywords
series data
causal relationship
time series
time
variables
Prior art date
Application number
PCT/JP2022/020680
Other languages
English (en)
French (fr)
Inventor
良行 乗松
Original Assignee
三菱電機株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 三菱電機株式会社 filed Critical 三菱電機株式会社
Priority to PCT/JP2022/020680 priority Critical patent/WO2023223461A1/ja
Priority to JP2024513679A priority patent/JP7483180B2/ja
Publication of WO2023223461A1 publication Critical patent/WO2023223461A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N99/00—Subject matter not provided for in other groups of this subclass

Definitions

  • the present disclosure relates to a causal relationship estimation device, a causal relationship estimation method, and a causal relationship estimation program.
  • Patent Document 1 calculates the maximum mean discrepancy (MMD), which is the value that maximizes the difference between the kernel averages of bivariate time series data (also called “two-dimensional time series data”), and performs supervised learning. Discloses a device for clarifying causal relationships between variables.
  • MMD maximum mean discrepancy
  • Patent Document 1 has a problem in that it is not possible to estimate causal relationships between time-series data of an arbitrary number of variables (for example, three or more variables).
  • the present disclosure has been made to solve the above problems, and provides a causal relationship estimation device, an estimation method, and an estimation program that enable estimation of causal relationships between time-series data of an arbitrary number of variables. With the goal.
  • the causality estimation device of the present disclosure includes a data acquisition unit that acquires learning data including a set of time series data of a plurality of state variables and a set of time series data of a plurality of observed variables; Calculate a causal relationship parameter indicating a causal relationship between the time series data and the time series data of the plurality of observed variables, calculate a variance-covariance matrix of a Gaussian process from the learning data and the causal relationship parameter, a calculation unit that expresses the causality parameter using a multitask Gaussian process model; an optimization unit that calculates an optimization function based on the variance-covariance matrix and updates the causality parameter based on the optimization function; It is characterized by having the following.
  • Another causal relationship estimation device of the present disclosure reads a causal relationship parameter indicating a causal relationship for each layer between time series data of a plurality of state variables and time series data of a plurality of observation variables from a causal relationship parameter database.
  • a causal graph construction unit that rearranges the time-series data of the plurality of state variables and the time-series data of the plurality of observation variables for each layer based on the causal relationship parameter;
  • a data acquisition unit that acquires verification data including a set of data and a set of time series data of a plurality of observed variables; and a rearranged time series data of the plurality of state variables and a time series of the plurality of observed variables.
  • the present invention is characterized by comprising a causal graph verification unit that performs one or both of verification using the verification data of Granger causality and verification using the verification data of pseudo correlation on data.
  • FIG. 1 is a block diagram showing the configuration of a causal relationship estimation device according to Embodiment 1.
  • FIG. 2 is a diagram illustrating an example of the hardware configuration of the causal relationship estimation device in FIG. 1.
  • FIG. 2 is a block diagram showing the configuration of a preprocessing section of the data acquisition section in FIG. 1.
  • FIG. 2 is a diagram showing an example of time information before being dimensionally compressed by the preprocessing section of the data acquisition section of FIG. 1 and time information after being dimensionally compressed.
  • FIG. 2 is a diagram showing an example of angle information before being dimensionally expanded by a preprocessing unit of the data acquisition unit in FIG. 1 and angle information after being dimensionally expanded.
  • 2 is a diagram illustrating an example of processing performed by a calculation unit of the learning unit in FIG.
  • FIG. FIG. 2 is a diagram illustrating an example of processing performed by a time shift operator of the calculation unit of the learning unit in FIG. 1;
  • FIG. 2 is a diagram showing, in a table format, an example of initial values and numbers of causality parameters created by a calculation unit of a learning unit in FIG. 1;
  • 2 is a flowchart showing the operation of the causal relationship estimation device of FIG. 1.
  • FIG. FIG. 2 is a block diagram showing the configuration of a causal relationship estimation device according to a second embodiment.
  • 11 is a diagram illustrating an example of the hardware configuration of the causal relationship estimation device in FIG. 10.
  • FIG. 11 is a diagram showing a process of constructing a causal graph by a causal graph constructing unit of the constructing unit in FIG. 10.
  • FIG. 10 is a diagram showing a process of constructing a causal graph by a causal graph constructing unit of the constructing unit in FIG. 10.
  • 10A and 10B are diagrams illustrating an example of rearrangement processing of state information and observation information performed by the causal graph construction unit of the construction unit in FIG. 10.
  • FIG. (A) and (B) are diagrams illustrating an example of verification processing performed by the causal graph verification unit of the construction unit in FIG. 10.
  • (A) and (B) are diagrams showing other examples of the verification process performed by the causal graph verification unit of the construction unit in FIG. 10.
  • 11 is a flowchart showing the operation of the causality estimation device of FIG. 10.
  • 10A and 10B are diagrams showing a process of predicting unobserved time series data from relationships between time series data by the causality estimation device of FIG. 10.
  • FIG. 11 is a diagram illustrating an operation when the causal relationship estimation device of FIG. 10 uses a model that has learned change points or failures of sensor data by introducing a Change Point Kernel.
  • the causal relationship estimation device uses time-series data (X) regarding various economic indicators (for example, yen-dollar exchange rate, oil price, public investment, etc.) and time-series data (Y) on company stock prices. ) is given as a sample, the causal relationship between the time series data (X) and the time series data (Y) can be expressed as "X ⁇ Y", that is, "time series data ( This is a device for estimating that "X) is the cause and time series data (Y) is the result.”
  • Time series data (X) is called a state variable
  • time series data (Y) is called an observed variable.
  • a state variable is also called an explanatory variable or a latent variable.
  • Observed variables are also called objective variables, dependent variables, or explained variables.
  • both the time series data (X) and the time series data (Y) do not need to be one-dimensional time series data, but can be multidimensional time series data of two or more variables (i.e., multivariate time series data). It's okay.
  • the time series data (X) that is a state variable may be an observed variable of time series data (X') of another state variable.
  • the time series data (Y) that is an observation variable may be a state variable of time series data (Y') of another observation variable.
  • the causality estimation device is a machine learning device that generates a learned model (including causality parameters) from learning time series data.
  • the causality estimation device uses the causality parameters of the generated trained model (for example, the causality parameters generated by the causality estimation device according to the first embodiment) and the verification time series data. This is a device that constructs and outputs a causal graph.
  • the causality estimation device is a separate device from the causality estimation device according to the first embodiment, it may have the configuration of the causality estimation device according to the first embodiment.
  • the learning unit of the causal relationship estimation device models multidimensional time series data including correlations between data series and lag information between data series using, for example, a Gaussian process model.
  • a Gaussian process model known methods such as a Gauss Process Dynamical Model (GPDM) and a Multi Task Gauss Process (MTGP) model can be used.
  • GPDM is described in, for example, Non-Patent Document 1.
  • MTGP model is described in, for example, Non-Patent Document 2.
  • the causal relationship estimation device stores the causal relationship obtained through learning as a causal relationship parameter in a causal relationship parameter database (causal relationship parameter DB) of a storage device.
  • the causal relationship estimation device constructs a causal graph from the causal relationship parameters stored in the causal relationship parameter DB, and calculates the causal relationship using the time series data for verification stored in the time series data DB. Verify the graph and output the verified causal graph.
  • FIG. 1 is a block diagram showing the configuration of a causal relationship estimation device 100 according to the first embodiment.
  • the causal relationship estimation device 100 is a device that can implement the causal relationship estimation method (that is, the learning method) according to the first embodiment.
  • the causal relationship estimation device 100 includes a data acquisition section 10 and a learning section 20.
  • the causality estimation device 100 is, for example, a computer.
  • the causal relationship parameters obtained through learning are stored in the causal relationship parameter DB 26 in the storage section.
  • the causality parameter DB 26 may be stored in a storage unit (storage unit 103 in FIG.
  • causality estimation device 100 may be stored in an external device other than the causality estimation device 100 (for example, It may be stored in a storage unit of a server on a network that can communicate with the causal relationship estimation device 100.
  • the causality estimation device 100 acquires learning data including a set X (0:t) of time series data of a plurality of state variables and a set Y (0:t) of time series data of a plurality of observed variables, Calculate the causal relationship parameter ⁇ that indicates the causal relationship between the time series data of multiple state variables and the time series data of multiple observed variables, and calculate the variance-covariance matrix of the Gaussian process from the training data and the causal relationship parameter ⁇ .
  • the causality parameter ⁇ is a correlation coefficient that indicates the correlation between time series data of multiple state variables and time series data of multiple observed variables, and a correlation coefficient that indicates the correlation between time series data of multiple state variables and time series data of multiple observed variables.
  • the causal relationship estimation device 100 expresses the correlation as a linear correlation of an LMC (Linear Model of Corregion) kernel of the multi-task Gaussian process model.
  • the data acquisition unit 10 includes an input unit 14 that receives time series data from a time series data DB 12 and outputs the time series data to a preprocessing unit 16, and performs preprocessing on the time series data output from the input unit 14. It has a preprocessing unit 16 that outputs preprocessed time series data to the learning unit 20.
  • the time series data DB 12 may be stored in a storage unit (the storage unit 103 in FIG. 2 described later) as a part of the causality estimation device 100, but may be stored in an external device other than the causality estimation device 100 (for example, It may be stored in a storage unit of a server on a network that can communicate with the causal relationship estimation device 100.
  • the learning section 20 includes a calculation section 22 and an optimization section 24.
  • the calculation unit 22 initializes causality parameters and calculates a variance-covariance matrix.
  • the optimization unit 24 calculates an optimization function based on the variance-covariance matrix, and updates the causality parameter based on the optimization function. Details of each part will be described later.
  • FIG. 2 is a diagram showing an example of the hardware configuration of the causality estimation device 100.
  • the causality estimation device 100 includes a processor 101, a memory 102, and a storage unit 103 that is a nonvolatile storage device.
  • the causal relationship estimation device 100 may include a communication unit that communicates with other devices via a network.
  • the processor 101 is a CPU (Central Processing Unit) or the like.
  • the memory 102 is, for example, a volatile semiconductor memory such as RAM (Random Access Memory).
  • the storage unit 103 is a storage device such as a hard disk drive (HDD) or a solid state drive (SSD).
  • the storage unit 103 stores information (eg, various databases) and programs.
  • Each function of the causal relationship estimation device 100 is realized by a processing circuit.
  • the processing circuit may be dedicated hardware or may be the processor 101 that executes a program stored in the memory 102.
  • the processor 101 may be any one of a processing device, an arithmetic device, a microprocessor, a microcomputer, and a DSP (Digital Signal Processor).
  • the processing circuit may be, for example, a single circuit, a composite circuit, a programmed processor, an ASIC (Application Specific Integrated Circuit), an FPGA (Field-Programmable Gate Array), or Of these It is a combination of either.
  • ASIC Application Specific Integrated Circuit
  • FPGA Field-Programmable Gate Array
  • the causal relationship estimation program (that is, the learning program) according to the first embodiment is realized by software, firmware, or a combination of software and firmware.
  • Software and firmware are written as programs and stored in memory 102.
  • the processor 101 can realize the functions of each part shown in FIG. 1 by reading and executing a program stored in the memory 102.
  • the causal relationship estimation program is provided by downloading via a network or by a recording medium that records information such as an optical disk (i.e., a computer-readable storage medium), and is installed in the causal relationship estimation device 100. Ru.
  • the causal relationship estimation device 100 may be partially implemented using dedicated hardware and partially implemented using software or firmware. In this way, the processing circuit can implement the functions of each functional block shown in FIG. 1 using hardware, software, firmware, or any combination thereof.
  • Input section 14 The input unit 14 of the data acquisition unit 10 receives a set of time-series data of state variables X(0:t) (i.e., a plurality of time-series data) and a set of time-series data of observed variables Y( 0:t) (that is, a plurality of time series data) and outputs them to the preprocessing unit 16.
  • a set of time-series data of state variables X(0:t) i.e., a plurality of time-series data
  • Y( 0:t) that is, a plurality of time series data
  • the input unit 14 inputs time-series data x 1 (where Q is a positive integer) of state variables that are considered to have a cause-and-effect relationship from the time-series data DB 12. 0:t), x 2 (0:t),..., x Q (0:t) and time series data of D observed variables (D is a positive integer) y 1 (0:t) , y 2 (0:t), ..., y D (0:t), respectively.
  • time-series data of D observed variables are observed due to time-series data of Q state variables.
  • select the time series data of the angle ⁇ which is the state information as the time series data x q (0:t) of the state variable expand the dimension of the angle ⁇ which is the state information as described later, and create the time series data of sin ⁇ .
  • time series data of It may also be converted into two-dimensional time series data.
  • the input unit 14 inputs the time series data of the Q state variables selected from the time series data DB 12 and the time series data of the D observation variables to a length T (T is a positive integer). ) and a set of time-series data of state variables X(0:t) and a set of time-series data of observed variables Y(0:t) are obtained and passed to the preprocessing unit 16. Note that in X(0:t) and Y(0:t), 0:t in parentheses represents time series data from time 0 to time t. However, the time series data does not need to start with 0, and time series data having a length T at an arbitrary location may be selected from the time series data stored in the time series data DB 12.
  • x q (0:t) represents time-series data of state variables of state q.
  • y d (0:t) represents time-series data of the d-th dimension observed variable.
  • the preprocessing unit 16 of the data acquisition unit 10 acquires a set X (0:t) of time series data of state variables output from the input unit 14 and a set Y (0:t) of time series data of observed variables. Then, preprocessing is performed on these, and a set of preprocessed time series data of state variables X (0: t) and a set of preprocessed time series data of observed variables Y (0: t) are learned. It is output to the calculation section 22 of the section 20.
  • FIG. 3 is a block diagram showing the configuration of the preprocessing section 16 of the data acquisition section 10. As shown in FIG. 3, the preprocessing section 16 includes a dimension changing section 17 and a normalizing section 18.
  • the dimension change unit 17 performs dimension reduction processing (i.e. dimension compression processing) or dimension reduction processing appropriate for each state information on the time-series data x q (0:t) of state information whose dimensions need to be changed. Perform expansion processing.
  • dimension reduction processing i.e. dimension compression processing
  • dimension reduction processing i.e. dimension reduction processing
  • dimension reduction processing appropriate for each state information on the time-series data x q (0:t) of state information whose dimensions need to be changed.
  • Perform expansion processing When it is necessary to change the dimensions of the time series data x q (0:t) of the state information, for example, when the state information is periodic information such as time information or angle information, and when the dimensions of each side of a rectangle are For example, when the diagonal length is more effective than the length. Examples of dimension compression processing and dimension expansion processing are shown below. Note that it is also possible to perform the inverse processing of each of the following examples.
  • FIG. 4 is a diagram showing an example of time information before being dimensionally compressed by the dimension changing unit 17 of the preprocessing unit 16 (table on the left) and time information after being dimensionally compressed (table on the right).
  • the time information represents a periodicity of every 24 hours. Therefore, by integrating two pieces of information, ⁇ minutes'' and ⁇ hours,'' and compressing them into one piece of information, ⁇ hours and minutes,'' time information can be converted from 3D time series data to 2D time series data. Can be compressed.
  • FIG. 5 is a diagram showing an example of angle information before being dimensionally expanded by the dimension changing unit 17 of the preprocessing unit 16 and angle information after being dimensionally expanded.
  • the angle information is expanded from one piece of information, ⁇ angle ⁇ °,'' to a combination of two pieces of information, ⁇ sin ⁇ '' and ⁇ cos ⁇ .'' Data can be expanded to two-dimensional time series data.
  • the dimension changing unit 17 may use another dimension compression method (for example, principal component analysis, etc.) or another dimension expansion method.
  • the normalization unit 18 divides the set X (0:t) of time series data of state variables whose dimension has been changed and the set Y (0:t) of time series data of observed variables whose dimension has been changed into a set with an average of 0. Normalize so that the variance is 1.
  • the calculation unit 22 of the learning unit 20 receives from the preprocessing unit 16 a set of time series data of preprocessed state variables X(0:t) and a set of time series data of preprocessed observed variables Y(0:t). ), calculates the variance-covariance matrix K(X, X') of the Gaussian process, and outputs the variance-covariance matrix K(X, X') to the optimization unit 24.
  • FIG. 6 is a diagram showing an example of processing performed by the calculation unit 22 of the learning unit 20.
  • the correlation between time-series data of observed variables expressed in GPDM is expressed by the linear correlation of the LMC kernel of the MTGP model.
  • a specific example of the lag effect is that when the price of crude oil rises, the price of gasoline does not increase in tandem, but after a certain period of time (for example, until the next purchase is made), the price of gasoline increases. be.
  • another specific example of the lag effect is to raise gasoline prices first in anticipation of higher oil prices in the future.
  • the calculation unit 22 can take time series or periodicity into consideration. However, the calculation unit 22 may not consider the lag effect regarding time information.
  • X(t) represents a set of state variables ⁇ x 1 (t), x 2 (t), . . . , x Q (t) ⁇ .
  • u q (t) and v d (t) represent white Gaussian noise.
  • the state transition function f(x) and state function g(x) of the model are modeled by a Gaussian process, and using Gaussian process notation, the following equation (3) is generally used. It is expressed as (4).
  • Equation (3) represents the nonlinear time evolution of the state
  • Equation (4) represents the transformation from the state function to the observation function.
  • Kg and Kf represent Gram matrices.
  • the state function g q (x) in Equation (1) represents a Gaussian process model gp(0, k(x q , x q ′)) generated from the state information x q .
  • k q represents a positive definite kernel.
  • the positive definite kernel to be used one suitable for the data, such as an RBF kernel (Radial basis function kernel), is selected.
  • the RBF kernel is given by the following equation (5).
  • a d,q represents the linear correlation of the LMC of the conventional method, and represents the correlation coefficient from the state variable q to the observed variable d.
  • Equation (7) of Embodiment 1 is a time-shift operator newly introduced in the first embodiment, and L d,q represents a lag coefficient from state variable q to observation variable d.
  • F L represents a time-shift operator (lag operator) that shifts the state information of the state function to the future (or past) along the time axis, and is defined as in the following equation (8).
  • FIG. 7 is a diagram illustrating an example of processing performed by the time-shift operator of the calculation unit 22 of the learning unit 20.
  • FIG. 7 shows the operation of the time-shift operator FL .
  • the example in FIG. 7 shows a case where information is shifted into the future along the time axis t (L>0).
  • the direction of shift is the opposite direction to that in the case of FIG.
  • equation (7) When the observed data follows a Gaussian process, f is expressed by the multidimensional Gaussian distribution of equation (9) below.
  • K(X, X') is called a variance-covariance matrix (or Gram matrix), and is a matrix representing the similarity between state variables (that is, X ⁇ X').
  • the variance-covariance matrix K(X,X') is expressed as the following equations (11) and (12).
  • B q is called a coregionalization matrix, represents a linear transformation from a state function to an observation function, and is expressed as follows.
  • FIG. 8 is a diagram showing an example of the initial value and number of the causality parameter ⁇ in a table format. These causal relationship parameters are optimized by the optimization unit 24 of the learning unit 20.
  • the optimization unit 24 receives the variance-covariance matrix K(X, The process is performed and the optimized causal relationship parameter ⁇ is stored in the causal relationship parameter DB.
  • the marginal likelihood can be obtained by calculating as follows.
  • the probability that observation information is observed can be obtained from equation (13) as follows.
  • K ⁇ (X, X') represents the variance-covariance matrix K (X, X') calculated using the causality parameters.
  • N is the length of the feature vector X
  • D is the number of output dimensions of y.
  • the optimization unit 24 updates ⁇ so that the optimization function E in equation (15) is minimized. During optimization, if the causal relationship parameter ⁇ is updated, it is also necessary to update K ⁇ (X, X'), so the calculation unit 22 calculates K ⁇ (X, X').
  • the optimization unit 24 can use a known technique such as stochastic gradient descent during optimization. For example, L d,q can be optimized using a grid search, and the remaining causality parameters can be optimized using stochastic gradient descent.
  • each of the state variables and observation variables can have one or more dimensions, so it is possible to estimate the causal relationship between time series data of any number of variables, and it is possible to It is possible to estimate the causal relationship of time series data of variables (Q+D).
  • Embodiment 1 it becomes possible to understand the causal relationships between multiple state variables and multiple observed variables, and it is possible to estimate not only the causal relationships of bivariate time series data but also time series data of three or more variables. It is possible to infer causal relationships in data.
  • FIG. 9 is a flowchart showing the operation (ie, learning method) of the causal relationship estimation device 100 according to the first embodiment.
  • the input unit 14 obtains a set of time-series data of state variables X(0:t) and a set of time-series data of observed variables Y(0:t) from the time-series data DB 12.
  • step S102 the preprocessing unit 16 performs a dimension change (i.e., dimension compression) of the state variable x q (0:t) that requires dimension change in the state variable time series data set X (0:t). or dimension expansion).
  • a dimension change i.e., dimension compression
  • step S103 the preprocessing unit 16 normalizes the set of time-series data of state variables X(0:t) and the set of time-series data of observed variables Y(0:t).
  • step S104 the calculation unit 22 sets the following causal relationship parameter ⁇ to an initial value.
  • step S105 the calculation unit 22 calculates the variance-covariance matrix K ⁇ (X, ⁇ ) Calculate.
  • step S106 the optimization unit 24 calculates the optimization function E of equation (15).
  • step S107 the optimization unit 24 optimizes (that is, updates) the causality parameter ⁇ so that the optimization function E is minimized.
  • Steps S105 to S107 are repeated until the optimization function E converges.
  • the correlation coefficients a d, q and the lag coefficients L d, q, which are causal parameters are stored in the causal parameter DB 26 .
  • Embodiment 1 Effects According to Embodiment 1, it becomes possible to understand the causal relationships between multiple state variables and multiple observed variables, and it is possible to estimate the causal relationships of time-series data of any number of variables. Can be done.
  • Embodiment 2 ⁇ 2-1>> Configuration ⁇ 2-1-1>>
  • Causal relationship estimation device 200 The causal relationship estimation device 200 reads a causal relationship parameter ⁇ indicating a causal relationship for each layer between the time series data of a plurality of state variables and the time series data of a plurality of observed variables from the causal relationship parameter DB, and calculates the causal relationship.
  • the time series data of multiple state variables and the time series data of multiple observation variables are rearranged for each layer, and the set X(0:t) of time series data of multiple state variables and the time series data of multiple observed variables are Obtain verification data including a set Y(0:t) of time-series data of observed variables, and apply Granger causality to the rearranged time-series data of multiple state variables and time-series data of multiple observed variables.
  • One or both of the verification data using the verification data and the verification using the verification data of pseudo correlation are performed.
  • the causal relationship estimation device 200 is a causal graph construction device that has a causal graph construction function that constructs a causal graph.
  • a causal graph is a graph structure of a data item list with a "cause ⁇ effect" relationship based on the causal relationship information obtained in the first embodiment.
  • the state functions and observation functions that are elements of the constructed causal graph are rearranged based on the correlation coefficients a d, q and the lag coefficients L d, q stored in the causal relationship parameter DB.
  • the constructed causal graph is verified using Granger causality and pseudo-correlation.
  • the verified causal graph is output.
  • the causal relationship estimation device 200 constructs a causal graph using correlation coefficients a d, q and lag coefficients L d, q, which are causal relationship parameters.
  • the correlation coefficients a d, q and the lag coefficients L d, q, which are causal relationship parameters are, for example, causal relationship parameters of the learned model created by the causal relationship estimation device 100 according to the first embodiment.
  • the causality estimation device 200 is, for example, a computer.
  • Causal relationship estimation device 200 may be the same computer as the computer that constitutes causal relationship estimation device 100 according to Embodiment 1, or may be a different computer.
  • FIG. 10 is a block diagram showing the configuration of a causal relationship estimation device 200 according to the second embodiment.
  • the causal relationship estimation device 200 is a device that can perform the causal relationship estimation (that is, the causal graph construction method) according to the second embodiment.
  • the causality estimation device 200 includes a construction section 30, a data acquisition section 40, and an output section 90.
  • the construction unit 30 includes a causal graph construction unit 32 and a causal graph verification unit 34.
  • the construction unit 30 may include a causality parameter DB 80 that provides causality parameters to the causal graph construction unit 32.
  • the data acquisition section 40 includes an input section 44 and a preprocessing section 46.
  • the data acquisition unit 40 may include a time series data DB 42 that provides time series data to the input unit 44.
  • the output unit 90 outputs the causal graph constructed by the construction unit 30.
  • the time series data DB 42 and the causal relationship parameter DB 80 may be stored in a storage unit as part of the causal relationship estimation device 200 (storage unit 203 in FIG. It may be stored in a storage unit of an external device (for example, a server on a network that can communicate with the causal relationship estimation device 200).
  • FIG. 11 is a diagram showing an example of the hardware configuration of the causal relationship estimation device 200.
  • the causality estimation device 200 includes a processor 201, a memory 202, and a storage unit 203 that is a nonvolatile storage device.
  • the causal relationship estimation device 200 may include an interface with an external device, a communication unit that communicates with other devices via a network, and the like.
  • Processor 201 is a CPU or the like.
  • Memory 202 is, for example, a volatile semiconductor memory such as RAM.
  • the storage unit 203 is a storage device such as an HDD or an SSD.
  • the storage unit 203 stores information (eg, various databases) and programs.
  • Each function of the causality estimation device 200 is realized by a processing circuit.
  • the processing circuit may be dedicated hardware or may be the processor 201 that executes a program stored in the memory 202.
  • the estimation program (that is, the causal graph construction program) according to the second embodiment is realized by software, firmware, or a combination of software and firmware.
  • Software and firmware are written as programs and stored in memory 202.
  • the processor 201 can realize the functions of each section shown in FIG. 10 by reading and executing a program stored in the memory 202.
  • the program is installed in the causality estimation device 200 by downloading via a network or from a recording medium that records information, such as an optical disc.
  • the data acquisition unit 40 has the same functions as the data acquisition unit 10 of the first embodiment. However, the input unit 44 acquires data for verifying a causal graph using Granger causality and pseudo-correlation, which will be described later.
  • the causal graph construction unit 32 of the construction unit 30 uses the causal relationship parameters stored in the causal relationship parameter DB 80. is acquired, a causal graph is constructed using the causal relationship parameters, and the constructed causal graph is output to the causal graph verification unit 34.
  • FIG. 12 is a diagram showing the process of constructing a causal graph by the causal graph constructing unit 32 of the constructing unit 30.
  • y observation information
  • the correlation coefficient between x (state information) and y observation information
  • FIGS. 13A and 13B are diagrams illustrating an example of state information and observation information rearrangement processing performed by the causal graph construction unit 32 of the construction unit 30. As shown in FIGS. 13A and 13B, the state information and observation information are connected by an arrow pointing from the state information to the observation information.
  • FIG. 13(A) if the observation information has a lag (advancement) than the state function, the direction of the arrow is reversed.
  • FIG. 13(B) if the correlation coefficient is low, consider the possibility that it is a confounding factor or an intermediate factor, and rearrange the factor so as to move it to a higher position in the causal graph.
  • the causality parameter is calculated by the causality estimation device 100 according to the first embodiment using the state variables and observation variables after the rearrangement. Ask again.
  • FIGS. 14A and 14B are diagrams illustrating an example of verification processing performed by the causal graph verification unit 34 of the construction unit 30.
  • Granger causality is applied when observed variables y 1 and y 2 are predicted by state variables x 1 , x 2 , and x 3 including state variable x 2 (case 1).
  • the prediction accuracy is lower when the observed variables y 1 and y 2 are predicted by the state variables x 1 and x 3 (case 2) by deleting the state variable x 2. If the value decreases, the state variable x 2 is considered to have Granger causality with respect to the observed variables y 1 and y 2 .
  • the observed variables y 1 and y 2 were predicted by the state variables x 1 , x 2 , and x 3 (case 1)
  • the observed variables y 1 and y 2 were predicted by the state variables x 1 and x 3 .
  • the prediction accuracy increases in case 2
  • future predictions are made from the set of time-series data of state variables X(0:t) and the set of time-series data of observed variables Y(0:t) used for estimating the causal relationship in the causality estimation device 100.
  • a set of time-series data of state variables X (t+1:t+ ⁇ t) and a set of time-series data of future observed variables Y(t+1:t+ ⁇ t) may be used as prediction verification data (i.e., test data).
  • the prediction error can be evaluated using, for example, the following RMSE (Root Mean Squared Error).
  • FIGS. 15A and 15B are diagrams showing other examples of the verification process performed by the causal graph verification unit 34 of the construction unit 30. From Figures 15 (A) and (B), we can see how much a factor is affected in an unsteady state (a factor undergoes a sudden change due to an external factor (impulse response)) (or how much it remains constant without any change). It is also possible to verify this by estimating .
  • state variable x 2 affects both observed variables y 1 and y 2 , but state variables x 1 and x 3 are different from observed variables y 1 and y 2 , respectively.
  • x 2 does not change, such as taking a fixed value, x 2 does not change and does not affect observed variables y 1 and y 2 , so there is no correlation between observed variables y 1 and y 2 . It is thought that it will disappear.
  • the correlation between observed variables y 1 and y 2 remains even if the influence of x 2 disappears, or the correlation between x 1 and observed variables y 1 and between x 2 and y 2 is high. If a correlation is shown, the state variable x 2 is rearranged or deleted from the causal graph.
  • a causal graph By constructing a causal graph, it becomes possible to predict sensor values in places that cannot be directly measured or search for the cause of anomalies. Possible applications include predicting road damage based on traffic volume or weather information, predicting river water levels based on nearby rainfall and water level data, and predicting power demand using weather or economic status data in surrounding areas.
  • FIG. 16 is a flowchart showing the operation (ie, inference operation) of the causal relationship estimation device 200.
  • the construction unit 30 extracts the information for each layer from the causal relationship parameter DB. get.
  • step S202 the causal graph construction unit 32 arranges the state information and observation information in the order of "state information ⁇ observation information" for each layer.
  • step S203 the causal graph construction unit 32 corrects the causal direction (that is, the direction of the arrow).
  • step S204 the causal graph construction unit 32 rearranges the state information and observation information of the causal graph.
  • step S205 the input unit 44 acquires Granger causality verification data from the time series data DB 42.
  • step S206 the preprocessing unit 46 changes the dimension of the verification data.
  • step S207 the preprocessing unit 46 normalizes the verification data.
  • step S208 the causal graph verification unit 34 verifies the causal graph using Granger causality.
  • step S209 the input unit 44 acquires data for verification by pseudo correlation from the time series data DB 42.
  • step S210 the preprocessing unit 46 changes the dimension of the verification data.
  • step S211 the preprocessing unit 46 normalizes the verification data.
  • step S212 the causal graph verification unit 34 verifies the causal relationship using pseudo correlation.
  • the causal graph verified by the causal graph verification unit 34 is output to the output unit 90.
  • Embodiment 2 also includes, for example, predicting road surface damage from traffic volume or weather information, predicting river water level from nearby rainfall and water level data, predicting power demand using weather or economic status data of surrounding areas, It is applicable to
  • FIGS. 17A and 17B show a learning unit and an inference unit that perform a process of predicting unobserved time series data from the relationships between time series data by the causality estimation device 200.
  • the causal relationship estimation device 100 generates a learned model by performing multi-task learning on past congestion information of A1, A2, and A3 stations shown in FIG. 17(A).
  • the causal relationship estimation device 100 knows only the congestion information of A1 station and A2 station among A1 station, A2 station, and A3 station, and does not know the congestion information of A3 station (case of FIG. 17(B))
  • the inference unit uses the trained model generated by multi-task learning, it is possible for the inference unit to predict (infer) the degree of congestion at A3 station from the correlation or lag information of mutual congestion degrees.
  • FIG. 18 is a diagram illustrating an operation when changing points or failures in sensor data are learned by introducing a Change Point Kernel into the causal relationship estimation device 100. Failure prediction can be performed by using a trained model generated by learning sensor data change points or failures.
  • causal relationship estimation device 100 causal relationship estimation device, 200 causal relationship estimation device (causal graph construction device), 10 data acquisition unit, 12, 42 time series data DB, 14 input unit, 16 preprocessing unit, 17 dimension change unit, 18 normalization unit, 20 learning unit, 22 calculation unit, 24 optimization unit, 26, 80 causality parameter DB, 30 construction unit, 32 causal graph construction unit, 34 causal graph verification unit, 40 data acquisition unit, 44 input unit, 46 preprocessing unit , 90 output section.

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

因果関係推定装置（１００）は、複数の状態変数の時系列データの集合（Ｘ（０：ｔ））と複数の観測変数の時系列データの集合（Ｙ（０：ｔ））とを含む学習用データを取得するデータ取得部（１０）と、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータ（θ）を計算し、学習用データと因果関係パラメータ（θ）とからガウス過程の分散共分散行列（Ｋ（Ｘ，Ｘ´））を計算し、因果関係パラメータ（θ）をマルチタスクガウス過程モデルで表現する計算部（２２）と、分散共分散行列に基づいて最適化関数を計算し、最適化関数に基づいて因果関係パラメータ（θ）を更新する最適化部（２４）とを有する。

Description

因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム
　本開示は、因果関係推定装置、因果関係推定方法、及び因果関係推定プログラムに関する。
　例えば、特許文献１は、２変量時系列データ（「２次元時系列データ」とも呼ばれる。）のカーネル平均の差が最大になる値であるｍａｘｉｍｕｍ　ｍｅａｎ　ｄｉｓｃｒｅｐａｎｃｙ（ＭＭＤ）を計算し、教師あり学習によって変量間の因果関係を明らかにする装置を開示している。
特開２０１７－２２８２５６号公報
Ｊａｃｋ　Ｍ．　Ｗａｎｇ，　Ｄａｖｉｄ　Ｊ．　Ｆｌｅｅｔ，　ａｎｄ　Ａａｒｏｎ　Ｈｅｒｔｚｍａｎｎ，　"Ｇａｕｓｓｉａｎ　Ｐｒｏｃｅｓｓ　Ｄｙｎａｍｉｃａｌ　Ｍｏｄｅｌｓ"，　ＮＩＰＳ´０５：　Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ｔｈｅ　１８ｔｈ　Ｉｎｔｅｒｎａｔｉｏｎａｌ　Ｃｏｎｆｅｒｅｎｃｅ　ｏｎ　Ｎｅｕｒａｌ　Ｉｎｆｏｒｍａｔｉｏｎ　Ｐｒｏｃｅｓｓｉｎｇ　Ｓｙｓｔｅｍｓ，　Ｄｅｃｅｍｂｅｒ　２００５，　ｐｐ．１４４１－１４４８ Ｅｄｗｉｎ　Ｖ．　Ｂｏｎｉｌｌａ，　Ｋｉａｎ　Ｍｉｎｇ　Ａ．　Ｃｈａｉ，　Ｃｈｒｉｓｔｏｐｈｅｒ　Ｋ．Ｉ．　Ｗｉｌｌｉａｍｓ，　"Ｍｕｌｔｉ－ｔａｓｋ　Ｇａｕｓｓｉａｎ　Ｐｒｏｃｅｓｓ　Ｐｒｅｄｉｃｔｉｏｎ"，　Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ｔｈｅ　Ａｄｖａｎｃｅｓ　ｉｎ　Ｎｅｕｒａｌ　Ｉｎｆｏｒｍａｔｉｏｎ　Ｐｒｏｃｅｓｓｉｎｇ　Ｓｙｓｔｅｍｓ　２０，　（２００８）
　しかしながら、特許文献１の装置では、任意数の変量（例えば、３変量以上）の時系列データ間の因果関係を推定することができないという課題がある。
　本開示は、上記課題を解決するためになされたものであり、任意数の変量の時系列データ間の因果関係の推定を可能にする因果関係推定装置、推定方法、及び推定プログラムを提供することを目的とする。
　本開示の因果関係推定装置は、複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するデータ取得部と、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現する計算部と、前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新する最適化部と、を有することを特徴とする。
　本開示の他の因果関係推定装置は、因果関係パラメータデータベースから、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の階層ごとの因果関係を示す因果関係パラメータを読み出し、前記因果関係パラメータに基づいて、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとを前記階層ごとに配置換えする因果グラフ構築部と、複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む検証用データを取得するデータ取得部と、配置換えされた前記複数の状態変数の時系列データと前記複数の観測変数の時系列データに対し、グレンジャー因果の前記検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う因果グラフ検証部とを有することを特徴とする。
　本開示によれば、任意数の変量の時系列データ間の因果関係の推定を可能にすることができる。
実施の形態１に係る因果関係推定装置の構成を示すブロック図である。 図１の因果関係推定装置のハードウェア構成の例を示す図である。 図１のデータ取得部の前処理部の構成を示すブロック図である。 図１のデータ取得部の前処理部によって次元圧縮される前の時刻情報と次元圧縮された後の時刻情報との例を示す図である。 図１のデータ取得部の前処理部によって次元拡張される前の角度情報と次元拡張された後の角度情報との例を示す図である。 図１の学習部の計算部によって行われる処理の例を示す図である。 図１の学習部の計算部の時間シフトオペレータによって行われる処理の例を示す図である。 図１の学習部の計算部で作成された因果関係パラメータの初期値と個数の例を表形式で示す図である。 図１の因果関係推定装置の動作を示すフローチャートである。 実施の形態２に係る因果関係推定装置の構成を示すブロック図である。 図１０の因果関係推定装置のハードウェア構成の例を示す図である。 図１０の構築部の因果グラフ構築部による因果グラフの構築処理を示す図である。 （Ａ）及び（Ｂ）は、図１０の構築部の因果グラフ構築部によって行われる状態情報及び観測情報の配置替え処理の例を示す図である。 （Ａ）及び（Ｂ）は、図１０の構築部の因果グラフ検証部によって行われる検証処理の例を示す図である。 （Ａ）及び（Ｂ）は、図１０の構築部の因果グラフ検証部によって行われる検証処理の他の例を示す図である。 図１０の因果関係推定装置の動作を示すフローチャートである。 （Ａ）及び（Ｂ）は、図１０の因果関係推定装置によって時系列データ間の関係性から観測されていない時系列データを予測する処理を示す図である。 図１０の因果関係推定装置が、Ｃｈａｎｇｅ　Ｐｏｉｎｔ　Ｋｅｒｎｅｌを導入することでセンサデータの変化点又は故障を学習したモデルを使用する場合の動作を示す図である。
　以下に、実施の形態に係る因果関係推定装置、因果関係推定方法、及び因果関係推定プログラムを、図面を参照しながら説明する。以下の実施の形態は、例にすぎず、実施の形態を適宜組み合わせること及び各実施の形態を適宜変更することが可能である。
　実施の形態に係る因果関係推定装置は、例えば、各種経済指標（例えば、円ドルの為替レート、石油価格、公共投資、など）に関する時系列データ（Ｘ）と企業の株価の時系列データ（Ｙ）とからなる多次元時系列データがサンプルとして与えられた場合に、時系列データ（Ｘ）と時系列データ（Ｙ）との間の因果関係を、「Ｘ→Ｙ」すなわち「時系列データ（Ｘ）が原因であり、時系列データ（Ｙ）が結果である。」というように推定するための装置である。時系列データ（Ｘ）は状態変数と呼ばれ、時系列データ（Ｙ）は観測変数と呼ばれる。また、状態変数は、説明変数又は潜在変数とも呼ばれる。観測変数は、目的変数、従属変数、又は被説明変数とも呼ばれる。また、時系列データ（Ｘ）及び時系列データ（Ｙ）のいずれも、１次元時系列データである必要はなく、２変量以上の多次元時系列データ（すなわち、多変量時系列データ）であってもよい。また、状態変数である時系列データ（Ｘ）は、別の状態変数の時系列データ（Ｘ´）の観測変数であってもよい。また、観測変数である時系列データ（Ｙ）は、別の観測変数の時系列データ（Ｙ´）の状態変数であってもよい。
　実施の形態１に係る因果関係推定装置は、学習用時系列データから学習済みモデル（因果関係パラメータを含む）を生成する機械学習装置である。実施の形態２に係る因果関係推定装置は、生成された学習済みモデルの因果関係パラメータ（例えば、実施の形態１に係る因果関係推定装置で生成された因果関係パラメータ）と検証用時系列データとから因果グラフを構築して出力する装置である。実施の形態２に係る因果関係推定装置は、実施の形態１に係る因果関係推定装置と別個の装置であるが、実施の形態１に係る因果関係推定装置の構成を有していてもよい。
　実施の形態１に係る因果関係推定装置の学習部は、例えば、ガウス過程モデルを用いて、データ系列間の相関とデータ系列間のラグ情報とを含む多次元時系列データをモデリングする。ガウス過程モデルとしては、ガウス過程動的力学モデル（Ｇａｕｓｓ　Ｐｒｏｃｅｓｓ　Ｄｙｎａｍｉｃａｌ　Ｍｏｄｅｌ：ＧＰＤＭ）及びマルチタスクガウス過程（Ｍｕｌｔｉ　Ｔａｓｋ　Ｇａｕｓｓ　Ｐｒｏｃｅｓｓ：ＭＴＧＰ）モデル、などの公知の方法を用いることができる。ＧＰＤＭは、例えば、非特許文献１に記載されている。ＭＴＧＰモデルは、例えば、非特許文献２に記載されている。
　ＧＰＤＭを用いることで、多次元時系列データの非線形な状態の時間発展を表現することができる。また、ＭＴＧＰモデルを用いることで、複数の状態関数と複数の観測関数との間の相関及び複数の状態関数と複数の観測関数との間のラグを、カーネル関数の因果関係パラメータとして表すことができる。実施の形態１に係る因果関係推定装置は、学習によって得られた因果関係を、記憶装置の因果関係パラメータデータベース（因果関係パラメータＤＢ）に因果関係パラメータとして保存する。実施の形態２に係る因果関係推定装置は、因果関係パラメータＤＢに保存されている因果関係パラメータから因果グラフを構築し、時系列データＤＢに保存されている検証用の時系列データを用いて因果グラフを検証して、検証済みの因果グラフを出力する。
《１》実施の形態１
《１－１》構成
《１－１－１》因果関係推定装置１００
　図１は、実施の形態１に係る因果関係推定装置１００の構成を示すブロック図である。因果関係推定装置１００は、実施の形態１に係る因果関係推定方法（すなわち、学習方法）を実施することができる装置である。図１に示されるように、因果関係推定装置１００は、データ取得部１０と、学習部２０とを有してる。因果関係推定装置１００は、例えば、コンピュータである。学習によって得られた因果関係パラメータは、記憶部の因果関係パラメータＤＢ２６に記憶される。因果関係パラメータＤＢ２６は、因果関係推定装置１００の一部としての記憶部（後述の図２の記憶部１０３）に格納されてもよいが、因果関係推定装置１００とは別の外部装置（例えば、因果関係推定装置１００と通信可能なネットワーク上のサーバ）の記憶部に格納されてもよい。
　因果関係推定装置１００は、複数の状態変数の時系列データの集合Ｘ（０：ｔ）と複数の観測変数の時系列データの集合Ｙ（０：ｔ）とを含む学習用データを取得し、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータθを計算し、学習用データと因果関係パラメータθとからガウス過程の分散共分散行列Ｋ（Ｘ，Ｘ´）を計算し、因果関係パラメータθをマルチタスクガウス過程モデルで表現し、分散共分散行列に基づいて最適化関数を計算し、最適化関数に基づいて因果関係パラメータθを更新する。因果関係パラメータθは、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の相関を示す相関係数と、複数の状態変数の時系列データと複数の観測変数の時系列データとの間のラグを示すラグ係数とを含み、因果関係推定装置１００は、相関を、前記マルチタスクガウス過程モデルのＬＭＣ（Ｌｉｎｅａｒ　Ｍｏｄｅｌ　ｏｆ　Ｃｏｒｅｇｉｏｎ）カーネルの線形相関で表現する。
　データ取得部１０は、時系列データＤＢ１２から時系列データが入力され、前処理部１６に時系列データを出力する入力部１４と、入力部１４から出力された時系列データに前処理を施して学習部２０に前処理済みの時系列データを出力する前処理部１６とを有している。時系列データＤＢ１２は、因果関係推定装置１００の一部としての記憶部（後述の図２の記憶部１０３）に格納されてもよいが、因果関係推定装置１００とは別の外部装置（例えば、因果関係推定装置１００と通信可能なネットワーク上のサーバ）の記憶部に格納されてもよい。
　学習部２０は、計算部２２と、最適化部２４とを有している。計算部２２は、因果関係パラメータを初期化し、分散共分散行列を計算する。最適化部２４は、分散共分散行列に基づいて最適化関数を計算し、最適化関数に基づいて因果関係パラメータを更新する。各部の詳細については後述する。
　図２は、因果関係推定装置１００のハードウェア構成の例を示す図である。因果関係推定装置１００は、プロセッサ１０１と、メモリ１０２と、不揮発性の記憶装置である記憶部１０３とを有している。因果関係推定装置１００は、ネットワークを介して他の装置と通信を行う通信部を備えてもよい。プロセッサ１０１は、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）などである。メモリ１０２は、例えば、ＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）などの、揮発性の半導体メモリである。記憶部１０３は、ハードディスクドライブ（ＨＤＤ）又はソリッドステートドライブ（ＳＳＤ）などの記憶装置である。記憶部１０３は、情報（例えば、各種のデータベース）及びプログラムを記憶する。
　因果関係推定装置１００の各機能は、処理回路により実現される。処理回路は、専用のハードウェアであってもよいし、メモリ１０２に格納されるプログラムを実行するプロセッサ１０１であってもよい。プロセッサ１０１は、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、及びＤＳＰ（Ｄｉｇｉｔａｌ　Ｓｉｇｎａｌ　Ｐｒｏｃｅｓｓｏｒ）のいずれであってもよい。
　処理回路が専用のハードウェアである場合、処理回路は、例えば、単一回路、複合回路、プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、又はこれらのうちのいずれかを組み合わせたものである。
　処理回路がプロセッサ１０１である場合、実施の形態１に係る因果関係推定プログラム（すなわち、学習プログラム）は、ソフトウェア、ファームウェア、又はソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェア及びファームウェアは、プログラムとして記述され、メモリ１０２に格納される。プロセッサ１０１は、メモリ１０２に記憶されたプログラムを読み出して実行することにより、図１に示される各部の機能を実現することができる。因果関係推定プログラムは、ネットワークを介してのダウンロードにより、又は、光ディスクなどのような情報を記録する記録媒体（すなわち、コンピュータで読取可能な記憶媒体）によって提供され、因果関係推定装置１００にインストールされる。なお、因果関係推定装置１００は、一部を専用のハードウェアで実現し、一部をソフトウェア又はファームウェアで実現するようにしてもよい。このように、処理回路は、ハードウェア、ソフトウェア、ファームウェア、又はこれらのうちのいずれかの組み合わせによって、図１に示される各機能ブロックの機能を実現することができる。
《１－１－２》入力部１４
　データ取得部１０の入力部１４は、時系列データＤＢ１２から状態変数の時系列データの集合Ｘ（０：ｔ）（すなわち、複数個の時系列データ）及び観測変数の時系列データの集合Ｙ（０：ｔ）（すなわち、複数個の時系列データ）を取得し、これらを前処理部１６に出力する。
　具体的に言えば、入力部１４は、時系列データＤＢ１２から、原因と結果の関係にあると考えられる、Ｑ個（Ｑは正の整数である。）の状態変数の時系列データｘ１（０：ｔ），ｘ２（０：ｔ），…，ｘＱ（０：ｔ）と、Ｄ個（Ｄは正の整数である。）の観測変数の時系列データｙ１（０：ｔ），ｙ２（０：ｔ），…，ｙＤ（０：ｔ）とを、それぞれ選択する。このことは、Ｑ個の状態変数の時系列データが原因となって、Ｄ個の観測変数の時系列データが観測されたことを意味する。なお、観測変数の時系列データの各々ｙｄ（０：ｔ）（ｄ＝１，２，…，Ｄ）は、１次元時系列データである必要があるが、状態変数の時系列データの各々ｘｑ（０：ｔ）（ｑ＝１，２，…，Ｑ）は、１次元時系列データである必要はなく、多次元時系列データであってもよい。例えば、状態変数の時系列データｘｑ（０：ｔ）として状態情報である角度αの時系列データを選択し、後述するように状態情報である角度αを次元拡張してｓｉｎαの時系列データとｃｏｓαの時系列データとした場合、ｓｉｎαの時系列データとｃｏｓαの時系列データとのそれぞれを１次元時系列データの状態変数とせずに、これらのセットである（ｓｉｎα，ｃｏｓα）を１つの２次元時系列データに変換してもよい。
　続いて、入力部１４では、時系列データＤＢ１２から選択したＱ個の状態変数の時系列データと、Ｄ個の観測変数の時系列データとから、長さＴ（Ｔは正の整数である。）の状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）を取得し、これらを前処理部１６に渡す。なお、Ｘ（０：ｔ）及びＹ（０：ｔ）において、括弧内の０：ｔは、時刻０から時刻ｔまでの時系列データを表す。ただし、時系列データは、０で始まる必要はなく、時系列データＤＢ１２に格納されている時系列データから任意の個所の長さＴの時系列データを選択してもよい。
　Ｘ（０：ｔ）は、状態変数の時系列データの集合を表す。すなわち、
Ｘ（０：ｔ）＝｛ｘ１（０：ｔ），ｘ２（０：ｔ），…，ｘＱ（０：ｔ）｝である。
　Ｙ（０：ｔ）は、観測変数の時系列データの集合を表す。すなわち、
Ｙ（０：ｔ）＝｛ｙ１（０：ｔ），ｙ２（０：ｔ），…，ｙＤ（０：ｔ）｝である。
　なお、ｘｑ（０：ｔ）は、状態ｑの状態変数の時系列データを表す。また、ｙｄ（０：ｔ）は、ｄ次元目の観測変数の時系列データを表す。
《１－１－３》前処理部１６
　データ取得部１０の前処理部１６は、入力部１４から出力された状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）とを取得し、これらに前処理を施し、前処理済みの状態変数の時系列データの集合Ｘ（０：ｔ）と前処理済みの観測変数の時系列データの集合Ｙ（０：ｔ）とを、学習部２０の計算部２２に出力する。
　図３は、データ取得部１０の前処理部１６の構成を示すブロック図である。図３に示されるように、前処理部１６は、次元変更部１７と、正規化部１８とを有している。
　次元変更部１７は、次元の変更が必要な状態情報の時系列データｘｑ（０：ｔ）に対して、各々の状態情報に対して適切な次元削減処理（すなわち、次元圧縮処理）又は次元拡張処理を行う。状態情報の時系列データｘｑ（０：ｔ）の次元の変更が必要な場合は、例えば、状態情報が時刻情報又は角度情報などの周期性のある情報である場合、及び長方形の各辺の長さよりも対角線の長さのほうが有効である場合、などである。以下に、次元圧縮処理と次元拡張処理の例を示す。なお、以下の例のそれぞれの逆の処理を行うことも可能である。
　図４は、前処理部１６の次元変更部１７によって次元圧縮される前の時刻情報（左側の表）と次元圧縮された後の時刻情報（右側の表）との例を示す図である。図４に示されるように、時刻情報は、観測データが１日周期の場合には、２４時間ごとの周期性を表す。このため、２つの情報である「分」と「時」とを統合して１つの情報である「時分」に圧縮することで、時刻情報を３次元時系列データから２次元時系列データに圧縮することができる。
　図５は、前処理部１６の次元変更部１７によって次元拡張される前の角度情報と次元拡張された後の角度情報との例を示す図である。図５に示されるように、角度情報は、１つの情報である「角度α°」から２つの情報の組合せである「ｓｉｎα」と「ｃｏｓα」に拡張することで、角度情報を１次元時系列データから２次元時系列データに拡張することができる。なお、次元変更部１７は、別の次元圧縮手法（例えば、主成分分析等）又は別の次元拡張方法を用いてもよい。
　正規化部１８は、次元変更された状態変数の時系列データの集合Ｘ（０：ｔ）と次元変更された観測変数の時系列データの集合Ｙ（０：ｔ）とを、平均が０であり分散が１であるように正規化する。
《１－１－４》計算部２２
　学習部２０の計算部２２は、前処理部１６から前処理済の状態変数の時系列データの集合Ｘ（０：ｔ）と前処理済の観測変数の時系列データの集合Ｙ（０：ｔ）を受け取り、ガウス過程の分散共分散行列Ｋ（Ｘ，Ｘ´）の計算を行い、最適化部２４に分散共分散行列Ｋ（Ｘ，Ｘ´）を出力する。
　図６は、学習部２０の計算部２２によって行われる処理の例を示す図である。図６の例では、ＧＰＤＭで表される観測変数の時系列データ間の相関を、ＭＴＧＰモデルのＬＭＣカーネルの線形相関で表現する。
　図６においては、ｔｉｍｅ－ｓｈｉｆｔ　ｏｐｅｒａｔｏｒ（時間シフトオペレータ）（「Ｌａｇ　Ｏｐｅｒａｔｏｒ（ラグオペレータ）」とも呼ばれる。）ＦＬを導入することで、状態関数ｇ１（ｘ１（ｔ））～ｇＱ（ｘＱ（ｔ））が遅れて（又は、進んで）観測関数ｆ１（ｘ１（ｔ））～ｆＤ（ｘＤ（ｔ））に影響を及ぼすラグ効果を表現することができる。ラグ効果とは、原因が結果に同時刻に影響を及ぼすのではなく、「原因が遅れて結果に影響を及ぼす」又は「原因が速く結果に影響を及ぼす」ことを表す。ラグ効果の具体例は、原油価格が高くなって、ガソリン価格が連動して高くなるわけではなく、ある期間（例えば、次の仕入れが行われるまでの期間）を経てガソリン価格が高くなることである。あるいは、ラグ効果の他の具体例は、将来、原油価格が高くなることを見越して、先にガソリン価格を値上げすることである。
　状態情報の一つに時刻情報を加えることで、計算部２２は、時系列性又は周期性を考慮することができる。ただし、計算部２２は、時刻情報に関してのラグ効果を考慮しないことも可能である。
　以下、実施の形態１における分散共分散行列Ｋ（Ｘ，Ｘ´）の計算方法について説明する。ＧＰＤＭでは、状態方程式は式（１）のように定義され、観測方程式は式（２）のように定義される。
Figure JPOXMLDOC01-appb-M000001
　ここで、Ｘ（ｔ）は、状態変数の集合｛ｘ１（ｔ），ｘ２（ｔ），…，ｘＱ（ｔ）｝を表す。ｕｑ（ｔ）及びｖｄ（ｔ）は、ホワイトガウスノイズを表す。
　ガウス過程状態空間モデルでは、モデルの状態遷移関数ｆ（ｘ）及び状態関数ｇ（ｘ）のガウス過程によるモデル化が行われ、ガウス過程の表記を用いて、一般に、以下の式（３）、（４）のように表される。
　ｇ（ｘ）～ｇｐ（０，Ｋｇ）　　　（３）
　ｆ（ｘ）～ｇｐ（０，Ｋｆ）　　　（４）
式（３）は、状態の非線形な時間発展を表し、式（４）は、状態関数から観測関数への変換を表す。また、Ｋｇ、Ｋｆは、グラム行列を表す。
　式（１）の状態関数ｇｑ（ｘ）は、状態情報ｘｑから生成されるガウス過程のモデルｇｐ（０，ｋ（ｘｑ，ｘｑ´））を表す。
　ｋｑは、正定値カーネルを表す。使用する正定値カーネルとしては、ＲＢＦカーネル（Ｒａｄｉａｌ　ｂａｓｉｓ　ｆｕｎｃｔｉｏｎ　ｋｅｒｎｅｌ）等のデータに適したものを選択する。例えば、ＲＢＦカーネルは、以下の式（５）で与えられる。
Figure JPOXMLDOC01-appb-M000002
　式（５）において、
Figure JPOXMLDOC01-appb-M000003
 は、ＲＢＦカーネルの因果関係パラメータである。
　従来手法のＬＭＣの観測関数ｆｄ（Ｘ（ｔ））は、式（６）で表される。
Figure JPOXMLDOC01-appb-M000004
　これに対し、実施の形態１で提案するＬＭＣの観測関数ｆｄ（Ｘ（ｔ））は、以下の式（７）で表される。
Figure JPOXMLDOC01-appb-M000005
　従来手法の式（６）において、ａｄ，ｑは、従来手法のＬＭＣの線形相関を表し、状態変数ｑから観測変数ｄへの相関係数を表す。
　実施の形態１の式（７）において、
Figure JPOXMLDOC01-appb-M000006
 は、実施の形態１で新たに導入したｔｉｍｅ－ｓｈｉｆｔ　ｏｐｅｒａｔｏｒであり、Ｌｄ，ｑは、状態変数ｑから観測変数ｄへのラグ係数を表す。
　ＦＬは、時間軸に沿って状態関数の状態情報を未来（あるいは、過去）にずらすｔｉｍｅ－ｓｈｉｆｔ　ｏｐｅｒａｔｏｒ（ラグオペレータ）を表し、以下の式（８）のように定義される。
Figure JPOXMLDOC01-appb-M000007
　図７は、学習部２０の計算部２２のｔｉｍｅ－ｓｈｉｆｔ　ｏｐｅｒａｔｏｒによって行われる処理の例を示す図である。図７は、ｔｉｍｅ－ｓｈｉｆｔ　ｏｐｅｒａｔｏｒ　ＦＬの動作を示す。図７の例は、情報を時間軸ｔに沿って未来にずらす場合（Ｌ＞０）を示している。情報を時間軸ｔに沿って過去にずらす場合（Ｌ＜０）は、ずらす方向は図７の場合の方向の逆方向である。
　図７の例のように、状態情報のうち、時刻情報についてはｔｉｍｅ－ｓｈｉｆｔ　ｏｐｅｒａｔｏｒを適用しないことも可能であり、適用しない場合にはＬ＝０を適用する。
　図７の例では、計算部２２は、例えば、３次元時系列データである状態情報＃２にＬ＝１を適用しており、この場合、３次元時系列データのいずれも時間軸ｔに沿ってデータをシフトさせる。計算部２２は、シフト後、未来にはみ出した部分（すなわち、ｔ＝Ｔ＋１の部分）を削除し、空いた部分（ｔ＝０の部分）は、直前値（すなわち、ｔ＝１の部分の値）で穴埋めする。
　図７の例では、計算部２２は、例えば、２次元時系列データである状態情報＃３にＬ＝２を適用しており、この場合、２次元時系列データのいずれも時間軸ｔに沿ってデータをシフトさせる。計算部２２は、シフト後、未来にはみ出した部分（すなわち、ｔ＝Ｔ＋１、ｔ＝Ｔ＋２の部分）を削除し、空いた部分（すなわち、ｔ＝０、ｔ＝１の部分）は直前値（ｔ＝２）で穴埋めする。
　図７の例では、計算部２２は、例えば、３次元時系列データである状態情報＃ＱにＬ＝０を適用しており、この場合、３次元時系列データのシフトはなく、３次元時系列データは変更されない。
〈提案モデルの分散共分散行列の計算の仕方〉
　既存手法により、式（７）の
Figure JPOXMLDOC01-appb-M000008
 の観測データがガウス過程に従うと、ｆは、以下の式（９）の多次元ガウス分布で表される。
Figure JPOXMLDOC01-appb-M000009
　式（９）において、
Figure JPOXMLDOC01-appb-M000010
 である。
　Ｋ（Ｘ，Ｘ´）は、分散共分散行列（又は、グラム行列）と呼ばれ、状態変数間（すなわち、Ｘ⇔Ｘ´）の類似度を表す行列である。
　式（９）において、μ（Ｘ）は、状態変数の平均行列を表す。正規化部１８で状態変数を正規化しているので、μ（Ｘ）＝０である。
　分散共分散行列Ｋ（Ｘ，Ｘ´）の各成分（Ｋ（Ｘ，Ｘ´））ｄ，ｄは、式（７）のｆｄ（Ｘ）を用いて以下の式（１０）のように計算できる。
Figure JPOXMLDOC01-appb-M000011
　なお、式（１０）における２行目から３行目への式変換では、異なる状態間（ｑ≠ｑ´）は、独立であるため、
Figure JPOXMLDOC01-appb-M000012
 となる。
　分散共分散行列Ｋ（Ｘ，Ｘ´）は、以下の式（１１）及び式（１２）のように表される。
Figure JPOXMLDOC01-appb-M000013
　ここで、Ｂｑは、ｃｏｒｅｇｉｏｎａｌｉｚａｔｉｏｎ　ｍａｔｒｉｘと呼ばれ、状態関数から観測関数への線形変換を表し、以下のように表記される。
Figure JPOXMLDOC01-appb-M000014
　分散共分散行列Ｋ（Ｘ，Ｘ´）の計算は、状態変数の時系列データの集合Ｘ（０：ｔ）と、以下に示される因果関係パラメータθを用いて行うことができる。
Figure JPOXMLDOC01-appb-M000015
《１－１－５》最適化部２４
　図８は、因果関係パラメータθの初期値と個数の例を表形式で示す図である。これらの因果関係パラメータは、学習部２０の最適化部２４で最適化される。
　最適化部２４は、計算部２２によって計算された分散共分散行列Ｋ（Ｘ，Ｘ´）を受け取り、周辺尤度の計算と周辺尤度が最小化するように因果関係パラメータθを最適化する処理を行い、最適化した因果関係パラメータθを因果関係パラメータＤＢに保存する。
　周辺尤度は、以下のように、計算することで得られる。観測情報が観測される確率は、以下のように式（１３）で得ることができる。
Figure JPOXMLDOC01-appb-M000016
　式（１３）の両辺の対数をとると、以下の周辺尤度を計算することができる。ただし、Ｋθ（Ｘ，Ｘ´）は、因果関係パラメータを用いて計算した分散共分散行列Ｋ（Ｘ，Ｘ´）を表す。
Figure JPOXMLDOC01-appb-M000017
　ただし、Ｎは、特徴量ベクトルＸの長さ、Ｄは、ｙの出力次元数である。
　因果関係パラメータθの最適化には、周辺尤度ｌｏｇｐ（ｙ｜Ｘ，θ）が最大になるようにすればよい。一般的な最適化問題と形を合わせるため、式（１４）の両辺にマイナスをかけた以下の式（１５）で最適化関数Ｅを最小化するようにする。
Figure JPOXMLDOC01-appb-M000018
　最適化部２４は、式（１５）の最適化関数Ｅが最小化するようにθを更新する。最適化の際に、因果関係パラメータθを更新すると、Ｋθ（Ｘ，Ｘ´）の更新も必要になるため計算部２２でＫθ（Ｘ，Ｘ´）を計算する。
　最適化部２４は、最適化に際し、既知技術である確率的勾配降下法などを用いることができる。例えば、Ｌｄ，ｑは、グリッドサーチ用いて最適化することができ、残りの因果関係パラメータは、確率的勾配降下法を用いて最適化することができる。
　以上の動作を実行することにより、最適化された多変量因果関係の因果関係パラメータθを求めることができる。従来は、状態変数（Ｑ＝１）、観測変数（Ｄ＝１）にした場合の２変量の時系列データから「状態情報→観測情報」の因果関係を推定していた。これに対し、実施の形態１では、状態変数及び観測変数の各々を１次元以上とすることが可能であるため、任意数の変量の時系列データ間の因果関係の推定が可能であり、多変量（Ｑ＋Ｄ）の時系列データの因果関係を推定することができる。つまり、実施の形態１では、複数の状態変数と複数の観測変数との間の因果関係がわかるようになり、２変量の時系列データの因果関係の推定だけでなく、３変量以上の時系列データの因果関係の推定を行うことができる。
《１－２》動作
　図９は、実施の形態１に係る因果関係推定装置１００の動作（すなわち、学習方法）を示すフローチャートである。まず、ステップＳ１０１で、入力部１４は、時系列データＤＢ１２から状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）を取得する。
　次に、ステップＳ１０２で、前処理部１６は、状態変数の時系列データの集合Ｘ（０：ｔ）において次元変更が必要な状態変数ｘｑ（０：ｔ）の次元変更（すなわち、次元圧縮又は次元拡張）を行う。
　次に、ステップＳ１０３で、前処理部１６は、状態変数の時系列データの集合Ｘ（０：ｔ）と観測変数の時系列データの集合Ｙ（０：ｔ）の正規化を行う。
　次に、ステップＳ１０４で、計算部２２は、以下の因果関係パラメータθを初期値に設定する。
Figure JPOXMLDOC01-appb-M000019
　次に、ステップＳ１０５で、計算部２２は、状態変数の時系列データの集合Ｘ（０：ｔ）と因果関係パラメータθを用いて、式（１１）の分散共分散行列Ｋθ（Ｘ，Ｘ´）を計算する。
　次に、ステップＳ１０６で、最適化部２４は、式（１５）の最適化関数Ｅを計算する。
　次に、ステップＳ１０７で、最適化部２４は、最適化関数Ｅが最小になるように因果関係パラメータθを最適化（すなわち、更新）する。
　分散共分散行列Ｋθ（Ｘ，Ｘ´）の更新の際には、ステップＳ１０５で、更新された因果関係パラメータθを用いて、式（１１）のＫ（Ｘ，Ｘ´）が計算される。
　最適化関数Ｅが収束するまでステップＳ１０５～Ｓ１０７が繰り返される。最適化関数Ｅが収束したと判定されたときに、因果関係パラメータである相関係数ａｄ，ｑ及びラグ係数Ｌｄ，ｑが因果関係パラメータＤＢ２６に保存される。
《１－３》効果
　実施の形態１によれば、複数の状態変数と複数の観測変数との間の因果関係がわかるようになり、任意数の変量の時系列データの因果関係を推定することができる。
《２》実施の形態２
《２－１》構成
《２－１－１》因果関係推定装置２００
　因果関係推定装置２００は、因果関係パラメータＤＢから、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の階層ごとの因果関係を示す因果関係パラメータθを読み出し、因果関係パラメータに基づいて、複数の状態変数の時系列データと複数の観測変数の時系列データとを階層ごとに配置換えし、複数の状態変数の時系列データの集合Ｘ（０：ｔ）と複数の観測変数の時系列データの集合Ｙ（０：ｔ）とを含む検証用データを取得し、配置換えされた複数の状態変数の時系列データと複数の観測変数の時系列データに対し、グレンジャー因果の検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う。
　具体的に言えば、実施の形態２に係る因果関係推定装置２００は、因果グラフを構築する因果グラフ構築機能を備えた因果グラフ構築装置である。因果グラフとは、実施の形態１で得られた因果関係の情報に基づいて、「原因→結果」の関係でデータ項目一覧をグラフ構造にしたものである。構築された因果グラフの要素である状態関数及び観測関数は、因果関係パラメータＤＢに記憶された相関係数ａｄ，ｑ及びラグ係数Ｌｄ，ｑに基づいて、配置換えされる。構築された因果グラフは、グレンジャー因果と疑似相関とを用いて検証される。検証済みの因果グラフは、出力される。
　実施の形態２に係る因果関係推定装置２００は、因果関係パラメータである相関係数ａｄ，ｑ及びラグ係数Ｌｄ，ｑを用いて因果グラフを構築する。因果関係パラメータである相関係数ａｄ，ｑ及びラグ係数Ｌｄ，ｑは、例えば、実施の形態１に係る因果関係推定装置１００で作成された学習済みモデルの因果関係パラメータである。因果関係推定装置２００は、例えば、コンピュータである。因果関係推定装置２００は、実施の形態１に係る因果関係推定装置１００を構成するコンピュータと同じコンピュータであってもよく、又は、異なるコンピュータであってもよい。
　図１０は、実施の形態２に係る因果関係推定装置２００の構成を示すブロック図である。因果関係推定装置２００は、実施の形態２に係る因果関係推定（すなわち、因果グラフ構築方法）を実施することができる装置である。図１０に示されるように、因果関係推定装置２００は、構築部３０と、データ取得部４０と、出力部９０とを有している。構築部３０は、因果グラフ構築部３２と、因果グラフ検証部３４とを有している。構築部３０は、因果グラフ構築部３２に因果関係パラメータを提供する因果関係パラメータＤＢ８０を有してもよい。データ取得部４０は、入力部４４と、前処理部４６とを有している。データ取得部４０は、入力部４４に時系列データを提供する時系列データＤＢ４２を有してもよい。出力部９０は、構築部３０によって構築された因果グラフを出力する。時系列データＤＢ４２及び因果関係パラメータＤＢ８０は、因果関係推定装置２００の一部としての記憶部（後述の図１１の記憶部２０３）に格納されてもよいが、因果関係推定装置１００とは別の外部装置（例えば、因果関係推定装置２００と通信可能なネットワーク上のサーバ）の記憶部に格納されてもよい。
　図１１は、因果関係推定装置２００のハードウェア構成の例を示す図である。因果関係推定装置２００は、プロセッサ２０１と、メモリ２０２と、不揮発性の記憶装置である記憶部２０３とを有している。因果関係推定装置２００は、外部の装置とのインタフェース、ネットワークを介して他の装置と通信を行う通信部、などを備えてもよい。プロセッサ２０１は、ＣＰＵなどである。メモリ２０２は、例えば、ＲＡＭなどの、揮発性の半導体メモリである。記憶部２０３は、ＨＤＤ又はＳＳＤなどの記憶装置である。記憶部２０３は、情報（例えば、各種のデータベース）及びプログラムを記憶する。
　因果関係推定装置２００の各機能は、処理回路により実現される。処理回路は、専用のハードウェアであってもよいし、メモリ２０２に格納されるプログラムを実行するプロセッサ２０１であってもよい。
　処理回路がプロセッサ２０１である場合、実施の形態２に係る推定プログラム（すなわち、因果グラフ構築プログラム）は、ソフトウェア、ファームウェア、又はソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェア及びファームウェアは、プログラムとして記述され、メモリ２０２に格納される。プロセッサ２０１は、メモリ２０２に記憶されたプログラムを読み出して実行することにより、図１０に示される各部の機能を実現することができる。プログラムは、ネットワークを介してのダウンロードにより、又は、光ディスクなどのような情報を記録する記録媒体から、因果関係推定装置２００にインストールされる。
《２－１－２》データ取得部４０
　データ取得部４０は、実施の形態１のデータ取得部１０と同様の機能を有している。ただし、入力部４４は、後述するグレンジャー因果と疑似相関を用いた因果グラフの検証用のデータを取得する。
《２－１－３》因果グラフ構築部３２
　構築部３０の因果グラフ構築部３２は、因果関係パラメータＤＢ８０に保存されている因果関係パラメータ
Figure JPOXMLDOC01-appb-M000020
 を取得し、因果関係パラメータを用いて因果グラフの構築を行い、構築された因果グラフを因果グラフ検証部３４に出力する。
　図１２は、構築部３０の因果グラフ構築部３２による因果グラフの構築処理を示す図である。図１２において、ｙ（観測情報）よりもｘ（状態情報）がラグ（遅れ）があって、ｘ（状態情報）とｙ（観測情報）との間の相関係数が高い場合には、「ｘ´を原因としてｙ´が変動する」という因果関係（すなわち、ｘ´→ｙ´）がある可能性が高いと考えられる。
　図１２おいて、因果関係パラメータは、以下のとおりである。
Figure JPOXMLDOC01-appb-M000021
　これらの因果関係パラメータは、階層ごとに「状態情報→観測情報」の順番に並べられる。なお、ｈは、階層番号を示す正の整数である。
　図１３（Ａ）及び（Ｂ）は、構築部３０の因果グラフ構築部３２によって行われる状態情報及び観測情報の配置替え処理の例を示す図である。図１３（Ａ）及び（Ｂ）に示されるように、状態情報と観測情報との間は、状態情報から観測情報に向かう矢印で結ばれる。
　図１３（Ａ）においてより、状態関数よりも観測情報のほうがラグ（進み）がある場合には、矢印の向きを逆にする。図１３（Ｂ）より、相関係数が低い場合には、交絡因子又は中間因子である可能性を考え、因果グラフの上位などに持っていくような配置換えを行う。
　図１３（Ａ）及び（Ｂ）で配置換えを行った際には、配置換えを行った後の状態変数と観測変数を用いて実施の形態１に係る因果関係推定装置１００で因果関係パラメータを再度求める。
《２－１－４》因果グラフ検証部３４
　図１４（Ａ）及び（Ｂ）は、構築部３０の因果グラフ検証部３４によって行われる検証処理の例を示す図である。グレンジャー因果は、例えば、図１４（Ａ）に示されるように、状態変数ｘ２を含めて状態変数ｘ１，ｘ２，ｘ３によって観測変数ｙ１，ｙ２を予測した場合（ケース１）に比べて、図１４（Ｂ）に示されるように、状態変数ｘ２を削除して状態変数ｘ１，ｘ３によって観測変数ｙ１，ｙ２を予測した場合（ケース２）に予測精度が下がると、状態変数ｘ２は観測変数ｙ１，ｙ２に対してグレンジャー因果があると考える。また、状態変数ｘ１，ｘ２，ｘ３によって観測変数ｙ１，ｙ２を予測した場合（ケース１）に比べて、状態変数ｘ１，ｘ３によって観測変数ｙ１，ｙ２を予測した場合（ケース２）に予測精度が上がると、状態変数ｘ２は観測変数ｙ１，ｙ２に対してグレンジャー因果がないと考える。グレンジャー因果がない場合には、状態変数ｘ２の配置換え又は状態変数ｘ２の因果グラフからの削除を行う。
　予測には、因果関係推定装置１００で因果関係の推定に使用した状態変数の時系列データの集合Ｘ（０：ｔ）及び観測変数の時系列データの集合Ｙ（０：ｔ）から、未来の状態変数の時系列データの集合Ｘ（ｔ＋１：ｔ＋Δｔ）及び未来の観測変数の時系列データの集合Ｙ（ｔ＋１：ｔ＋Δｔ）を、予測の検証用データ（すなわち、テストデータ）として用いてもよい。
　また、訓練データで因果関係パラメータ
Figure JPOXMLDOC01-appb-M000022
 を求めて、既知技術のガウス過程回帰で、以下の予測値
Figure JPOXMLDOC01-appb-M000023
 を求めてもよい。
　予測誤差は、例えば、以下のＲＭＳＥ（Ｒｏｏｔ　Ｍｅａｎ　Ｓｑｕａｒｅｄ　Ｅｒｒｏｒ）で評価することができる。
Figure JPOXMLDOC01-appb-M000024
　図１５（Ａ）及び（Ｂ）は、構築部３０の因果グラフ検証部３４によって行われる検証処理の他の例を示す図である。図１５（Ａ）及び（Ｂ）より、非定常状態（ある因子が外的要因などで急激な変化が起こる（インパルス応答））で、どのくらい因子に影響するか（又は、どれくらい変化がなく一定値をとるか）、を見積もることで検証することも可能である。
　例えば、図１５（Ａ）に示される例では、状態変数ｘ２を介して疑似相関がある観測変数ｙ１とｙ２を用いる方法があると考えられる。図１５（Ｂ）に示される例では、状態変数ｘ２が観測変数ｙ１とｙ２の両方に影響を及ぼすが、状態変数ｘ１とｘ３は、観測変数ｙ１とｙ２のそれぞれ以外に影響を及ぼさない場合を考える。ｘ２が固定値をとるなど変化がない場合には、ｘ２が変化しないため観測変数ｙ１とｙ２に影響を及ばさなくなるため、観測変数ｙ１とｙ２の間には、相関がなくなると考えられる。
　疑似相関を用いた因果関係の検証において、ｘ２の影響がなくなっても観測変数ｙ１とｙ２の相関が残っている、あるいは、ｘ１と観測変数ｙ１、ｘ２とｙ２が高い相関を示されるようになる場合には、状態変数ｘ２の配置換え又は因果グラフからの削除を行う。
　因果グラフを構築することで、直接測定できない場所のセンサ値の予測又は異常原因の探索等を可能となる。適用事例として交通量又は気象情報からの路面損傷の予測又は、付近の雨量と水位データからの河川の水位予測、天候又は周辺地域の経済状況データを活用した電力需要予測などが考えられる。
《２－２》動作
　図１６は、因果関係推定装置２００の動作（すなわち、推論動作）を示すフローチャートである。まず、ステップＳ２０１で、構築部３０は、因果関係パラメータＤＢから各階層の
Figure JPOXMLDOC01-appb-M000025
 を取得する。
　次に、ステップＳ２０２で、因果グラフ構築部３２は、状態情報及び観測情報を、階層ごとに「状態情報→観測情報」の順に並べる。
　次に、ステップＳ２０３で、因果グラフ構築部３２は、因果方向（すなわち、矢印の向き）を修正する。
　次に、ステップＳ２０４で、因果グラフ構築部３２は、因果グラフの状態情報及び観測情報の配置換えを行う。
　次に、ステップＳ２０５で、入力部４４は、時系列データＤＢ４２からグレンジャー因果の検証用データを取得する。
　次に、ステップＳ２０６で、前処理部４６は、検証用データの次元変更を行う。
　次に、ステップＳ２０７で、前処理部４６は、検証用データの正規化を行う。
　次に、ステップＳ２０８で、因果グラフ検証部３４は、グレンジャー因果を用いて因果グラフを検証する。
　次に、ステップＳ２０９で、入力部４４は、時系列データＤＢ４２から疑似相関による検証用のデータを取得する。
　次に、ステップＳ２１０で、前処理部４６は、検証用データの次元変更を行う。
　次に、ステップＳ２１１で、前処理部４６は、検証用データの正規化を行う。
　次に、ステップＳ２１２で、因果グラフ検証部３４は、疑似相関を用いた因果関係の検証をする。
　以上に示した処理に基づいて、因果グラフ検証部３４が検証した因果グラフは出力部９０に出力される。
《２－３》効果
　実施の形態２によれば、因果グラフを構築することで、直接測定できない場所のセンサ値の予測又は異常原因の探索等が可能となる。
　また、実施の形態２は、例えば、交通量又は気象情報からの路面損傷の予測、付近の雨量と水位データからの河川の水位予測、天候又は周辺地域の経済状況データを活用した電力需要予測、などに適用可能である。
《２－４》適用例
［観測値の予測］
　図１７（Ａ）及び（Ｂ）は、因果関係推定装置２００によって時系列データ間の関係性から観測されていない時系列データを予測する処理を行う学習部と推論部とを示す。この場合、因果関係推定装置１００は、図１７（Ａ）に示される過去のＡ１駅、Ａ２駅、Ａ３駅の混雑情報をマルチタスク学習することで、学習済みモデルを生成する。因果関係推定装置１００は、Ａ１駅、Ａ２駅、Ａ３駅のうちの、Ａ１駅とＡ２駅の混雑情報のみがわかり、Ａ３駅の混雑情報がわからない場合（図１７（Ｂ）の場合）に、マルチタスク学習で生成された学習済みモデルを用いて、お互いの混雑度の相関又はラグ情報からＡ３駅の混雑度を推論部で予測（推論）することが可能である。
［異常検知又は故障予知］
　図１８は、因果関係推定装置１００にＣｈａｎｇｅ　Ｐｏｉｎｔ　Ｋｅｒｎｅｌを導入することでセンサデータの変化点又は故障を学習した場合の動作を示す図である。センサデータの変化点又は故障の学習によって生成された学習済みモデルを用いることで、故障予知を行うことができる。
Figure JPOXMLDOC01-appb-M000026
［Ｃａｕｓａｌ　Ｉｍｐａｃｔへの適用］
　公知の時系列因果の推論フレームワーク（例えば、Ｃａｕｓａｌ　Ｉｍｐａｃｔ）において、標準の線形状態空間モデルを、実施の形態の方法（ＧＰＤＭ＋ＭＴＧＰモデル）で置き換えることも可能である。この場合、非線形又は非定常の相関も考慮できるようになり、予測精度を向上させることができる。
［計算の高速化］
　上述したＭＴＧＰモデルの計算は、計算量Ｏ（Ｄ３Ｎ３）、メモリＯ（Ｄ２Ｎ２）と計算負荷又はメモリコストが高い。このため、例えば、乱択化フーリエ特徴（Ｒａｎｄｏｍ　Ｆｏｕｒｉｅｒ　ｆｅａｔｕｒｅｓ）又は変分化フーリエ特徴（Ｖａｒｉａｔｉｏｎａｌ　Ｆｏｕｒｉｅｒ　ｆｅａｔｕｒｅｓ）を用いることで、マルチタスク学習を高速化し、メモリコストを削減することが可能である。
　１００　因果関係推定装置、　２００　因果関係推定装置（因果グラフ構築装置）、　１０　データ取得部、　１２、４２　時系列データＤＢ、　１４　入力部、　１６　前処理部、　１７　次元変更部、　１８　正規化部、　２０　学習部、　２２　計算部、　２４　最適化部、　２６、８０　因果関係パラメータＤＢ、　３０　構築部、　３２　因果グラフ構築部、　３４　因果グラフ検証部、　４０　データ取得部、　４４　入力部、　４６　前処理部、　９０　出力部。

Claims (12)

  1. 　複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するデータ取得部と、
    　前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現する計算部と、
    　前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新する最適化部と、
    　を有することを特徴とする因果関係推定装置。
  2. 　前記因果関係パラメータは、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の相関を示す相関係数と、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間のラグを示すラグ係数とを含み、
    　前記計算部は、前記相関を、前記マルチタスクガウス過程モデルのＬＭＣカーネルの線形相関で表現する
    　ことを特徴とする請求項１に記載の因果関係推定装置。
  3. 　前記データ取得部は、前記複数の状態変数の時系列データの次元変更を行い、次元変更された前記状態変数の時系列データを前記計算部に提供する
    　ことを特徴とする請求項１又は２に記載の因果関係推定装置。
  4. 　前記複数の状態変数の時系列データは、時刻情報の時系列データを含む
    　ことを特徴とする請求項１から３のいずれか１項に記載の因果関係推定装置。
  5. 　前記複数の状態変数の時系列データは、角度情報の時系列データを含む
    　ことを特徴とする請求項１から３のいずれか１項に記載の因果関係推定装置。
  6. 　前記最適化部は、更新された前記因果関係パラメータを因果関係パラメータデータベースに保存する
    　ことを特徴とする請求項１から５のいずれか１項に記載の因果関係推定装置。
  7. 　前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の階層ごとの因果関係を示す前記因果関係パラメータに基づいて、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとを前記階層ごとに配置換えする因果グラフ構築部と、
    　複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む検証用データを取得する他のデータ取得部と、
    　配置換えされた前記複数の状態変数の時系列データと前記複数の観測変数の時系列データに対し、グレンジャー因果の前記検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う因果グラフ検証部と、
    　を有することを特徴とする請求項１から６のいずれか１項に記載の因果関係推定装置。
  8. 　因果関係パラメータデータベースから、複数の状態変数の時系列データと複数の観測変数の時系列データとの間の階層ごとの因果関係を示す因果関係パラメータを読み出し、前記因果関係パラメータに基づいて、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとを前記階層ごとに配置換えする因果グラフ構築部と、
    　複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む検証用データを取得するデータ取得部と、
    　配置換えされた前記複数の状態変数の時系列データと前記複数の観測変数の時系列データに対し、グレンジャー因果の前記検証用データを用いた検証と疑似相関の前記検証用データを行いた検証との一方又は両方を行う因果グラフ検証部と、
    　を有することを特徴とする因果関係推定装置。
  9. 　前記因果関係パラメータは、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の相関を示す相関係数と、前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間のラグを示すラグ係数とを含む
    　ことを特徴とする請求項８に記載の因果関係推定装置。
  10. 　予測対象に関する複数の状態変数の時系列データと複数の観測変数の時系列データとを取得し、前記予測対象に関する前記複数の状態変数の時系列データと前記複数の観測変数の時系列データから未観測の次元の観測情報を予測するための、前記因果関係パラメータに基づく学習済みモデルを用いて、前記未観測の次元の観測情報を予測する推論部をさらに有する
    　ことを特徴とする請求項１から９のいずれか１項に記載の因果関係推定装置。
  11. 　因果関係推定装置によって実施される因果関係推定方法であって、
    　複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するステップと、
    　前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現するステップと、
    　前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新するステップと、
    　を有することを特徴とする因果関係推定方法。
  12. 　複数の状態変数の時系列データの集合と複数の観測変数の時系列データの集合とを含む学習用データを取得するステップと、
    　前記複数の状態変数の時系列データと前記複数の観測変数の時系列データとの間の因果関係を示す因果関係パラメータを計算し、前記学習用データと前記因果関係パラメータとからガウス過程の分散共分散行列を計算し、前記因果関係パラメータをマルチタスクガウス過程モデルで表現するステップと、
    　前記分散共分散行列に基づいて最適化関数を計算し、前記最適化関数に基づいて前記因果関係パラメータを更新するステップと、
    　をコンピュータに実行させることを特徴とする因果関係推定プログラム。
PCT/JP2022/020680 2022-05-18 2022-05-18 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム WO2023223461A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2022/020680 WO2023223461A1 (ja) 2022-05-18 2022-05-18 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム
JP2024513679A JP7483180B2 (ja) 2022-05-18 2022-05-18 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2022/020680 WO2023223461A1 (ja) 2022-05-18 2022-05-18 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム

Publications (1)

Publication Number Publication Date
WO2023223461A1 true WO2023223461A1 (ja) 2023-11-23

Family

ID=88834904

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2022/020680 WO2023223461A1 (ja) 2022-05-18 2022-05-18 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム

Country Status (2)

Country Link
JP (1) JP7483180B2 (ja)
WO (1) WO2023223461A1 (ja)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2019057198A (ja) * 2017-09-22 2019-04-11 株式会社神戸製鋼所 油圧システムのパラメータ推定方法
US20200241171A1 (en) * 2016-05-12 2020-07-30 The Climate Corporation Statistical blending of weather data sets
JP2022054018A (ja) * 2020-09-25 2022-04-06 株式会社Ｇｓユアサ 推定装置、推定方法、及びコンピュータプログラム

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200241171A1 (en) * 2016-05-12 2020-07-30 The Climate Corporation Statistical blending of weather data sets
JP2019057198A (ja) * 2017-09-22 2019-04-11 株式会社神戸製鋼所 油圧システムのパラメータ推定方法
JP2022054018A (ja) * 2020-09-25 2022-04-06 株式会社Ｇｓユアサ 推定装置、推定方法、及びコンピュータプログラム

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
FUKUMIZU KENJI, BACH FRANCIS R, JORDAN MICHAEL I: "Dimensionality reduction for supervised learning with reproducing kernel hilbert spaces", THE JOURNAL OF MACHINE LEARNING RESEARCH, vol. 5, 1 January 2004 (2004-01-01), pages 73 - 99, XP093110202 *
YOSHIYUKI KOBORI, MASATOSHI SATO, MAMORU TANAKA: "Function resolution and synthesis of sound by Covariance Structure Analysis -- Modeling of a synthesis sound of musical note --", IEICE TECHNICAL REPORT, NLP, IEICE, JP, vol. 111, no. 106, NLP2011-43, 23 June 2011 (2011-06-23), JP, pages 107 - 112, XP009550667 *

Also Published As

Publication number Publication date
JPWO2023223461A1 (ja) 2023-11-23
JP7483180B2 (ja) 2024-05-14

Similar Documents

Publication Publication Date Title
US6928398B1 (en) System and method for building a time series model
KR102044205B1 (ko) 빅데이터와 기계학습을 이용한 타겟 정보 예측 시스템 및 예측 방법
CN110008301B (zh) 基于机器学习的区域性地质灾害易发性预测方法及装置
Lall et al. A nearest neighbor bootstrap for resampling hydrologic time series
JP2022548654A (ja) 機械学習モデルにおいて動的外れ値偏り低減を実装するように構成されるコンピュータベースシステム、コンピュータコンポーネント及びコンピュータオブジェクト
US20060217939A1 (en) Time series analysis system, time series analysis method, and time series analysis program
EP1760657A2 (en) Methods and systems for assessing loss severity for commercial loans
US20160004800A1 (en) Methods and systems for reservoir history matching for improved estimation of reservoir performance
CN107016571A (zh) 数据预测方法及其系统
JP6845126B2 (ja) 故障確率算出装置、故障確率算出方法及びプログラム
KR102342321B1 (ko) 딥러닝 기반의 주식 스크리닝과 이를 이용한 포트폴리오 자동화 및 고도화 방법 및 장치
WO2016149906A1 (en) Analyzing equipment degradation for maintaining equipment
CN114580263A (zh) 基于知识图谱的信息系统故障预测方法及相关设备
KR102330423B1 (ko) 이미지 인식 딥러닝 알고리즘을 이용한 온라인 부도 예측 시스템
AU2019371339A1 (en) Finite rank deep kernel learning for robust time series forecasting and regression
CN112598248A (zh) 负荷预测方法、装置、计算机设备和存储介质
US20200265307A1 (en) Apparatus and method with multi-task neural network
JP4443619B2 (ja) ポートフォリオの信用リスクの計算方法および装置
CN111444649A (zh) 基于强度折减法的边坡系统可靠度分析方法
EP4009239A1 (en) Method and apparatus with neural architecture search based on hardware performance
WO2023223461A1 (ja) 因果関係推定装置、因果関係推定方法、及び因果関係推定プログラム
US20210042820A1 (en) Extending finite rank deep kernel learning to forecasting over long time horizons
WO2021168346A1 (en) Systems and methods for anonymizing sensitive data and simulating accelerated schedule parameters using the anonymized data
Leong et al. Control variates for efficient long-term extreme analysis of mooring lines
JP6615892B2 (ja) 物理システムの経時変化プロファイリングエンジン

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22942666

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2024513679

Country of ref document: JP

Kind code of ref document: A