WO2021157669A1 - 回帰分析装置、回帰分析方法及びプログラム - Google Patents

回帰分析装置、回帰分析方法及びプログラム Download PDF

Info

Publication number
WO2021157669A1
WO2021157669A1 PCT/JP2021/004178 JP2021004178W WO2021157669A1 WO 2021157669 A1 WO2021157669 A1 WO 2021157669A1 JP 2021004178 W JP2021004178 W JP 2021004178W WO 2021157669 A1 WO2021157669 A1 WO 2021157669A1
Authority
WO
WIPO (PCT)
Prior art keywords
coefficient
regression
positive
training data
constraint
Prior art date
Application number
PCT/JP2021/004178
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 CN202180012527.4A priority Critical patent/CN115053216A/zh
Priority to JP2021575868A priority patent/JPWO2021157669A1/ja
Priority to EP21751487.6A priority patent/EP4102420A4/en
Priority to US17/797,141 priority patent/US20230059056A1/en
Publication of WO2021157669A1 publication Critical patent/WO2021157669A1/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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Definitions

  • This disclosure relates to a regression analyzer, a regression analysis method and a program.
  • Non-Patent Document 1 a method of giving a constraint condition called an L1 norm has been proposed (for example, Non-Patent Document 1).
  • LASSO Least Absolute Shrinkage and Selection Operator
  • selection of an explanatory variable suitable for explaining the objective variable and determination of a coefficient are also performed.
  • the regression analyzer defines in advance the training data used as the objective variable and the explanatory variable of the regression model, and whether the explanatory variable should be changed to positive or negative in order to change the objective variable in the positive or negative direction.
  • the training data is stored so as to minimize the data acquisition unit that reads the training data and the constraint condition from the storage device that stores the constraint condition, and the cost function including the regularization term that increases the cost when the constraint condition is violated. It is provided with a coefficient update unit that repeatedly updates the coefficient of the explanatory variable in the regression model.
  • the coefficient that violates the constraint condition is not selected, and the explanatory variable should be changed to positive or negative in order to change the objective variable in the positive or negative direction.
  • the regularization term may increase the cost according to the sum of the absolute values of the coefficients in the interval where the coefficient is positive or negative according to the constraint condition.
  • a regression model using L1 regularization may be constructed on one side with positive or negative coefficients.
  • the regularization term increases the cost according to the sum of the absolute values of the coefficients in one of the intervals where the coefficient is positive or negative according to the constraint condition, and makes the cost infinite on the other side. May be good.
  • the coefficient update unit may set the coefficient to zero when the coefficient does not converge to a value satisfying the constraint condition. In this way, the explanatory variables that do not contribute to the objective variable can be deleted from the regression model under the above-mentioned constraints, and sparse modeling is realized.
  • the coefficient updating unit may update the coefficient by the proximity gradient method. In this way, it is possible to avoid passing through the non-differentiable point of the regularization term in the convergence calculation. Therefore, the time required for convergence can be shortened.
  • the contents described in the means for solving the problems can be combined as much as possible without departing from the problems and technical ideas of the present disclosure. Further, the content of the means for solving the problem can be provided as a device such as a computer or a system including a plurality of devices, a method executed by the computer, or a program executed by the computer. A recording medium for holding the program may be provided.
  • FIG. 1 is a diagram showing an example of training data used for creating a regression equation.
  • FIG. 2A is a schematic diagram for explaining the constraints imposed on the regression coefficients.
  • FIG. 2B is a schematic diagram for explaining the constraints imposed on the regression coefficients.
  • FIG. 3 is a diagram for explaining the update of the parameter w.
  • FIG. 4 is a diagram for explaining the update of the parameter ⁇ .
  • FIG. 5 is a block diagram showing an example of the configuration of the regression analyzer 1 that performs the regression analysis described above.
  • FIG. 6 is a processing flow diagram showing an example of the regression analysis process executed by the regression analyzer.
  • FIG. 7A is a diagram showing the relationship between the parameter ⁇ representing the strength of the constraint and the correlation coefficient r.
  • FIG. 7A is a diagram showing the relationship between the parameter ⁇ representing the strength of the constraint and the correlation coefficient r.
  • FIG. 7B is a diagram showing the relationship between the parameter ⁇ representing the strength of the constraint and the correlation coefficient r.
  • FIG. 8 is a diagram showing the relationship between the parameter ⁇ representing the strength of the constraint and the coefficient of determination E.
  • FIG. 9 is a diagram showing the relationship between the number of data T used for learning and the correlation coefficient r.
  • FIG. 10 is a diagram showing the relationship between the number of data T used for learning and the coefficient of determination E.
  • FIG. 11A is a schematic diagram for explaining the constraints imposed on the regression coefficients.
  • FIG. 11B is a schematic diagram for explaining the constraints imposed on the regression coefficients.
  • FIG. 12 is a diagram showing the relationship between the parameter ⁇ and the correlation coefficient r.
  • Figure 13 is a diagram showing the relationship between the coefficient of determination R 2 and parameter beta.
  • FIG. 14 is a diagram showing the relationship between the parameter ⁇ and RMSE.
  • the regression analyzer constructs a regression equation (regression model) representing the relationship between one or more explanatory variables (independent variables) and one objective variable (dependent variable). At this time, at least one of the explanatory variables is constrained to have a certain correspondence between the positive or negative direction of the fluctuation of the explanatory variable and the positive or negative direction of the fluctuation of the objective variable ( Create a regression equation by imposing a "sign constraint").
  • FIG. 1 is a diagram showing an example of observation values (training data) used for creating a regression equation.
  • the table of FIG. 1 includes columns of K types of inputs x (x 1 to x K ) and columns of outputs y.
  • the input x corresponds to the explanatory variable and the output y corresponds to the objective variable.
  • a regression equation is created using T records among a plurality of records representing data points t (t 1 to t T , ...) Which are individual samples of training data.
  • at least a part of the input x of the K type is associated with a positive or negative code (information representing the constraint condition according to the present embodiment, which is referred to as a "constraint code").
  • the constraint code associated with each input x determines in which direction the input x should be changed, positive or negative, in order to change the output y in the positive direction in the regression equation to be constructed. This is information for defining in advance.
  • the regression equation is represented by, for example, the following equation (1). Note that w k is a regression coefficient and w 0 is a constant term. Further, w k is determined according to a predetermined constraint code.
  • the cost function represented by the following equation (2) can be used to determine the regression coefficient and the constant term.
  • the regression equation is determined by selecting a coefficient w k that minimizes the cost function E (w).
  • ⁇ R is a regularization term (penalty term), and its coefficient ⁇ is a parameter representing the strength of the constraint. If the constraint sign of x k in the table of FIG. 1 is positive, a value of R + (w), if the constraint code is negative, R - takes the value of (w).
  • the regularization term ⁇ R according to the present embodiment imposes a sign constraint by L1 type regularization on one side of positive or negative. That is, the regularization term increases the cost according to the sum of the absolute values of the coefficients in either the positive or negative interval in which the coefficient w k corresponds to the constraint code.
  • FIGS. 2A and 2B are schematic diagrams for explaining the constraint imposed on one regression coefficient w.
  • the vertical axis represents R + (w) and the horizontal axis represents w.
  • the arrow schematically indicates that the regularization term is defined so that the value of R + (w) becomes larger as the value of ⁇ becomes larger in the interval where w is negative.
  • the predicted value ⁇ by the regression equation decreases as the input x k of the regression equation shown in the equation (1) increases. That is, if the constraint code associated with the input x k is negative, small regularization term when the value of the input x k is the value of the predicted value ⁇ decreases as increasing the value of the input x k is increased.
  • the cost function is defined so that the regularization term becomes larger when the value of the predicted value ⁇ also increases.
  • regression is imposed by imposing a constraint that has a certain correspondence between the positive or negative direction of the fluctuation of the explanatory variable and the positive or negative direction of the fluctuation of the objective variable. Perform an analysis.
  • the parameter w that minimizes E (w) may be updated by using the following equation (4), for example, by the gradient method.
  • FIG. 3 is a diagram for explaining the update of the parameter w. Based on the gradient of the variable w of the cost function E (w) in a certain step s, the variable w in the later step s + 1 is updated, and such a process is repeated until w converges.
  • a value corresponding to the constraint code may be calculated for each input x k , and the sum of the values may be used as a regularization term to perform regression by the steepest descent method, but the calculation becomes unstable. Therefore, for example, the proximity gradient method may be used. Also in the proximity gradient method, for example, w that minimizes the above-mentioned equation (2) is obtained. Assuming that the sum-of-squares error of equation (2) is f (x) and the regularization term is g (w), the update equation of w is expressed by the following equation (5).
  • ⁇ is a step width that determines the size of updating the coefficient w in one step (one iteration).
  • ⁇ f (W (t)) is a gradient. The update is repeated until the gradient approaches zero sufficiently, and when the gradient approaches zero sufficiently, it is judged that the update has converged and the update is terminated.
  • the update formula of w is represented by the following formula (6).
  • the constraint sign When the constraint sign is positive, it can be calculated by the following equation (7).
  • the constraint sign When the constraint sign is negative, it can be calculated by the following equation (8).
  • the coefficient w can be determined by the above processing. The coefficient w converges to a value that satisfies the code constraint and contributes to the objective variable, and if there is no such value, the coefficient w approaches zero. That is, when there is no value that satisfies the code constraint, as shown in FIGS. 2A and 2B, the penalty effect of regularization works and the value contrary to the code constraint is pulled back, and as a result, it converges to zero. Therefore, a part of the regression coefficient can be estimated to be zero as in the so-called LASSO.
  • FIG. 4 shows an example of a schematic code for searching for a suitable ⁇ .
  • ⁇ 0 is a predetermined initial value.
  • ⁇ is, for example, a positive value less than 1, and is updated to decrease ⁇ .
  • FIG. 5 is a block diagram showing an example of the configuration of the regression analyzer 1 that performs the regression analysis described above.
  • the regression analyzer 1 is a general computer, and includes a communication interface (I / F) 11, a storage device 12, an input / output device 13, and a processor 14.
  • the communication I / F 11 may be, for example, a network card or a communication module, and communicates with another computer based on a predetermined protocol.
  • the storage device 12 includes a main storage device such as a RAM (Random Access Memory) and a ROM (Read Only Memory), and an auxiliary storage device (secondary) such as an HDD (Hard-Disk Drive), an SSD (Solid State Drive), and a flash memory. It may be a storage device).
  • the main storage device temporarily stores the program read by the processor 14 and the information processed by the program.
  • the auxiliary storage device stores a program executed by the processor 14, information processed by the program, and the like.
  • the storage device 12 temporarily or permanently stores training data and information representing the constraint conditions.
  • the input / output device 13 is, for example, a user interface such as an input device such as a keyboard and a mouse, an output device such as a monitor, and an input / output device such as a touch panel.
  • the processor 14 is an arithmetic processing unit such as a CPU (Central Processing Unit), and executes each process according to the present embodiment by executing a program.
  • a functional block is shown in the processor 14. That is, the processor 14 functions as a data acquisition unit 141, a coefficient update unit 142, a convergence determination unit 143, a verification processing unit 144, and an operation processing unit 145 by executing a predetermined program.
  • the data acquisition unit 141 acquires training data and information representing the constraint conditions from the storage device 12.
  • the coefficient updating unit 142 updates the coefficient of the regression equation under the above-mentioned constraint conditions. Further, the convergence test unit 143 determines whether or not the updated coefficient values have converged. If it is determined that the coefficients have not converged, the coefficient updating unit 142 repeats updating the coefficients. When it is determined that the coefficients have converged, for example, the coefficient updating unit 142 stores the finally generated coefficient in the storage device 12. Further, the verification processing unit 144 evaluates the regression equation created based on the predetermined evaluation index. The operation processing unit 145 calculates the predicted value by using the created regression equation and, for example, the newly acquired observed value.
  • the operation processing unit 145 may calculate a predicted value when the condition is changed by using the created regression equation and an arbitrary value.
  • the arbitrary value may be a value input by the user via, for example, the communication I / F 11 or the input / output device 13. Since the regression equation created in this embodiment has a certain correspondence between the direction of fluctuation of the explanatory variable and the direction of fluctuation of the objective variable, for example, the input value is increased in order to bring the predicted value closer to a desired value. The user can easily estimate whether to reduce or decrease. Therefore, for example, the regression equation according to the present embodiment is effective when some control is performed based on the estimated value.
  • FIG. 6 is a processing flow diagram showing an example of the regression analysis process executed by the regression analyzer.
  • the data acquisition unit 141 of the regression analysis device 1 reads out from the training data, the information representing the constraint condition, and the storage device 12 (FIG. 6: S11).
  • the values of the input x and the output y as shown in FIG. 1 are read out as training data. It is assumed that the input x is treated as an explanatory variable and the output y is treated as an objective variable.
  • a positive or negative code registered in association with the input x is read out as information representing a constraint condition.
  • the regression analyzer 1 uses the code to be read as the constraint code described above. In this embodiment, a regression equation as shown in the equation (1) is used.
  • the coefficient updating unit 142 of the regression analyzer 1 updates the regression coefficient under the above-mentioned code constraint (FIG. 6: S12).
  • the coefficient updating unit 142 updates the coefficient w so as to minimize the cost function E (w) shown in the equation (2), for example, as shown by the arrow on the upper side in FIG.
  • the coefficient updating unit 142 can update the coefficient w based on the equations (6) to (8).
  • the regularization term of the cost function E (w) is defined so that the cost increases when the constraint condition acquired in S11 is not satisfied. That is, the regularization term is a cost function E (w) when the positive or negative direction of the variation of the explanatory variable and the positive or negative direction of the variation of the objective variable have a predetermined correspondence relationship. It reduces the value of. Further, the coefficient updating unit 143 sets the coefficient to zero when the coefficient does not converge to a value satisfying the constraint condition.
  • the convergence test unit 143 of the regression analyzer 1 determines whether the coefficient w has converged or the coefficient w has been set to zero (FIG. 6: S13). In this step, the convergence test unit 143 determines that the coefficient w has converged when the gradient of the coefficient w to be updated approaches zero sufficiently. Specifically, the convergence test unit 143 determines that the coefficient has converged when the value of the coefficient w does not change in the equation (7) or the equation (8).
  • the convergence determination unit 143 stores the regression equation in the storage device 12 (FIG. 6: S14). In this step, the convergence test unit 143 stores the updated coefficient w in the storage device 12.
  • the verification processing unit 144 of the regression analyzer 1 may verify the accuracy of the created regression equation (FIG. 6: S20). In this step, the verification processing unit 144 verifies the accuracy of the regression equation using test data, for example, by cross-validation. Further, the verification processing unit 144 can perform verification based on a predetermined evaluation index such as a correlation coefficient and a predetermined coefficient of determination. As will be described later, this step may be omitted.
  • the operation processing unit 145 of the regression analyzer 1 performs an operation process using the created regression equation (FIG. 6: S30).
  • the operation processing unit 145 calculates the predicted value of the output y for the new input x, for example, the record whose data number is t T + 1 shown in FIG. Note that this step may be performed by an apparatus (not shown) other than the regression analyzer 1 using the regression equation stored in S14.
  • Example> A regression equation was constructed using the sensing data obtained from the production plant, and the accuracy was evaluated.
  • the output values of different sensors were used as the inputs and outputs shown in FIG. Further, for the sensing data continuously output from the sensor, the latest data number T was set as the learning interval.
  • the constraint code was preset based on the knowledge about the production plant.
  • the correlation coefficient r used as an evaluation index is obtained by the following equation (9). That is, the numerator of the formula (9) is the covariance of the predicted value ⁇ and the measured value y of the training data.
  • the denominator of equation (9) is the product of the standard deviation of the predicted value ⁇ and the standard deviation of the measured value y of the training data.
  • the coefficient of determination E used as another evaluation index is obtained by the following equation (10).
  • the coefficient of determination E is a value representing the size of the distribution of predicted values with respect to the distribution of observed values.
  • FIGS. 7A and 7B are diagrams showing the relationship between the parameter ⁇ representing the strength of the constraint and the correlation coefficient r for the model constructed by the plurality of methods.
  • FIG. 8 is a diagram showing the relationship between the parameter ⁇ representing the strength of the constraint and the coefficient of determination E for the model constructed by the plurality of methods.
  • the horizontal axis represents the parameter ⁇ and the vertical axis represents the correlation coefficient r. 7A and 7B have different scales on the horizontal axis.
  • the horizontal axis represents ⁇ and the vertical axis represents the coefficient of determination E.
  • the solid line is the method disclosed in the embodiment, the broken line is a comparative example in which a part of the sign constraint of the embodiment is randomly selected and the positive and negative are reversed. Represents each result of.
  • the model was constructed with the number of data T as 40.
  • the constraint code is preset based on the knowledge about the production plant, but in general, it may include an inappropriate setting. It can be said that the comparative example simulates an erroneous code constraint.
  • the value of the correlation coefficient r was higher in the order of the method of the present disclosure, the comparative example, Lasso, and no constraint. Further, as shown in FIG. 8, the coefficient of determination E was close to 1 in the order of Lasso, the method and comparative example of the present disclosure, and no constraint.
  • the coefficient of determination E was close to 1 in the order of Lasso, the method and comparative example of the present disclosure, and no constraint.
  • ⁇ is a so-called hyperparameter and needs to be adjusted by cross-validation.
  • the accuracy can be improved by setting ⁇ sufficiently large. This has the effect that manual parameter adjustment can be eliminated.
  • the correlation coefficient r is lower than that of the method according to the embodiment. That is, the method according to the embodiment is particularly applicable when the data to be analyzed has a certain correspondence relationship between the fluctuation of the explanatory variable and the fluctuation of the objective variable and a code constraint matching the correspondence is given. It can be said that a good model can be created. Further, as can be seen from FIGS. 7A and 7B, even in the case of the comparative example in which the code constraint is randomly given as shown by the broken line, the correlation coefficient is higher than in the case of no regularity shown by the alternate long and short dash line.
  • FIG. 9 is a diagram showing the relationship between the number of data T used for learning and the correlation coefficient r for a model constructed by a plurality of methods.
  • FIG. 10 is a diagram showing the relationship between the number of data T used for learning and the coefficient of determination E for a model constructed by a plurality of methods.
  • the coefficient of determination E was close to 1 in the order of Lasso, the method and comparative example of the present disclosure, and no constraint.
  • the method of the present disclosure can be said to be effective when the training data is relatively small. In other words, it is also useful when sufficient data cannot be collected, or when the prediction model changes from moment to moment but there are changes in the state that cannot be observed from the data alone, so only the latest data can be used. be.
  • a regression equation that satisfies a constraint that has a certain correspondence between the positive or negative direction of the fluctuation of the explanatory variable and the positive or negative direction of the fluctuation of the objective variable is generated. can do. Therefore, the user can use the regression equation to know whether the value of the input x k should be changed to positive or negative in order to bring the predicted value ⁇ closer to the desired value. Further, as described with reference to FIGS. 7A and 7B, there is an advantage that it is not necessary to adjust the parameter ⁇ indicating the strength of the constraint. Further, as described with reference to FIGS. 9 and 10, the method of the present disclosure is particularly effective when the training data is relatively small.
  • w k is obtained as follows. Moreover, when this is solved again, it is obtained as follows.
  • w k can be expressed by the following equation (11) without considering the case of the lower stage.
  • the coefficient w k is set to zero as shown in the lower part of the equation (11).
  • the explanatory variable x k that cannot satisfy the constraint is not used in the created regression equation. Therefore, it is possible to generate a regression equation that satisfies the constraint that there is a certain correspondence between the positive or negative direction of the fluctuation of the explanatory variable and the positive or negative direction of the fluctuation of the objective variable. Further, the value of the parameter ⁇ can be set to a sufficiently large value, and it can be said that adjustment is unnecessary.
  • w k is calculated as follows, for example. That is, it is estimated to be biased so as to be smaller by ⁇ from the value that should originally converge. Such a bias acts to increase the squared error. On the other hand, according to the technique of the present disclosure, such a bias does not occur, so that it can be said that the accuracy of the regression equation is improved.
  • the Oracle property (Oracle property, Fan and Li, 2001) is satisfied. That is, when the sample size increases, the probability that the explanatory variables used in the model are correctly selected converges to 1 (consistency of variable selection). Also, the estimator for the explanatory variable has asymptotic normality.
  • the regression coefficient is subject to the above-mentioned code constraint, and the sparsification performance can be improved.
  • the parameter ⁇ for controlling the strength of regularization is a so-called hyperparameter. That is, in addition to the processing shown in FIG. 6, the optimum value of the coefficient is determined by a method using existing cross-validation.
  • the cost function shown in the following equation (12) is used instead of the cost function shown in the equation (2).
  • the regression equation is the same as that shown in equation (1).
  • ⁇ is a parameter for controlling the strength of regularization, and takes a value of zero or more.
  • the optimum value of ⁇ is determined by an existing method using cross-validation.
  • the regularization term ⁇ R SL (w) also imposes a sign constraint on either the positive or negative side. Specifically, restriction sign of x k in the table of FIG. 1 for positive, a value of R SL + (w), if the constraint code is negative, taking the value of R SL- (w). That is, in the regularization term, the coefficient w k increases the cost according to the sum of the absolute values of the coefficients in either the positive or negative interval according to the constraint sign, and infinitely increases the cost in the other interval. Make it big. In other words, not only the cost is set to infinity when it does not match the constraint code (that is, it corresponds to the case where ⁇ in Eq. (2) is set to infinity), but also when it matches the constraint code, ⁇ and w are set. Increase costs accordingly.
  • FIG. 11A and 11B are schematic diagrams for explaining the constraints imposed on the regression coefficient w.
  • the vertical axis represents ⁇ R SL + (w) and the horizontal axis represents w.
  • the coefficient w k is zero or more, the predicted value ⁇ by the regression equation increases as the input x k of the regression equation shown in the equation (1) increases.
  • the constraint code associated with x k is the case of the positive, the regularization term when the value of the input x k is also increased value of the predicted value ⁇ as increased small, the value of the input x k is increased
  • the cost function is defined so that the regularization term becomes larger as the value of the predicted value ⁇ decreases.
  • the vertical axis represents ⁇ R SL ⁇ (w) and the horizontal axis represents w.
  • FIG. 12 is a diagram showing the relationship between the parameter ⁇ and the correlation coefficient r.
  • the horizontal axis represents the parameter ⁇ and the vertical axis represents the correlation coefficient r.
  • the solid line represents the result of the method according to the present embodiment
  • the broken line represents the result of the existing L1 regularization (LASSO).
  • the correlation coefficient r was higher in the results of the method according to the present embodiment than in the existing LASSO results, especially in the range where ⁇ was smaller than 0.001.
  • Figure 13 is a diagram showing the relationship between the coefficient of determination R 2 and parameter beta.
  • the solid line represents the result of the method according to the present embodiment, and the broken line represents the result of the existing L1 regularization (LASSO).
  • LASSO L1 regularization
  • the coefficient of determination R 2 the result of the method according to the present embodiment was higher than the result of the existing LASSO, particularly in the range where ⁇ was smaller than 0.001.
  • FIG. 14 is a diagram showing the relationship between the parameter ⁇ and RMSE (Root Mean Square Error). In the line graph of FIG. 14, the horizontal axis represents the parameter ⁇ and the vertical axis represents the RMSE.
  • the solid line represents the result of the method according to the present embodiment
  • the broken line represents the result of the existing L1 regularization (LASSO).
  • LASSO existing L1 regularization
  • the regression coefficient can be determined even when the number of explanatory variables is larger than the number of training data, and the regression coefficient can be determined even when the number of explanatory variables is larger than the number of training data.
  • Prediction performance generalization performance
  • LASSO LASSO
  • the computer configuration shown in FIG. 5 is an example, and is not limited to such an example.
  • at least a part of the functions of the regression analysis device 1 may be distributed to a plurality of devices to be realized, or the same function may be provided by a plurality of devices in parallel.
  • at least a part of the functions of the regression analyzer 1 may be provided on the so-called cloud.
  • the regression analyzer 1 does not have to have a part of the configuration such as the verification processing unit 144.
  • Equation (2) The cost function shown in equation (2) is L1 regularized on one side of positive or negative, but it also works with L2 norm and other convex functions. That is, instead of the sum of the absolute values of the coefficients, a term that imposes a sum of squares of the coefficients or other penalties on one side of the positive or negative may be used.
  • the content of the data analyzed by the regression analyzer 1 is not particularly limited. In addition to the prediction of characteristic values such as quality in the manufacturing industry described in the examples, it can be applied to the non-manufacturing industry and various other fields.
  • the present disclosure also includes a method for executing the above-mentioned processing, a computer program, and a computer-readable recording medium on which the program is recorded.
  • the recording medium on which the program is recorded can perform the above-mentioned processing by causing the computer to execute the program.
  • the computer-readable recording medium means a recording medium in which information such as data and programs is stored by electrical, magnetic, optical, mechanical, or chemical action and can be read from a computer.
  • recording media those that can be removed from the computer include flexible disks, magneto-optical disks, optical disks, magnetic tapes, memory cards, and the like.
  • HDD High Density Digital
  • SSD Solid State Drive
  • ROM Read Only Memory
  • Regression analyzer 11 Communication I / F 12: Storage device 13: Input / output device 14: Processor 141: Data acquisition unit 142: Coefficient update unit 143: Convergence judgment unit 144: Verification processing unit 145: Operation processing unit

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)
  • Stored Programmes (AREA)

Abstract

説明変数の変動と目的変数の変動とに対応関係を有する回帰モデルを構築する。回帰分析装置は、回帰モデルの目的変数及び説明変数として用いられる訓練データと、目的変数を正又は負の方向に変動させるために、説明変数を正及び負のいずれに変動させるべきかを予め定義する制約条件とを格納する記憶装置から、訓練データ及び制約条件を読み出すデータ取得部と、制約条件に反する場合にコストを増大させる正則化項を含むコスト関数を最小化させるように、訓練データを用いて、回帰モデルにおける説明変数の係数を繰り返し更新する係数更新部とを備える。

Description

回帰分析装置、回帰分析方法及びプログラム
　本開示は、回帰分析装置、回帰分析方法及びプログラムに関する。
　従来、回帰モデルのパラメータを最小二乗法で推定するとき、例えばデータのサンプル数が少ないと最小二乗推定量が求められないという問題があった。そこで、Ｌ１ノルムと呼ばれる制約条件を与える手法が提案されていた（例えば、非特許文献１）。Ｌ１ノルムを制約条件とするパラメータ推定手法であるＬＡＳＳＯ（Least Absolute Shrinkage and Selection Operator）によれば、目的変数を説明するために適した説明変数の選択及び係数の決定が併せて行われる。
　また、ＬＡＳＳＯに関して、相関の高い説明変数を予めグループ化したり、クラスタリングしたりするような、様々な改良手法が提案されている。
Robert Tibshirani, "Regression Shrinkage and Selection via the Lasso", Journal of the Royal Statistical Society. Series B (Methodological) Vol. 58, No. 1 (1996), pp. 267-288
　従来、例えば所望の結果が得られるように制御を行う場合、予測モデルを用いて逆問題を解いても適切な結果が得られないことがあった。すなわち、予測モデルによる推定値を所望の値に近づけるために、説明変数の値をどのように変更すべきかがわからない。しかしながら、説明変数の組合せを変更してシミュレーションを繰り返す手法では計算コストがかかる。そこで、本技術は、説明変数の変動と目的変数の変動とに対応関係を有する回帰モデルを構築することを目的とする。
　回帰分析装置は、回帰モデルの目的変数及び説明変数として用いられる訓練データと、目的変数を正又は負の方向に変動させるために、説明変数を正及び負のいずれに変動させるべきかを予め定義する制約条件とを格納する記憶装置から、訓練データ及び制約条件を読み出すデータ取得部と、制約条件に反する場合にコストを増大させる正則化項を含むコスト関数を最小化させるように、訓練データを用いて、回帰モデルにおける説明変数の係数を繰り返し更新する係数更新部とを備える。
　上記のような正則化項により、制約条件に反するような係数は選択されず、目的変数を正又は負の方向に変動させるために、説明変数を正及び負のいずれに変動させればよいかがわかる回帰モデルを作成することができる。すなわち、説明変数の変動と目的変数の変動とに対応関係を有する回帰モデルを構築することができる。
　また、正則化項は、係数が、制約条件に応じた正又は負である区間において、係数の絶対値の和に応じてコストを増大させるようにしてもよい。例えば、係数が正又は負の片側において、Ｌ１正則化を用いた回帰モデルの構築を行ってもよい。また、正則化項は、係数が、制約条件に応じた正又は負である区間の一方において、係数の絶対値の和に応じてコストを増大させ、他方においてコストを無限大にするようにしてもよい。
　また、係数更新部は、係数が制約条件を満たす値に収束しない場合、係数をゼロにするようにしてもよい。このようにすれば、上述した制約条件の下で目的変数に寄与しない説明変数を回帰モデルから削除でき、スパースモデリングが実現される。
　また、係数更新部は、近接勾配法により前記係数を更新するようにしてもよい。このようにすれば、収束計算において、正則化項の微分不可能点を通過することが回避される。したがって、収束に要する時間を短縮することができる。
　なお、課題を解決するための手段に記載の内容は、本開示の課題や技術的思想を逸脱しない範囲で可能な限り組み合わせることができる。また、課題を解決するための手段の内容は、コンピュータ等の装置若しくは複数の装置を含むシステム、コンピュータが実行する方法、又はコンピュータに実行させるプログラムとして提供することができる。なお、プログラムを保持する記録媒体を提供するようにしてもよい。
　開示の技術によれば、説明変数の変動と目的変数の変動とに対応関係を有する回帰モデルを構築することができる。
図１は、回帰式の作成に用いる訓練データの一例を示す図である。 図２Ａは、回帰係数に課せられる制約を説明するための模式的な図である。 図２Ｂは、回帰係数に課せられる制約を説明するための模式的な図である。 図３は、パラメータｗの更新を説明するための図である。 図４は、パラメータηの更新を説明するための図である。 図５は、上述した回帰分析を行う回帰分析装置１の構成の一例を示すブロック図である。 図６は、回帰分析装置が実行する回帰分析処理の一例を示す処理フロー図である。 図７Ａは、制約の強さを表すパラメータαと相関係数ｒとの関係を示す図である。 図７Ｂは、制約の強さを表すパラメータαと相関係数ｒとの関係を示す図である。 図８は、制約の強さを表すパラメータαと決定係数Ｅとの関係を示す図である。 図９は、学習に用いるデータ数Ｔと相関係数ｒとの関係を示す図である。 図１０は、学習に用いるデータ数Ｔと決定係数Ｅとの関係を示す図である。 図１１Ａは、回帰係数に課せられる制約を説明するための模式的な図である。 図１１Ｂは、回帰係数に課せられる制約を説明するための模式的な図である。 図１２は、パラメータβと相関係数ｒとの関係を示す図である。 図１３は、パラメータβと決定係数Ｒ２との関係を示す図である。 図１４は、パラメータβとＲＭＳＥとの関係を示す図である。
　以下、図面を参照しつつ回帰分析装置の実施形態について説明する。
＜実施形態＞
　本実施形態に係る回帰分析装置は、１以上の説明変数（独立変数）と、１つの目的変数（従属変数）との関係を表す回帰式（回帰モデル）を構築する。このとき、説明変数の少なくともいずれかには、当該説明変数の変動の、正又は負の方向と、目的変数の変動の、正又は負の方向とに、一定の対応関係を有するような制約（「符号制約」と呼ぶ）を課して回帰式を作成する。
　図１は、回帰式の作成に用いる観測値（訓練データ）の一例を示す図である。図１の表は、Ｋ種類の入力ｘ（ｘ１～ｘＫ）の列と、出力ｙの列とを含む。入力ｘは説明変数に相当し、出力ｙは目的変数に相当する。また、訓練データの個々の標本であるデータポイントｔ（ｔ１～ｔＴ，・・・）を表す複数のレコードのうち、Ｔ個のレコードを用いて回帰式を作成するものとする。また、Ｋ種の入力ｘの少なくとも一部に対して正又は負の符号（本実施形態に係る制約条件を表す情報であり、「制約符号」と呼ぶものとする）が対応付けられているものとする。各入力ｘに対応付けられた制約符号は、構築する回帰式において、出力ｙを正の方向に変動させるために、当該入力ｘを正又は負のうちいずれの方向に変動させればよいかを予め定義するための情報である。
　回帰式は、例えば次の式（１）で表される。
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
以上のような処理によって、係数ｗを決定することができる。係数ｗは、符号制約を満たし、且つ目的変数に寄与する値に収束し、そのような値がなければ係数ｗはゼロに近づいてゆく。すなわち、符号制約を満たす値がない場合は、図２Ａ及び図２Ｂに示したように正則化による罰則効果がはたらいて符号制約に反する値を引き戻すことにより、結果的にゼロに収束してゆく。よって、いわゆるＬＡＳＳＯと同様に回帰係数の一部をゼロと推定し得る。
　なお、ηの値も、係数を更新する処理において繰り返される各ステップにおいて適宜更新するようにしてもよい。図４は、適切なηを探索するための模式的なコードの一例を示す。例えば図４に示すような処理が、各ステップにおいて実行される。η０は予め定められた初期値である。βは、例えば１より小さい正の値であり、ηを減少させるように更新する。このように係数ｗを更新するステップ幅であるηを調整することで、係数ｗを適切に収束させることができる。
＜装置構成＞
　図５は、上述した回帰分析を行う回帰分析装置１の構成の一例を示すブロック図である。回帰分析装置１は、一般的なコンピュータであり、通信インターフェース（Ｉ／Ｆ）１１と、記憶装置１２と、入出力装置１３と、プロセッサ１４とを備えている。通信Ｉ／Ｆ１１は、例えばネットワークカードや通信モジュールであってもよく、所定のプロトコルに基づき、他のコンピュータと通信を行う。記憶装置１２は、ＲＡＭ（Random Access Memory）やＲＯＭ（Read Only Memory）等の主記憶装置、及びＨＤＤ（Hard-Disk Drive）やＳＳＤ（Solid State Drive）、フラッシュメモリ等の補助記憶装置（二次記憶装置）であってもよい。主記憶装置は、プロセッサ１４が読み出すプログラムや当該プログラムが処理する情報を一時的に記憶する。補助記憶装置は、プロセッサ１４が実行するプログラムや当該プログラムが処理する情報等を記憶する。本実施形態では、記憶装置１２には、訓練データ及び制約条件を表す情報が、一時的に又は永続的に記憶されているものとする。入出力装置１３は、例えば、キーボード、マウス等の入力装置、モニタ等の出力装置、タッチパネルのような入出力装置等のユーザインターフェースである。プロセッサ１４は、ＣＰＵ（Central Processing Unit）等の演算処理装置であり、プログラムを実行することにより本実施形態に係る各処理を行う。図１の例では、プロセッサ１４内に機能ブロックを示している。すなわち、プロセッサ１４は、所定のプログラムを実行することにより、データ取得部１４１、係数更新部１４２、収束判定部１４３、検証処理部１４４及び運用処理部１４５として機能する。
　データ取得部１４１は、記憶装置１２から、訓練データ及び制約条件を表す情報を取得する。係数更新部１４２は、上述した制約条件の下で回帰式の係数を更新する。また、収束判定部１４３は、更新された係数の値が収束したか判定する。なお、収束していないと判定された場合、係数更新部１４２は、係数の更新を繰り返す。収束したと判定された場合、例えば係数更新部１４２は、最終的に生成される係数を記憶装置１２に記憶させる。また、検証処理部１４４は、所定の評価指標に基づいて作成された回帰式を評価する。運用処理部１４５は、作成された回帰式と例えば新たに取得される観測値とを用いて、予測値を算出する。また、運用処理部１４５は、作成された回帰式と任意の値とを用いて、条件を変更した場合の予測値を算出してもよい。ここで、任意の値は、例えば通信Ｉ／Ｆ１１又は入出力装置１３を介してユーザが入力する値であってもよい。本実施形態において作成される回帰式は、説明変数の変動の方向と目的変数の変動の方向とに一定の対応関係を有するため、例えば予測値を所望の値に近づけるために入力値を増加させればよいか減少させればよいか、ユーザは容易に推定できる。したがって、たとえば推定値に基づいて何らかの制御を行う場合に、本実施形態に係る回帰式は有効である。
　以上のような構成要素が、バス１５を介して接続されている。
＜回帰分析処理＞
　図６は、回帰分析装置が実行する回帰分析処理の一例を示す処理フロー図である。回帰分析装置１のデータ取得部１４１は、訓練データと制約条件を表す情報と記憶装置１２から読み出す（図６：Ｓ１１）。本ステップでは、例えば図１に示したような入力ｘ及び出力ｙの値が訓練データとして読み出される。なお、入力ｘを説明変数として扱い、出力ｙを目的変数として扱うものとする。また、図１において、入力ｘに対応付けて登録されている正又は負の符号が、制約条件を表す情報として読み出される。回帰分析装置１は、読み出される符号を、上述した制約符号として用いる。なお、本実施形態では、式（１）に示したような回帰式を用いる。
　また、回帰分析装置１の係数更新部１４２は、上述した符号制約の下で回帰係数を更新する（図６：Ｓ１２）。本ステップでは、係数更新部１４２は、例えば図３において上側の矢印で示したように、式（２）に示したコスト関数Ｅ（ｗ）を最小化するように係数ｗを更新する。具体的には、係数更新部１４２は、式（６）～式（８）に基づいて係数ｗを更新することができる。
　本実施形態に係るコスト関数Ｅ（ｗ）の正則化項は、Ｓ１１で取得した制約条件を満たさない場合にコストが増加するように定義されている。すなわち、正則化項は、説明変数の変動の、正又は負の方向と、目的変数の変動の、正又は負の方向とが、予め定められた対応関係を有するときにコスト関数Ｅ（ｗ）の値を減少させるものである。また、係数更新部１４３は、係数が制約条件を満たす値に収束しない場合、係数をゼロにする。
　また、回帰分析装置１の収束判定部１４３は、係数ｗが収束したか又は係数ｗがゼロにされたか判定する（図６：Ｓ１３）。本ステップでは、収束判定部１４３は、更新される係数ｗの勾配が充分ゼロに近づいた場合に収束したと判断する。具体的には、収束判定部１４３は、式（７）又は式（８）において係数ｗの値が変化しなくなったときに、収束したと判断する。
　係数ｗが収束しておらず、ゼロにされてもいないと判定された場合（Ｓ１３：ＮＯ）、Ｓ１２に戻って処理を繰り返す。一方、係数ｗが収束した、又はゼロにされたと判定された場合（Ｓ１３：ＹＥＳ）、収束判定部１４３は、回帰式を記憶装置１２に格納する（図６：Ｓ１４）。本ステップでは、収束判定部１４３は、更新後の係数ｗを記憶装置１２に記憶させる。
　また、回帰分析装置１の検証処理部１４４は、作成された回帰式の精度を検証するようにしてもよい（図６：Ｓ２０）。本ステップでは、検証処理部１４４は、例えば交差検証により、テストデータを用いて回帰式の精度を検証する。また、検証処理部１４４は、相関係数や所定の決定係数等、所定の評価指標に基づいて検証することができる。なお、後述するように、本ステップは省略してもよい。
　そして、回帰分析装置１の運用処理部１４５は、作成された回帰式を用いて、運用処理を行う（図６：Ｓ３０）。本ステップでは、運用処理部１４５は、例えば図１に示したデータ番号がｔＴ＋１のレコードのように、新たな入力ｘに対する出力ｙの予測値を算出する。なお、本ステップは、Ｓ１４で記憶された回帰式を用いて、回帰分析装置１以外の装置（図示せず）が行うようにしてもよい。
＜実施例＞
　生産プラントから得られるセンシングデータを用いて回帰式を構築し、精度を評価した。図１に示した入力及び出力の各々として、異なるセンサの出力値を用いた。また、センサから継続的に出力されるセンシングデータについて、直近のデータ数Ｔを学習区間とした。また、制約符号は、生産プラントに関する知見に基づいて予め設定された。
　評価指標として用いる相関係数ｒは、次の式（９）で求められる。
Figure JPOXMLDOC01-appb-M000009
すなわち、式（９）の分子は、予測値μと訓練データの実測値ｙとの共分散である。式（９）の分母は、予測値μの標準偏差と訓練データの実測値ｙの標準偏差との積である。
　また、他の評価指標として用いる決定係数Ｅは、次の式（１０）で求められる。
Figure JPOXMLDOC01-appb-M000010
決定係数Ｅは、観測値の分布に対する予測値の分布の大きさを表す値である。標準化により、観測値の分布と予測値の分布と一致する場合、Ｅ＝１となる。また、観測値の分布に対して予測値の分布が狭い場合、Ｅ＜１となる。そして、観測値の分布に対して予測値の分布が広い場合、Ｅ＞１となる。
　図７Ａ及び図７Ｂは、複数の手法で構築されたモデルについて、制約の強さを表すパラメータαと相関係数ｒとの関係を示す図である。図８は、複数の手法で構築されたモデルについて、制約の強さを表すパラメータαと決定係数Ｅとの関係を示す図である。図７Ａ及び図７Ｂの折れ線グラフは横軸がパラメータαを表し、縦軸が相関係数ｒを表す。図７Ａ及び図７Ｂは、横軸のスケールが異なる。また、図８の折れ線グラフは、横軸がαを表し、縦軸が決定係数Ｅを表す。実線は実施形態に開示の手法、破線は実施形態の符号制約の一部をランダムに選択して正負を逆にした比較例、一点鎖線はＬ１正則化（ＬＡＳＳＯ）、二点鎖線は正則化なしの各結果を表す。なお、各手法においてデータ数Ｔを４０としてモデルの構築を行った。また、上述の通り制約符号は生産プラントに関する知見に基づいて予め設定されたものであるが、一般的には不適切な設定を含み得るものである。比較例は、誤りのある符号制約をシミュレートしたものといえる。
　図７Ａ及び図７Ｂに示すように、相関係数ｒは、本開示の手法、比較例、ＬＡＳＳＯ、制約なしの順に値が高かった。また、図８に示すように、決定係数Ｅは、ＬＡＳＳＯ、本開示の手法及び比較例、制約なしの順に値が１に近かった。図７Ａ及び図７Ｂからもわかるように、一般的なＬＡＳＳＯでは、パラメータαを大きくし過ぎると精度は低下する。すなわち、ＬＡＳＳＯにおいてαはいわゆるハイパーパラメータであり、交差検証による調整が必要である。一方、本開示の手法によれば、αを充分大きくとることで精度を向上させることができた。これは、人手によるパラメータ調整を不要にし得るという効果がある。また、比較例のように符号制約をランダムに与えた場合、例えば相関係数ｒは、実施形態に係る手法よりも低下した。すなわち、実施形態に係る手法は、分析対象とするデータが、説明変数の変動と目的変数の変動とに一定の対応関係を有し、これに合致した符号制約を与えた場合に、特に当てはまりのよいモデルを作成し得るといえる。また、図７Ａおよび図７Ｂからわかるように、破線で示す、符号制約をランダムに与えた比較例の場合であっても、二点鎖線で示す正則なしの場合よりも相関係数が高い。このことは、一部の説明変数について適切でない符号制約が課されたとしても、依然当てはまりのよいモデルを作成し得ることを示す。現実には、説明変数の変動と目的変数の変動との対応関係に関する知識が必ずしも完全でない場合が往々にしてある。そのような場合においても、実施形態に係る手法によれば、正則化なしの場合よりも当てはまりのよいモデルを作成し得るという効果を発揮する。
　図９は、複数の手法で構築されたモデルについて、学習に用いるデータ数Ｔと相関係数ｒとの関係を示す図である。図１０は、複数の手法で構築されたモデルについて、学習に用いるデータ数Ｔと決定係数Ｅとの関係を示す図である。図９に示すように、例えばデータ数Ｔが４０以下の場合においては、相関係数ｒは、本開示の手法、比較例、ＬＡＳＳＯ、制約なしの順に値が高かった。また、図１０に示すように、決定係数Ｅは、ＬＡＳＳＯ、本開示の手法及び比較例、制約なしの順に値が１に近かった。このように、本開示の手法は、訓練データが比較的少ない場合において有効といえる。すなわち、データが充分に収集できていない場合や、予測モデルは時々刻々と変化するがデータのみからは観測できない状態の変化がある等の理由で直近のデータしか使えないような場合にも有用である。
＜効果＞
　本開示の手法によれば、説明変数の変動の、正又は負の方向と、目的変数の変動の、正又は負の方向とに、一定の対応関係を有するような制約を満たす回帰式を生成することができる。したがって、ユーザは、回帰式を用いて、予測値μを所望の値に近づけるために、入力ｘｋの値を正又は負のいずれに変動させればよいかがわかるようになる。また、図７Ａ及び図７Ｂを用いて説明したように、制約の強さを表すパラメータαの調整が不要になるという利点もある。また、図９及び図１０を用いて説明したように、本開示の手法は、訓練データが比較的少ない場合において特に有効である。
　以下、効果について補足する。ここで、式（２）の正則化項に関して、次のことがいえる。
Figure JPOXMLDOC01-appb-M000011
　そして、例えば制約符号が正（Ｒ＋（ｗ））のとき、式（２）のコスト関数Ｅ（ｗ）のｗｋに関する劣微分は次のように求められる。
Figure JPOXMLDOC01-appb-M000012
なお、ここでは複数の入力ｘkの間には相関がないものと仮定し、δkk’は単位行列を表すものとする。
　そして、ｗｋは以下のように求められる。
Figure JPOXMLDOC01-appb-M000013
また、これを解き直すと、次のように求められる。
Figure JPOXMLDOC01-appb-M000014
ここで、αが充分大きいとすれば、下段の場合を考慮せずに、ｗｋは次の式（１１）で表すことができる。
Figure JPOXMLDOC01-appb-M000015
　式（１１）の上段の場合は、最小二乗法と同じ解である。一方、一般的な最小二乗法においては符号制約が課されないため、例えばデータ数Ｔが比較的小さい場合においては、式（１１）の下段に相当するケースにおいても式（１１）の上段と同じ解が得られることがある。この場合、回帰式の出力を所望の値に近づけるために、説明変数の値をどのように変更すべきかがわからないことになる。一方、このような場合、本開示の技術によれば、式（１１）の下段に示すように係数ｗｋをゼロにする。すなわち、制約を満たすことができない説明変数ｘｋについては、作成される回帰式に用いられない。よって、説明変数の変動の、正又は負の方向と、目的変数の変動の、正又は負の方向とに、一定の対応関係を有するような制約を満たす回帰式を生成することができる。また、パラメータαの値は充分に大きな値とすることができ、調整は不要といえる。
　また、一般的なＬＡＳＳＯにおいては、例えば以下のようにｗｋが求められる。
Figure JPOXMLDOC01-appb-M000016
すなわち、本来収束すべき値からαだけ小さくなるようバイアスして推定される。このようなバイアスは、二乗誤差を大きくするように作用する。一方、本開示の技術によればこのようなバイアスは生じないため、回帰式の精度が向上するといえる。
　また、式（１１）によれば、オラクル性質（Oracle property, Fan and Li, 2001）が満たされる。すなわち、標本サイズが大きくなるとき、モデルに用いられる説明変数が正しく選択される確率が１に収束する（変数選択の一致性）。また、説明変数に対する推定量は漸近正規性を有する。
＜実施形態２＞
　本実施形態は、回帰係数に上述した符号制約を課すとともに、スパース化の性能を向上させることができる。また、正則化の強さを制御するためのパラメータβは、いわゆるハイパーパラメータとする。すなわち、図６に示した処理に加え、既存の交差検証を用いた手法により係数の最適値を決定する。本実施形態では、式（２）に示したコスト関数に代えて、次の式（１２）に示すコスト関数を用いる。なお、回帰式は式（１）に示したものと同じである。
Figure JPOXMLDOC01-appb-M000017
βは正則化の強さを制御するためのパラメータであり、ゼロ以上の値をとる。また、βは、交差検証を用いた既存の手法により最適値が決定される。本実施形態に係る正則化項βＲＳＬ（ｗ）も、正又は負の片側で符号制約を課す。具体的には、図１のテーブルにおいてｘｋの制約符号が正の場合、ＲＳＬ＋（ｗ）の値をとり、制約符号が負の場合、ＲＳＬ－（ｗ）の値をとる。すなわち、正則化項は、係数ｗｋが、制約符号に応じた正及び負のいずれか一方の区間において、係数の絶対値の和に応じてコストを増大させ、他方の区間においてはコストを無限大にする。換言すれば、制約符号と一致しない場合はコストを無限大にする（すなわち、式（２）のαを無限大にする場合に相当）だけでなく、制約符号と一致する場合もβ及びｗに応じてコストを増大させる。
　図１１Ａ及び図１１Ｂは、回帰係数ｗに課せられる制約を説明するための模式的な図である。図１１Ａのグラフは、縦軸がβＲＳＬ＋（ｗ）を、横軸がｗを表す。上述の式（１２）は、入力ｘｋに対応付けられた制約符号が正の場合であって、入力ｘｋの係数ｗｋがゼロ以上のときはＲＳＬ＋（ｗ）＝ｗであり、ｗの増加に応じてＥ（ｗ）も増加させる。一方、入力ｘｋの係数ｗｋがゼロ未満のときはＲＳＬ＋（ｗ）＝＋∞でありコストを正の無限大に発散させる。これは、図２Ａに示したαが十分に大きな値である場合に予測性能が最大化されることに基づいて、無限大としたものである。すなわち、本実施形態の正則化項は、制約符号と一致しない区間においてはコストを無限大にし、制約符号と一致する区間においても回帰係数ｗとパラメータβの大きさに応じてコストを増大させる。ここで、係数ｗｋがゼロ以上のときは、式（１）に示した回帰式の入力ｘｋが増加するほど回帰式による予測値μも増加する。すなわち、ｘｋに対応付けられた制約符号が正の場合は、入力ｘｋの値が増加するほど予測値μの値も増加するときに正則化項が小さく、入力ｘｋの値が増加するほど予測値μの値が減少するときに正則化項が大きくなるように、コスト関数が定義されている。
　図１１Ｂのグラフは、縦軸がβＲＳＬ－（ｗ）を、横軸がｗを表す。上述の式（１２）は、入力ｘｋに対応付けられた制約符号が負の場合であって、入力ｘｋの係数ｗｋがゼロ以上のときはＲＳＬ－（ｗ）＝＋∞でありコストを正の無限大に発散させる。これは、図２Ｂに示したαが十分に大きな値である場合に予測性能が最大化されることに基づいて無限大としたものであり、十分大きな値を意図したものである。一方、入力ｘｋの係数ｗｋがゼロ未満のときはＲＳＬ－（ｗ）＝－ｗであり、ｗの減少に応じてＥ（ｗ）を増加させる。ここで、係数ｗｋがゼロ未満のときは、式（１）に示した回帰式の入力ｘｋが増加するほど回帰式による予測値μは減少する。すなわち、ｘｋに対応付けられた制約符号が負の場合は、入力ｘｋの値が増加するほど予測値μの値が減少するときに正則化項が小さく、入力ｘｋの値が増加するほど予測値μの値が増加するときに正則化項が大きくなるように、コスト関数が定義されている。
＜効果＞
　Leave-one-out法による交差検証により、本実施形態に係る手法と、既存のＬ１正則化（ＬＡＳＳＯ）性能評価を行った。学習データ数Ｎは１０であり、特徴数Ｋは１１とした。図１２は、パラメータβと相関係数ｒとの関係を示す図である。図１２の折れ線グラフは横軸がパラメータβを表し、縦軸が相関係数ｒを表す。また、実線は、本実施形態に係る手法による結果を表し、破線は、既存のＬ１正則化（ＬＡＳＳＯ）による結果を表す。相関係数ｒは、特にβが０．００１よりも小さい範囲において、本実施形態に係る手法の結果の方が既存のＬＡＳＳＯの結果よりも高くなった。図１３は、パラメータβと決定係数Ｒ２との関係を示す図である。図１３の折れ線グラフは、横軸がパラメータβを表し、縦軸が決定係数Ｒ２を表す。また、実線は、本実施形態に係る手法による結果を表し、破線は、既存のＬ１正則化（ＬＡＳＳＯ）による結果を表す。決定係数Ｒ２も、特にβが０．００１よりも小さい範囲において、本実施形態に係る手法の結果の方が既存のＬＡＳＳＯの結果よりも高くなった。図１４は、パラメータβとＲＭＳＥ（Root Mean Square Error）との関係を示す図である。図１４の折れ線グラフは、横軸がパラメータβを表し、縦軸がＲＭＳＥを表す。また、実線は、本実施形態に係る手法による結果を表し、破線は、既存のＬ１正則化（ＬＡＳＳＯ）による結果を表す。ＲＭＳＥも、特にβが０．００１よりも小さい範囲において、本実施形態に係る手法の結果の方が既存のＬＡＳＳＯの結果よりも低くなった。一般的に、説明変数の数が学習データの数よりも大きい場合には、方程式の数が解くべき変数の数より少なくなるため、何らかの正則化を施さなければ回帰係数を一意に定めることができない。図１２から図１４に示すように、本実施形態に係る手法により正則化を行えば、説明変数の数が学習データの数よりも大きい場合にも回帰係数を決定することができ、さらに、既存のＬＡＳＳＯと比較して予測性能（汎化性能）を向上させることができる。
＜変形例＞
　各実施形態における各構成及びそれらの組み合わせ等は、一例であって、本発明の主旨から逸脱しない範囲内で、適宜、構成の付加、省略、置換、及びその他の変更が可能である。本開示は、実施形態によって限定されることはなく、クレームの範囲によってのみ限定される。また、本明細書に開示された各々の態様は、本明細書に開示された他のいかなる特徴とも組み合わせることができる。
　図５に示したコンピュータの構成は一例であり、このような例には限定されない。例えば、回帰分析装置１の機能の少なくとも一部は、複数の装置に分散して実現するようにしてもよいし、同一の機能を複数の装置が並列に提供するようにしてもよい。また、回帰分析装置１の機能の少なくとも一部は、いわゆるクラウド上に設けるようにしてもよい。また、回帰分析装置１は、例えば検証処理部１４４等、一部の構成を備えていなくてもよい。
　また、式（２）に示したコスト関数は、正又は負の片側でＬ１正則化を行うものとしたが、Ｌ２ノルムやその他の凸関数によっても動作する。すなわち、係数の絶対値の和に代えて、正又は負の片側で係数の二乗和やその他のペナルティを課す項を用いるようにしてもよい。
　また、回帰分析装置１によって分析されるデータの内容は、特に限定されない。実施例で述べた製造業における品質等の特性値の予測のほか、非製造業やその他の様々な分野に適用できる。
　また、本開示は、上述した処理を実行する方法やコンピュータプログラム、当該プログラムを記録した、コンピュータ読み取り可能な記録媒体を含む。当該プログラムが記録された記録媒体は、プログラムをコンピュータに実行させることにより、上述の処理が可能となる。
　ここで、コンピュータ読み取り可能な記録媒体とは、データやプログラム等の情報を電気的、磁気的、光学的、機械的、または化学的作用によって蓄積し、コンピュータから読み取ることができる記録媒体をいう。このような記録媒体のうちコンピュータから取り外し可能なものとしては、フレキシブルディスク、光磁気ディスク、光ディスク、磁気テープ、メモリカード等がある。また、コンピュータに固定された記録媒体としては、ＨＤＤやＳＳＤ（Solid State Drive）、ＲＯＭ等がある。
１：　回帰分析装置
１１：　通信Ｉ／Ｆ
１２：　記憶装置
１３：　入出力装置
１４：　プロセッサ
１４１：　データ取得部
１４２：　係数更新部
１４３：　収束判定部
１４４：　検証処理部
１４５：　運用処理部

Claims (6)

  1. 　回帰モデルの目的変数及び説明変数として用いられる訓練データと、前記目的変数を正又は負の方向に変動させるために、前記説明変数を正及び負のいずれに変動させるべきかを予め定義する制約条件とを格納する記憶装置から、前記訓練データ及び前記制約条件を読み出すデータ取得部と、
    　前記制約条件に反する場合にコストを増大させる正則化項を含むコスト関数を最小化させるように、前記訓練データを用いて、前記回帰モデルにおける説明変数の係数を繰り返し更新する係数更新部と
    　を備える回帰分析装置。
  2. 　前記正則化項は、前記係数が、前記制約条件に応じた正又は負である区間において、前記係数の絶対値の和に応じて前記コストを増大させる
    　請求項１に記載の回帰分析装置。
  3. 　前記係数更新部は、前記係数が前記制約条件を満たす値に収束しない場合、前記係数をゼロにする
    　請求項１又は２に記載の回帰分析装置。
  4. 　前記係数更新部は、近接勾配法により前記係数を更新する
    　請求項１から３のいずれか一項に記載の回帰分析装置。
  5. 　コンピュータが、
    　回帰モデルの目的変数及び説明変数として用いられる訓練データと、前記目的変数を正又は負の方向に変動させるために、前記説明変数を正及び負のいずれに変動させるべきかを予め定義する制約条件とを格納する記憶装置から、前記訓練データ及び前記制約条件を読み出し、
    　前記制約条件に反する場合にコストを増大させる正則化項を含むコスト関数を最小化させるように、前記訓練データを用いて、前記回帰モデルにおける説明変数の係数を繰り返し更新する
    　回帰分析方法。
  6. 　コンピュータに、
    　回帰モデルの目的変数及び説明変数として用いられる訓練データと、前記目的変数を正又は負の方向に変動させるために、前記説明変数を正及び負のいずれに変動させるべきかを予め定義する制約条件とを格納する記憶装置から、前記訓練データ及び前記制約条件を読み出させ、
    　前記制約条件に反する場合にコストを増大させる正則化項を含むコスト関数を最小化させるように、前記訓練データを用いて、前記回帰モデルにおける説明変数の係数を繰り返し更新させる
    　プログラム。
PCT/JP2021/004178 2020-02-04 2021-02-04 回帰分析装置、回帰分析方法及びプログラム WO2021157669A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
CN202180012527.4A CN115053216A (zh) 2020-02-04 2021-02-04 回归分析装置、回归分析方法以及程序
JP2021575868A JPWO2021157669A1 (ja) 2020-02-04 2021-02-04
EP21751487.6A EP4102420A4 (en) 2020-02-04 2021-02-04 REGRESSION ANALYSIS DEVICE, REGRESSION ANALYSIS METHOD AND PROGRAM
US17/797,141 US20230059056A1 (en) 2020-02-04 2021-02-04 Regression analysis device, regression analysis method, and program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2020017477 2020-02-04
JP2020-017477 2020-02-04

Publications (1)

Publication Number Publication Date
WO2021157669A1 true WO2021157669A1 (ja) 2021-08-12

Family

ID=77199336

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2021/004178 WO2021157669A1 (ja) 2020-02-04 2021-02-04 回帰分析装置、回帰分析方法及びプログラム

Country Status (5)

Country Link
US (1) US20230059056A1 (ja)
EP (1) EP4102420A4 (ja)
JP (1) JPWO2021157669A1 (ja)
CN (1) CN115053216A (ja)
WO (1) WO2021157669A1 (ja)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012106198A (ja) * 2010-11-18 2012-06-07 Toshiba Corp 生物学的廃水処理装置
JP2013059241A (ja) * 2011-09-09 2013-03-28 Toshiba Corp 蓄電池貸出容量決定装置および蓄電池貸出容量決定方法
JP2018109805A (ja) * 2016-12-28 2018-07-12 みずほ第一フィナンシャルテクノロジー株式会社 説明変数を選択する装置、方法及びプログラム
JP2019144779A (ja) * 2018-02-19 2019-08-29 日本電信電話株式会社 因果推定装置、因果推定方法、及びプログラム

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6069460B1 (ja) * 2015-10-30 2017-02-01 みずほ第一フィナンシャルテクノロジー株式会社 説明変数を選択する装置、方法及びプログラム

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012106198A (ja) * 2010-11-18 2012-06-07 Toshiba Corp 生物学的廃水処理装置
JP2013059241A (ja) * 2011-09-09 2013-03-28 Toshiba Corp 蓄電池貸出容量決定装置および蓄電池貸出容量決定方法
JP2018109805A (ja) * 2016-12-28 2018-07-12 みずほ第一フィナンシャルテクノロジー株式会社 説明変数を選択する装置、方法及びプログラム
JP2019144779A (ja) * 2018-02-19 2019-08-29 日本電信電話株式会社 因果推定装置、因果推定方法、及びプログラム

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GO SATO, MAKIKO KORESHIMA, TAKUYA OWA, HIROTAKA TAMURA, JUN OHKUBO: "Numerical experiments of QUBO formulation for ReLU-type functions", IEICE TECHNICAL REPORT, NC, vol. 118, no. 470 (NC2018-52), 12 April 2019 (2019-04-12), JP, pages 49 - 54, XP009538982 *
ROBERT TIBSHIRANI: "Regression Shrinkage and Selection via the Lasso", JOURNAL OF THE ROYAL STATISTICAL SOCIETY. SERIES B (METHODOLOGICAL, vol. 58, no. 1, 1996, pages 267 - 288
See also references of EP4102420A4

Also Published As

Publication number Publication date
JPWO2021157669A1 (ja) 2021-08-12
US20230059056A1 (en) 2023-02-23
EP4102420A4 (en) 2024-03-06
CN115053216A (zh) 2022-09-13
EP4102420A1 (en) 2022-12-14

Similar Documents

Publication Publication Date Title
TWI444844B (zh) 模擬參數校正技術
Shimoni et al. Benchmarking framework for performance-evaluation of causal inference analysis
Beintema et al. Nonlinear state-space identification using deep encoder networks
US20220067588A1 (en) Transforming a trained artificial intelligence model into a trustworthy artificial intelligence model
US9626631B2 (en) Analysis device, analysis method, and program
CN111542792A (zh) 诊断装置和诊断方法
US20210232957A1 (en) Relationship analysis device, relationship analysis method, and recording medium
US20120310619A1 (en) Fast function extraction
Dubbs Test set sizing via random matrix theory
WO2021157669A1 (ja) 回帰分析装置、回帰分析方法及びプログラム
LeKander et al. Empirical evaluation of gradient methods for matrix learning vector quantization
JP7152938B2 (ja) 機械学習モデル構築装置および機械学習モデル構築方法
CN116993548A (zh) 基于增量学习的LightGBM-SVM的教育培训机构信用评估方法及系统
CN110991660A (zh) 基于蝗虫优化的lssvm-arima模型的态势分析方法、系统和存储介质
Araujo et al. Hybrid intelligent design of morphological-rank-linear perceptrons for software development cost estimation
US20240020531A1 (en) System and Method for Transforming a Trained Artificial Intelligence Model Into a Trustworthy Artificial Intelligence Model
JP2019053589A (ja) 強化学習プログラム、強化学習方法、および強化学習装置
CN114970657A (zh) 异常度计算系统及方法
JP6929260B2 (ja) 時系列特徴抽出装置、時系列特徴抽出方法及びプログラム
Liu et al. Model pursuit and variable selection in the additive accelerated failure time model
Hosein et al. A successive quadratic approximation approach for tuning parameters in a previously proposed regression algorithm
Sinha et al. LFT representation of a class of nonlinear systems: A data-driven approach
JP7298870B2 (ja) 分子動力学データ解析装置及びプログラム
US20220245379A1 (en) Information processing apparatus, information processing method, and non-transitory storage medium
Heinänen A Method for automatic tuning of PID controller following Luus-Jaakola optimization

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: 21751487

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2021575868

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2021751487

Country of ref document: EP

Effective date: 20220905