WO2016030938A1 - シミュレーションシステム、及びシミュレーション方法 - Google Patents

シミュレーションシステム、及びシミュレーション方法 Download PDF

Info

Publication number
WO2016030938A1
WO2016030938A1 PCT/JP2014/072132 JP2014072132W WO2016030938A1 WO 2016030938 A1 WO2016030938 A1 WO 2016030938A1 JP 2014072132 W JP2014072132 W JP 2014072132W WO 2016030938 A1 WO2016030938 A1 WO 2016030938A1
Authority
WO
WIPO (PCT)
Prior art keywords
behavior
vector
random number
parameter
random
Prior art date
Application number
PCT/JP2014/072132
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/JP2014/072132 priority Critical patent/WO2016030938A1/ja
Priority to US15/502,303 priority patent/US11003810B2/en
Priority to JP2016545097A priority patent/JP6219528B2/ja
Publication of WO2016030938A1 publication Critical patent/WO2016030938A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F7/00—Methods or arrangements for processing data by operating upon the order or content of the data handled
    • G06F7/58—Random or pseudo-random number generators
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present invention relates to a simulation system and a simulation method.
  • the behavior of the system is modeled, and the computer simulates the behavior of the system based on the model and predicts the future behavior.
  • a certain physical phenomenon or a social phenomenon in which many elements are entangled cannot be described by a simple deterministic model, but is described as a complex and probabilistic model.
  • a deterministic system has one future behavior of the system under a given initial condition.
  • the future behavior of the system is distributed according to a certain probability distribution, and the computer needs to perform a simulation, for example, to obtain a statistic such as an expected value of the index of interest. There is.
  • the Monte Carlo simulation method is known as a method for a computer to predict and simulate the statistical behavior of a stochastic system.
  • the Monte Carlo simulation method is a method for estimating an expected value of an index of interest by repeatedly performing a simulation simulating a stochastic behavior existing in a target system using random numbers. At this time, by increasing the number of simulations, it is considered that the expected value of the index of interest approaches the true value according to the law of large numbers.
  • Patent Document 1 includes “a data acquisition step 0201 for acquiring a long-term forecast published by the Japan Meteorological Agency for meteorological physical quantities, a meteorological time series model creating step 0202 for creating a time series model according to meteorological physical quantities, and an estimation point in an estimation period.
  • a meteorological physical quantity simulation step 0203 for simulating the meteorological physical quantity for the number of simulations, a meteorological physical quantity estimation result output step 0204 for outputting the simulation result of the meteorological physical quantity as a meteorological physical quantity estimation result, and a meteorological time series model creation step 0202 are further included. It is composed of a time series model parameter estimation step 02021 and a weather time series model parameter correction step 02022 ”(see summary).
  • a probabilistic model that is a target of Monte Carlo simulation such as a meteorological physical quantity that is a simulation target of the technique described in Patent Document 1
  • a meteorological physical quantity that is a simulation target of the technique described in Patent Document 1
  • the sensitivity of the expected value of the target index to the parameter value is a differential coefficient that represents how much the expected value of the target index changes when the given parameter value changes slightly.
  • the sensitivity under the given parameter value can be optimized by a known optimization method such as a gradient method or a Newton method.
  • the sensitivity of the expected value of the target index to the parameter is simply the difference between the expected value of the target index under the given parameter value and the expected value of the target index when the parameter is changed slightly. Is divided by the minute change amount of the parameter.
  • a simulation system that repeats a simulation of a stochastic system with each of a plurality of random number vectors as input to calculate a plurality of realization values of the behavior of the stochastic system, comprising a processor and a memory, the memory Is a behavior function for determining an actual value of the behavior of a stochastic system based on a random number vector consisting of one or more random numbers and a parameter vector consisting of one or more parameters, and realizing the behavior of the stochastic system
  • a behavior evaluation function that determines an evaluation value of the behavior of the stochastic system based on the value, and a condition that keeps the realization value of the behavior function corresponding to the random number vector constant.
  • the random vector for the parameter based on the first differential coefficient
  • the processor acquires a plurality of random number vectors and a parameter vector, and each of the acquired plurality of random number vectors, the acquired parameter vector, and the behavior And a realization value of the behavior of the stochastic system corresponding to each of the plurality of acquired random number vectors, and each of the acquired plurality of random number vectors, the acquired parameter vector, and the weight A weight of each of the acquired random number vectors for each parameter of the acquired parameter vector based on the function, and the plurality of acquired multiple values based on each of the calculated realization values and the behavior evaluation function
  • the evaluation value of the behavior of the stochastic system corresponding to each of the random number vectors of Based on the calculated evaluation value of the behavior corresponding to each of the acquired plurality of random number vectors and the weight of each of the acquired random number vectors for the parameter selected from the acquired parameter vector, for the selected parameter,
  • the simulation system which calculates the sensitivity of the expected value of
  • the sensitivity to an expected value parameter of an arbitrarily determined attention index can be obtained from a result obtained by a set of Monte Carlo simulations. .
  • FIG. 1 is a block diagram illustrating a configuration example of a Monte Carlo simulation system in Embodiment 1.
  • FIG. 3 is a block diagram illustrating a configuration example of a random number importance calculation unit according to the first embodiment.
  • FIG. 3 is a flowchart illustrating an example of an entire process of Monte Carlo simulation according to the first exemplary embodiment.
  • 3 is a flowchart illustrating an example of a simulation execution process in the first embodiment.
  • 6 is a flowchart illustrating another example of the simulation execution process in the first embodiment.
  • 2 is an example of a user interface of the Monte Carlo simulation system in the first embodiment.
  • 6 is an example of a circuit of an SRAM cell composed of six transistors in the second embodiment. It is an example of the balance plan table of a bank in Example 3.
  • Example 4 it is a figure which shows the example of the concept which calculates a boundary residue in case a random vector follows a two-dimensional uniform distribution.
  • This example describes a Monte Carlo simulation system. First, a theoretical aspect for executing the simulation and calculating the sensitivity by the Monte Carlo simulation system of the present embodiment will be described.
  • the random variable X represents the behavior of the stochastic system that is the simulation target by the Monte Carlo simulation system of the present embodiment. That is, the behavior of the stochastic system to be simulated is modeled by the random variable X.
  • the stochastic system to be simulated is also simply referred to as a system.
  • x represents an actual value of the system behavior X.
  • the system behavior X may be either a one-dimensional or multi-dimensional vector, and an infinite number of random variables (discrete-time stochastic processes and distinctions) distinguished by a variable t called “time”. May be called). Therefore, when the system behavior X is a discrete-time stochastic process, X should originally be expressed as X (t), but is also simply expressed as X unless otherwise required.
  • a value x (t 1 ) indicating the behavior of the system at a certain time t 1 is referred to as a system state at the time t 1 .
  • the distribution of the system behavior X depends on the system parameter z, and the probability density function of X is expressed in the form of f (x, z).
  • the parameter z of the system is simply expressed as parameter z.
  • a (x) is called an evaluation function of system behavior.
  • the system behavior evaluation function is an example of an attention index.
  • the Monte Carlo simulation system generates many real values x of the behavior of the system by Monte Carlo simulation, and simulates the distribution of X, that is, f (x, z), so that the system evaluation function A ( The value of x) can be calculated approximately.
  • Each of the generated values x of the system behavior is also called a path.
  • the Monte Carlo simulation includes a calculation for calculating an actual value x of the behavior of the system from random number (s) whose distributions are known. This is formulated below.
  • the W (simultaneous) probability density function is written as g (w).
  • g (w) is known and random numbers according to this distribution can be easily generated on a computer.
  • w is a vector composed of M independent random numbers according to a uniform distribution on [0, 1] or a standard normal distribution N (0, 1).
  • w is a vector composed of M independent random numbers according to a uniform distribution on [0, 1] or a standard normal distribution N (0, 1).
  • ⁇ represents the distribution range of the random number vector w, that is, the support (unit) of g (w).
  • One of the objects of the present embodiment is to provide a method for efficiently calculating a differential coefficient obtained by differentiating a system evaluation function by each of N parameters (z i ) expressed by the following equation (3). .
  • the differential coefficient is called sensitivity to the parameter of the system evaluation function.
  • the basic idea of this embodiment is to continue to use the behavior (realized value) x of the system generated under a given set of parameters z even after slightly changing the parameter z i .
  • x (w, z) is continuous.
  • a minute change ⁇ z i that is equal to or smaller than a predetermined threshold value of a parameter
  • a minute value that is equal to or smaller than a predetermined threshold value of a random number vector that cancels it.
  • a (z + ⁇ z i ) cannot be calculated by a simple average such as Equation 2.
  • a (x (w, z)) is weighted appropriately and averaged to obtain A (z + ⁇ z i ) Can be calculated is the basic idea of this embodiment.
  • the technique of the present embodiment can be considered as an example of an importance sampling technique, and the weight corresponds to the importance.
  • Equation 5 means that each element of the vector ( ⁇ w / ⁇ z i )
  • x const is represented on the right side. However, as described above, since x can be an infinite dimensional variable, this notation is formal. This point will be described later.
  • Equation 7 the differential coefficient for the parameter of the system evaluation function, that is, the sensitivity is calculated based on the expected value of the evaluation value of the behavior of the system weighted by the equivalent differential importance.
  • R i in Equation 7 is a term representing a subtraction caused by setting the integration range of Equation 7 to ⁇ instead of ⁇ ′, and is expressed as Equation 9 below.
  • R i is called the boundary residue.
  • ⁇ represents the boundary (surface) of ⁇
  • n represents the outward unit normal vector of the boundary ⁇ .
  • h i (w, z) is expressed by the following formula 12.
  • the evaluation function a (x) of the system behavior can be freely set after the completion of the Monte Carlo simulation. For example, when a complex system is targeted for simulation, it is often known that what should be evaluated and what phenomenon should be noticed only after actually simulating the behavior of the system. Features are desirable.
  • Equations 13 and 14 are in the form of normal expected value calculation, the variance of the estimated value of the differential coefficient by L simulations decreases by O (L ⁇ 1 ).
  • the variance of the estimated value of the derivative by the method of simply restarting the Monte Carlo simulation by the number of parameters is reduced only by O (L ⁇ 1/4 ) or O (L ⁇ 1/3 ) with respect to the number of simulations. do not do. That is, the convergence speed of the variance of the estimated value of the differential coefficient by the simulation of this embodiment is fast.
  • the Monte Carlo simulation system since it is not necessary to calculate the boundary residual R i, it is possible to reduce the calculation time and computing resources.
  • x can also be an infinite dimensional vector, that is, the dimension of x can be greater than the dimension M of w.
  • Equation 4 since the number of constraint conditions is larger than the number of variables ⁇ w (number of dimensions), ⁇ w that satisfies Equation 4 generally does not exist.
  • ⁇ w the number of variables satisfying Equation 4 for a given parameter slight change ⁇ z i , the method of this embodiment cannot be applied.
  • time series x (t) representing the behavior of the system can be expressed as the following Expression 15, for example, the method of the present embodiment is applicable.
  • Equation 15 indicates that when the value y (w, z) determined by the random number vector w and the parameter z is determined, the system state x (t) at the time t is determined from the times t and y. In this case, obviously, the following equation 16 holds.
  • time series x (t) representing the behavior of the system is described by a recurrence formula such as the following Expression 17, for example, the method of the present embodiment is applicable.
  • Equation 17 shows that the system state x (t + 1) at time t + 1 is newly generated at time t, the system behavior in the past (before time t), and at time t, and a random number w j (w j may be a random number vector including a plurality of random numbers) and a parameter z.
  • w j may be a random number vector including a plurality of random numbers
  • FIG. 1 shows a configuration example of a Monte Carlo simulation system 100 of the present embodiment.
  • the Monte Carlo simulation system 100 is configured on a computer including a CPU 110, a memory 120, an input / output interface 130, and a secondary storage device 140, for example.
  • the Monte Carlo simulation system is composed of one computer, but may be composed of a plurality of computers.
  • the CPU 110 includes a processor and / or a logic circuit that operates in accordance with a program, inputs / outputs data, reads / writes, and executes each program described below.
  • the memory 120 temporarily loads and stores a program and data executed by the CPU 110, and further holds each program and each data.
  • the input / output interface 130 is an interface for inputting data from an external device or the like and outputting data to the external device or the like.
  • the secondary storage device 140 is a non-volatile storage device, such as an HDD. Programs, data, and the like may be stored in the secondary storage device 140.
  • the program is executed by the CPU 110 and performs a predetermined process using the memory 120 and the input / output interface 130. Therefore, in the present embodiment and the other embodiments, the description with the program as the subject may be the description with the CPU 110 as the subject. Alternatively, the process executed by the program is a process performed by a computer and a computer system on which the program operates.
  • the CPU 110 operates as a functional unit that realizes a predetermined function by operating according to a program.
  • the CPU 110 functions as a random number generation unit by operating in accordance with a random number generation unit 210 described later, and functions as a model behavior update unit by operating in accordance with a model behavior update unit 220 described later.
  • the CPU 110 also operates as a functional unit that realizes each of a plurality of processes executed by each program.
  • a computer and a computer system are an apparatus and a system including these functional units.
  • the program may be stored in the secondary storage device 140.
  • the memory 120 includes one or more scenario simulators 200, a storage unit 300, a post-processing calculation unit 400, and a Monte Carlo simulation overall control unit 500.
  • the scenario simulator 200 performs L times Monte Carlo simulation using different random number vectors.
  • the storage unit 300 stores numerical values output from the scenario simulator 200.
  • the post-processing calculation unit 400 calculates an expected value and sensitivity after the simulation is completed.
  • the Monte Carlo simulation overall control unit 500 controls the entire Monte Carlo simulation system 100. Further, the Monte Carlo simulation overall control unit 500 may calculate and output the calculation accuracy of the system evaluation function and the estimated sensitivity value calculated by the number of simulations specified by the user or the like. In addition, the Monte Carlo simulation overall control unit 500 may calculate and output the number of simulations and the like necessary for guaranteeing the calculation accuracy of the system evaluation function and the estimated sensitivity value specified by the user or the like. The Monte Carlo simulation overall control unit 500 activates the scenario simulator 200, for example, for the designated or calculated number of simulations, for example, L times, and activates the post-processing calculation unit 400 when the L simulations are completed.
  • the scenario simulator 200 inputs a one-dimensional or multi-dimensional parameter vector z of a simulation model given by a user or the like, and executes simulation of the model while internally generating a random number vector necessary for the simulation. Further, the scenario simulator 200 outputs an evaluation function a (k) of the system behavior and an equivalent differential importance h i (k) necessary for sensitivity calculation.
  • the scenario simulator 200 includes a random number generation unit 210, a model behavior update unit 220, a random number importance calculation unit 230, an index calculation unit 260, an equivalent differential importance calculation unit 270, and a simulation control unit 280, which are programs.
  • the scenario simulator 200 includes a model behavior storage unit 240 and a random number importance storage unit 250, each of which stores data.
  • the random number generation unit 210 generates a random number vector according to a probability density function corresponding to a simulation model that is predetermined or specified by a user or the like.
  • the model behavior update unit 220 receives a parameter given by a user or the like and a random number vector generated by the random number generation unit 210 as input, and the input parameter and the random number vector are predetermined or designated by the user or the like. Applies to simulation models (functions) to simulate real values of system behavior.
  • the random number importance calculation unit 230 receives the parameters and the random number vector generated by the random number generation unit 210 as input, and for each input parameter, the random importance for each parameter of w j (k) , that is, the summation symbol of Equation 8 Calculate the term in.
  • the model behavior storage unit 240 stores the realized value of the system behavior calculated by the model behavior update unit 220.
  • the random number importance storage unit 250 stores the random number importance calculated by the random number importance calculation unit 230.
  • the index calculation unit 260 calculates an evaluation function of the system behavior by, for example, calculating the average of the realized values of the system behavior stored in the model behavior storage unit 240.
  • the equivalent differential importance calculation unit 270 calculates the equivalent differential importance by calculating the sum of the random number importance stored in the random number importance storage unit 250, for example.
  • the simulation control unit 280 controls the overall operation of the scenario simulator 200.
  • the storage unit 300 stores an index value storage unit 310 that stores a system behavior evaluation function a (k) calculated by the index calculation unit 260, and an equivalent storage that stores the equivalent differential importance calculated by the equivalent differential importance calculation unit 270.
  • a differential importance storage unit 320 stores the equivalent differential importance calculated by the equivalent differential importance calculation unit 270.
  • the post-processing calculation unit 400 includes an expected value calculation unit 410 and a sensitivity calculation unit 420, which are programs.
  • the expected value calculation unit 410 calculates an average value of the evaluation functions a (k) of the behaviors of the L systems stored in the index value storage unit 310, and calculates the average value as an expected value, that is, the system evaluation function A (z ).
  • the expected value calculation unit 410 may further calculate the variance of the evaluation function a (k) of the behavior of the L systems, and output the accuracy (confidence interval) of the estimated value of A (z) at the same time. .
  • the sensitivity calculation unit 420 includes an evaluation function a (k) of the behavior of the system stored in the index value storage unit 310 and an equivalent differential importance h i for each of the N parameters stored in the equivalent differential importance storage unit 320. For the group (k) , the product of the corresponding values is calculated. Further, the sensitivity calculation unit 420 calculates an average value of the calculated products and outputs the average value of the system evaluation function A (z) for each of the N parameters as sensitivity ⁇ A / ⁇ z i . At this time, the sensitivity calculation unit 420 may further calculate the variance of the product of a (k) and h i (k ) and simultaneously output the accuracy (confidence interval) of the estimated value of sensitivity ⁇ A / ⁇ z i. Good.
  • FIG. 2 shows a configuration example of the random number importance calculation unit 230.
  • the random number importance calculation unit 230 includes a sample path fixed differential calculation unit 231, a sample path fixed second-order differential calculation unit 232, and a random number importance calculation unit 233, all of which are programs.
  • the sample path fixed differential calculation unit 231 receives the parameter and the random number vector generated by the random number generation unit 210 and calculates ( ⁇ w j / ⁇ z i )
  • x const .
  • x const ).
  • the sample path fixed second-order derivative calculation unit 231 receives ( ⁇ w j / ⁇ z i )
  • x const calculated by the sample path fixed derivative calculation unit 231 and the parameter as ⁇ ⁇ D i, j / ⁇ w j may be calculated.
  • the sample path fixed second-order differential calculation unit 232 may calculate ( ⁇ w j / ⁇ z i )
  • x const , and at this time, the sample path fixed differential calculation unit 231 has ( ⁇ w j / ⁇ z i )
  • It is not necessary to calculate x const .
  • the random number importance calculation unit 233 calculates the random number vector input from the random number generation unit 210 and the sample path fixed differential calculation unit 231 or the sample path fixed second-order differential calculation unit 232 ( ⁇ w j / ⁇ z i )
  • With x const and ⁇ D i, j / ⁇ w j calculated by the sample path fixed second-order derivative calculation unit 232 as an input, for each parameter, the random importance for each parameter of w j (k) , that is, The term in the summation symbol of Equation 8 is calculated.
  • FIG. 3 shows an example of the entire processing of the Monte Carlo simulation in the present embodiment.
  • the overall processing of the Monte Carlo simulation includes a simulation execution process in which the scenario simulator 200 performs L simulations using different random number vectors, and an expected value / sensitivity calculation process performed by the post-processing calculation unit 400 after the simulation ends. .
  • the scenario simulator 200 performs a simulation execution process (S301).
  • the simulation execution process is a process of executing L times of simulation that constitutes a set of Monte Carlo simulations.
  • the scenario simulator 200 executes individual simulations using different random number vectors.
  • the Monte Carlo simulation system of the present embodiment uses N system behavior evaluation values a (k) obtained from the system behavior obtained by the k-th simulation in addition to the system behavior evaluation function a (k) .
  • Equivalent differential importance h i (k) is calculated for each parameter. Since each of the L simulations in the simulation execution process is independent of each other, the L simulations are not sequentially performed. For example, the simulations are executed in parallel using a large number of scenario simulators 200 or a large number of computers. May be.
  • a series of functions that define a process for calculating the equivalent differential importance h i (k) is called a weight function.
  • the post-processing calculation unit 400 performs expected value / sensitivity calculation processing (S302).
  • the post-processing calculation unit 400 performs the expected value of the target index and the expected value parameter of the target index. This is a process for calculating an estimated value of sensitivity.
  • the expected value of the target index that is, the estimated value of the system evaluation function A (z) is calculated by calculating the average of the L system behavior evaluation functions a (k) stored in the index value storage unit. It is obtained by calculating.
  • the Monte Carlo simulation system in the present embodiment calculates an estimated value of sensitivity for each of the N parameters of the expected value of the attention index.
  • the sensitivity calculation unit 420 calculates the sensitivity estimation value for each i-th parameter of the system evaluation function A (z), and the corresponding system behavior evaluation function a ( stored in the index value storage unit 310 ). k) is multiplied by the equivalent differential importance h i (k) for the i-th parameter stored in the equivalent differential importance storage unit 320 to obtain an average of weights.
  • the Monte Carlo simulation system may perform, for example, the equivalent differential importance calculation process in step S301 and the sensitivity estimation value calculation process in step S302 only for parameters specified by the user or the like.
  • FIG. 3 and FIG. 4 show the details of the k-th simulation execution process when the time series x (t) representing the behavior of the system is calculated by the model expressed in the form of Equations 15 and 17, respectively. Indicates the processing contents.
  • FIG. 4 shows that when a time series x (t) representing the behavior of the system is expressed in the form of Formula 15, that is, when a value y (w, z) determined by a random number vector w and a parameter z is determined, at time t.
  • a time series x (t) representing the behavior of the system is expressed in the form of Formula 15, that is, when a value y (w, z) determined by a random number vector w and a parameter z is determined, at time t.
  • An example of detailed processing contents of the simulation execution processing is shown for the case where the system state x (t) is determined deterministically only at times t and y.
  • a scenario simulator 200 generates a random number vector w (k), calculates the y (k) from the w (k) a parameter z (S401). Specifically, the random number generation unit 210 generates an M-dimensional random number vector w (k) according to a predetermined probability density function g (w), and uses the generated random number as a model behavior update unit 220 and a random number importance level. It transmits to the calculation part 230. Note that the random number generation unit 210 may generate a random number vector using a known means such as a known pseudo-random number generation algorithm or a physical random number generation device. The model behavior updating unit 220 calculates y (k) by substituting the received w (k) and the given N-dimensional parameter value z into a given relational expression y (w, z).
  • the scenario simulator 200 calculates the equivalent differential importance h i (k) for each of the N parameters (S402). This is a process necessary for estimating the sensitivity to the parameter of the expected value of the attention index, and the process defined by the weight function includes this process. Note that the scenario simulator 200 may perform the process in step S402 only for parameters specified by the user or the like. Specifically, the equivalent differential importance h i (k) for the i-th parameter is calculated by the following four steps, that is, step S402 includes the following four steps.
  • the first step is a relationship D i, of the small change of the j-th element of the random vector w (k) and the i-th parameter such that the sample path fixed differential calculation unit 231 keeps the system behavior x constant .
  • j (k) ( ⁇ w j (k) / ⁇ z i )
  • x const is a step of calculating.
  • the subscript j represents the calculation for the j-th element of the M-dimensional random number vector w (k) .
  • the sample path fixed differential calculation unit 231 has ( ⁇ w j / ⁇ z i )
  • (y ⁇ (K) const) may be obtained. That is, when the sample path fixed differential calculation unit 231 is given by a user or the like, for example, the function ⁇ ( ⁇ y (k) / ⁇ z i ) / ( ⁇ y (k) / ⁇ w j (k) ), D i, j (k) may be calculated by substituting the input parameter z i and the input random number vector w j (k) into the function. The sample path fixed differential calculation unit 231 can calculate D i, j (k) from the function at high speed.
  • the sample path fixed differential calculation unit 231 calculates ⁇ w / ⁇ z i by numerically searching for ⁇ w that keeps y constant, for example, for a minute change ⁇ z i of the parameter, w j / ⁇ z i )
  • (y ⁇ (k) const) , that is, D i, j (k) may be calculated approximately.
  • the sample path fixed differential calculation unit 231 can calculate D i, j (k) by the above method.
  • the sample path fixed differential calculation unit 231 transmits the calculated D i, j (k) to the sample path fixed second-order differential calculation unit 232 and the random number importance calculation unit 233.
  • the second step is a step in which the sample path fixed second-order differential calculation unit 232 differentiates the received D i, j (k) by w j (k) .
  • D i, j and (k) referred to the results obtained by differentiating with w j (k) DD i, and j (k).
  • Sample paths fixed second-order derivative calculation unit 232 for example, received D i, by differentiating j (k) of numerically, DD i, may be calculated j (k).
  • the sample path fixed second-order differential calculation unit 232 for example, when a function for calculating DD i, j (k) is given, the input parameter z i and the input random vector w j (k).
  • DD i, j (k) may be calculated by substituting for this function.
  • the sample path fixed second-order differential calculation unit 232 transmits the calculated DD i, j (k) to the random number importance calculation unit 233.
  • the third step is a step in which the random number importance calculation unit 233 calculates the equivalent differential importance H i, j (k) regarding the j-th element w j (k) of the random number vector w (k) .
  • the third step corresponds to calculating the term in the summation symbol in equation 8.
  • the randomness importance calculation unit 233 is unique from the received D i, j (k) , the received DD i, j (k), and a predetermined probability density function g (w).
  • a value LG j (w j (k) ) (1 / g (w j (k) )) ⁇ ( ⁇ g / ⁇ w j (k) ) obtained by substituting the random number vector w j (k) into the function LG determined by And the random number importance H i, j (k) may be calculated by substituting these into the relational expression shown in the figure.
  • the random number importance calculation unit 233 stores the calculated H i, j (k) in the random number importance storage unit 250.
  • the calculation in the third step is performed when the M-dimensional w is a uniformly distributed random vector or when w is an M-dimensional independent standard normal random vector. It is expressed as Equation 11 or Equation 12. Therefore, for example, when the random number w j (k) follows a uniform distribution, since the random number importance is represented by DD i, j (k) from the expression 11, the random number importance calculating unit 230 performs the first step. It is possible to calculate the importance of random numbers by performing only.
  • the fourth step is a step in which the equivalent differential importance calculation unit 270 calculates an equivalent differential importance h i (k) for the i-th parameter.
  • the fourth step corresponds to the sum calculation in Equation 8 above.
  • the equivalent differential importance calculation unit 270 acquires the random importance H i, j (k) stored in the random importance storage unit 250 and calculates the sum of j. Since the calculation of the equivalent differential importance h i (k) for each parameter z i is independent from each other, the equivalent differential importance of N parameters z 1 ,. Rather than performing the calculations sequentially, for example, the calculations may be performed in parallel using a large number of computers.
  • the model behavior updating unit 220 actually simulates and calculates the system behavior x (k) (t) (S403).
  • the model behavior update unit 220 can perform the process by sequentially calculating the state of the system at the time t while increasing the time t from the initial time toward the end time. Other known efficient simulation calculation methods may be used.
  • the model behavior update unit 220 stores the calculated x (k) (t) in the model behavior storage unit 240.
  • the process in step S403 may be performed in parallel with the process in step S402, for example, or may be performed before the process in step S402.
  • the index calculation unit 260 acquires the system behavior x (k) (t) from the model behavior storage unit 240 and calculates the evaluation value of the target index, that is, the evaluation function a (k) of the system behavior. (S404). The index calculation unit 260 stores the calculated a (k) in the index value storage unit 310.
  • the scenario simulator 200 performs the system behavior evaluation function a (k) , by the processing shown in FIG. And the equivalent differential importance h i (k) for each of the N parameters.
  • FIG. 5 shows the case where the time series x (t) representing the behavior of the system is expressed in the form of Expression 17, that is, the state x (t + 1) of the system at time t + 1 is the time t and the past (before time t). ) Of the simulation execution process when determined by the state of the system and the random number w j (which may be different for each time) newly generated at time t (which may be a random number vector including a plurality of random numbers) and the parameter z Detailed processing contents are shown.
  • the scenario simulator 200 increases the time t from the simulation start time t_start to the simulation end time t_end by a simulation time step ⁇ t, and the system state x (k) (t) at time t and N parameters.
  • the simulation is performed by repeatedly performing the process of calculating the random importance H i, j (k) for the j-th random number w j (k) for all.
  • the scenario simulator 200 uses the time series of the system behavior x (k) (t) obtained as a result of the simulation to evaluate the evaluation function a (k) of the system behavior and each of the N parameters. Equivalent differential importance h i (k) is calculated.
  • the simulation control unit 280 sets the time t to the simulation start time t_start, and notifies each part of the scenario simulator 200 of the start of the loop processing related to the time t (S501).
  • the loop processing related to time t includes processing (S502) in which the model behavior update unit 220 calculates the system state x (k) (t) at time t, and random number importance calculation unit 230 generates random numbers generated at time t. Processing (S503) for calculating the randomness importance H i, j (k) for each of the N parameters.
  • the random number generation unit 210 generates a random number w j (k) necessary for calculating the system state x (k) (t) at time t according to a predetermined probability density function g (w).
  • the generated random number is transmitted to the model behavior update unit 220 and the randomness importance calculation unit 230.
  • the model behavior updating unit 220 substitutes the received random number w j (k) , the past system behavior (before time t), the parameter z, and the time t into the given recurrence formula (Equation 17) to obtain the time Compute the state x (k) (t) of the system at t.
  • the random number importance calculation unit 230 calculates, for each of the N parameters z 1 ,..., Z N , the random number importance H i, j (k) related to the parameter of the random number w j (k) generated at time t. Calculate (S503). Note that the scenario simulator 200 may perform the process in step S503 only for parameters designated by the user or the like.
  • step S503 includes the following three steps. Since the calculation of the equivalent differential importance h i (k) for each parameter is independent of each other, the calculation of the equivalent differential importance for N parameters z 1 ,..., Z N is performed by one scenario simulator 200. Instead of sequentially performing the calculation, for example, the calculation may be simultaneously performed in parallel using a large number of computers.
  • the relationship D i, j (k) the random number w j (k) and the minute change of the i-th parameter such that the sample path fixed differential calculation unit 231 keeps the system behavior x constant.
  • the sample path fixed differential calculation unit 231 receives the input parameter z i and the input random number. D i, j (k) may be calculated by substituting w j (k) into the function. In addition, the sample path fixed differential calculation unit 231 calculates ⁇ w / ⁇ z i by numerically searching for ⁇ w that keeps f constant, for example, for a minute change ⁇ z i of the parameter, thereby calculating ⁇ w / ⁇ z i.
  • the sample path fixed differential calculation unit 231 transmits the calculated D i, j (k) to the sample path fixed second-order differential calculation unit 232 and the random number importance calculation unit 233.
  • the second step is a step of calculating the DD i, j (k) .
  • the second step is the same as the second step in step S402.
  • the third step is a step in which the random number importance calculation unit 233 calculates the random number importance H i, j (k) related to the random number w j (k) at the time t.
  • the third step is the same as the third step in step S402.
  • the simulation control unit 280 sets the time to t + ⁇ t (S504).
  • the time t + ⁇ t is before the simulation end time t_end (S505: Yes)
  • the Monte Carlo simulation overall control unit 500 restarts the scenario simulator 200, and the scenario simulator 200 performs the processes of steps S502 to S504. .
  • the index calculation unit 260 acquires the time series of the system behavior x (k) (t) from the model behavior storage unit 240, and the target index Evaluation value, i.e., an evaluation function a (k) of system behavior is calculated (S506).
  • the index calculation unit 260 stores the calculated a (k) in the index value storage unit 310.
  • the equivalent differential importance calculation unit 270 obtains H i, j (k) for all j from the random number importance storage unit 250 and calculates the sum for j, thereby obtaining the i th parameter.
  • Equivalent differential importance h i (k) is calculated (S507). Note that the equivalent differential importance calculation unit 270 may perform the process in step S507 only for parameters designated by the user or the like. Note that the process in step S506 and the process in step S507 may be performed in parallel, or may be performed in a reversed order.
  • the process defined by the weight function includes the processes in S503 and S507.
  • the scenario simulator 200 causes the system behavior evaluation function a (k) by the processing shown in FIG. ) And N parameters, the equivalent differential importance h i (k) can be calculated.
  • FIG. 6 is an example of a user interface of the Monte Carlo simulation system in the first embodiment.
  • the user interface 600 includes a simulation condition setting section 610, an expected value calculation section 620, and a sensitivity calculation section 630.
  • the simulation condition setting section 610 includes input reception sections 611 and 612 and a display section 613.
  • the input reception section 611 receives an input of a simulation model, that is, a function x (w, t) that defines an actual value x of system behavior.
  • the input reception section 612 receives an input of the number of simulations by the simulation model.
  • the display section 613 displays the calculation time required for executing the simulation of the contents received by the input receiving sections 611 to 612.
  • the expected value calculation section 620 includes an input reception section 621 and a display section 622.
  • the input receiving section 621 receives an input of an attention index, that is, an evaluation function a (x) of system behavior.
  • the display section 622 displays the calculation accuracy of the system evaluation function A (x).
  • the sensitivity calculation section 630 includes input reception sections 632 and 634, check boxes 631 and 633, and a display section 635.
  • the input reception section 632 receives an input of a parameter z i that is a target for calculating sensitivity ⁇ A / ⁇ z i , that is, a sensitivity calculation target parameter.
  • the check box 631 is a check box for selecting all the parameters z i as target parameters for sensitivity calculation.
  • the input receiving section 634 receives an input of a method for calculating Di, j .
  • the check box 633 is a check box for automatically selecting a method for calculating Di, j .
  • Display section 635 displays the calculation accuracy of the sensitivity ⁇ A / ⁇ z i.
  • the input section and the display section in the user interface 600 may be a display section and an input section, respectively, or may also serve as a display section and an input section, respectively.
  • the display section 622 accepts the input of the calculation accuracy of the system evaluation function A (x)
  • the display section 635 may receive an input of a calculation accuracy of sensitivity ⁇ A / ⁇ z i.
  • the input reception section 611 performs, for example, the number of simulations necessary to guarantee the calculation accuracy of the input system evaluation function A (x) and the input sensitivity ⁇ A / ⁇ z i. May be displayed.
  • the Monte Carlo simulation system of the present embodiment can calculate the sensitivity regarding each parameter of the system evaluation function by using the result of one set of Monte Carlo simulation, thereby reducing the calculation time and the calculation resources. be able to.
  • the Monte Carlo simulation system uses the calculated sensitivity to allow the user to set the parameter value that maximizes the system evaluation function within the allowable range of the parameter value, such as parameter optimization. Search and the like can be performed.
  • the time series x (t) representing the behavior of the system is expressed in the form of Formula 15, that is, when the value y (w, z) determined by the random number vector w and the parameter z is determined, at time t.
  • the system state x (t) is determined deterministically only at times t and y (corresponding to FIG. 3).
  • the probability that the characteristic falls within the desired range, and the like are estimated by Monte Carlo simulation.
  • a circuit designer for example, tunes design parameters and performs a Monte Carlo simulation until the characteristics that take into account the characteristics of the designed circuit satisfy the given specifications, and confirms the characteristics in consideration of the variations. Repeat.
  • the design parameters in each transistor include about 10 related to variations, power, layout area, etc., such as gate width W, gate length L, diffusion layer width LOD, and inter-transistor distance PDX. Therefore, in general, the number of design parameters of the circuit block is about 10 times the number of transistors constituting the circuit. For this reason, design parameters in a complex circuit block composed of a large number of transistors can be enormous, and it is difficult to manually tune individual design parameters as described above.
  • the sensitivity (differential coefficient) for all the many parameters of the expected value of an arbitrary attention index can be obtained by one set of Monte Carlo simulations.
  • the circuit designer can not only determine the characteristics (performance expected value and yield) and the power value of the entire circuit block under the current parameter settings, but also all the transistors. For all parameters, it is possible to simultaneously obtain the sensitivity of how they affect the characteristics (expected performance value and yield) and power of the entire circuit block.
  • the circuit designer can know without changing the parameters of which transistors, the characteristics can be improved without increasing the power consumption, etc. Can greatly improve the efficiency.
  • FIG. 7 shows an example of an SRAM cell circuit composed of six transistors.
  • the SRAM cell needs to satisfy mutually contradicting characteristics that the value stored in the SRAM cell is not lost in normal times and that the stored value of the SRAM cell is properly rewritten to the intended value when the memory is rewritten.
  • the parameter of each transistor included in the SRAM is different from the actual parameter value of each transistor in the manufactured integrated circuit due to variations at the time of manufacturing the integrated circuit. Therefore it is distributed. As a result, for example, a value that is normally stored in the SRAM cell may be lost, or when the memory is rewritten, the stored value of the SRAM cell may not be properly rewritten to the intended value.
  • the circuit designer performs the Monte Carlo simulation by changing the design parameter value to various values, and optimizes the design parameter so that the defective rate of the actually manufactured SRAM cell circuit is equal to or less than a given threshold value.
  • the simulation method of the Monte Carlo simulation system of the first embodiment is applied.
  • the parameter vector z in this embodiment is composed of transistor design parameters. For example, when there are 10 design parameters for each transistor, the parameter vector z in an SRAM cell circuit composed of 6 transistors is a 60-dimensional vector.
  • the characteristic variation of the SRAM circuit is caused by the fact that the actual parameter value of the transistor varies at the time of manufacture. That is, for example, a 60-dimensional parameter y is generated with respect to a 60-dimensional design parameter z that determines the characteristics of the transistor due to variations in manufacturing.
  • the random number generator 210 generates a multidimensional random vector w according to the probability density function g (w) for each simulation.
  • the model behavior update unit 220 can determine the behavior x (t) of the SRAM cell circuit deterministically by normal circuit simulation if the actual parameter y after manufacture is determined. Accordingly, the time series x (t) representing the behavior of the system in the present embodiment is expressed in the form of Equation 15, that is, when a value y (w, z) determined by the random number vector w and the parameter vector z is determined, The system state x (t) at time t is deterministically determined only by times t and y.
  • the evaluation function a (x) of the system behavior is defined as 0 if the behavior of the cell is within the normal operation range, and 1 if it is defective, the system evaluation function A (x), that is, a The expected value E [a (x)] of (x) indicates the defective rate of the SRAM cell.
  • the Monte Carlo simulation system can calculate the sensitivity of the defective rate of the SRAM cell to all 60 design parameters by performing the processing described with reference to FIG.
  • the Monte Carlo simulation system has 60 design parameter values that minimize the steepest descent SRAM cell defect rate from the calculated sensitivity using a known optimization method such as the gradient method or Newton method. Combinations may be calculated.
  • the Monte Carlo simulation system can efficiently obtain the optimum design parameter value in consideration of manufacturing variations in an electronic circuit of a semiconductor integrated circuit including a large number of transistors. .
  • the time series x (t) representing the behavior of the system is expressed in the form of Expression 17, that is, the system state x (t + 1) at time t + 1 is the time t and the past (before time t). ) System behavior, a random number w j newly generated at time t (different for each time), and a parameter z.
  • Financial institutions make plans for the future transaction volume and profits of the financial products they handle when making future business plans, and make management decisions on the appropriateness of the business plan.
  • Fig. 8 shows an example of the balance plan table in a bank.
  • the leftmost column of the income and expenditure schedule shows the items of financial products handled by the bank.
  • the financial products handled include, for example, three categories from A to C for business loans for companies, housing loans for individuals, card loans, and securities dividends held by banks themselves.
  • Actual banks may have several thousand financial product subjects.
  • Each row of the table shows the transaction amount for each item of the financial product in the period shown in the top row. For example, if the current time is the end of 2Q in 2014, the transaction amount for '14 / 1Q and '14 / 2Q will show the actual value, but the transaction value after '14 / 3Q will be the current value (as of the end of 2Q in 2014) Indicates the planned value.
  • the example in FIG. 8 shows a plan for a total of six periods in quarter units. However, in an actual bank, for example, a plan for about five years may be made on a monthly or week basis.
  • the table includes several hundred periods.
  • the relationship between transaction volume and revenue in each subject depends mainly on, for example, bank funding interest rates, that is, spot rates.
  • spot rates The relationship between the spot rate and the profit can be modeled with a certain degree of accuracy for each subject, whereas the future spot rate is uncertain and can only be probabilistically predicted.
  • the expected value and dispersion of earnings are simulated by Monte Carlo simulation to determine the appropriateness of the business plan. Specifically, for example, a random spot is used to generate a future spot rate realization path, calculate the revenue value based on the planned transaction volume of each item, and total the revenues of all subjects. The process of calculating the profit value is performed many times. Through this multiple processing, the expected value and variability of the overall bank profit under the given business plan is estimated.
  • the planned transaction value for each of several hundred future periods that is, the future actual value for each of the total of about 100,000 planned values is the planned value. It is necessary to calculate the impact on the overall profit of the bank when it deviates. If the actual value deviates slightly from the planned value and it is known that the plan will have a significant impact on the overall profits of the bank, reserve the appropriate allowance, or change the business plan itself. Measures such as rebuilding will be necessary.
  • the parameter vector z in the present embodiment is, for example, a planned value of the handling amount in each future period for all subjects. For example, when the number of subjects is about several thousand and the period is several hundred periods, the parameter vector z is a vector of about 100,000 dimensions.
  • the time series x (t) of the behavior of the system is, for example, a combination of a set of profit values of each subject in each future period (for example, a vector of several thousand dimensions) and an interest spot rate. It is a vector of dimensions.
  • the profit value of each subject in the t-th period depends on, for example, the spot rate and the parameter z that is the planned value of the transaction amount of each subject at that time.
  • the future spot rate is expressed by a probabilistic model such as a Full-White model, a CIR model, or a BDT model.
  • the spot rate at time t + 1 is determined based on the spot rate at time t and a random number (newly generated at time t).
  • the behavior x (t) of the system in the present embodiment is expressed in the form of Expression 17 as a whole. That is, the system state x (t + 1) at the time t + 1 includes the time t, the behavior of the system in the past (before the time t), a random number w j newly generated at the time t (different for each time), It is determined by the parameter z.
  • the evaluation function a (x) of the system behavior is, for example, the total value of profits for all subjects in all periods.
  • the system evaluation function A (x) that is, the expected value E [a (x)] of a (x) indicates the expected value of the future profit of the entire bank under the given plan.
  • the Monte Carlo simulation system performs the processing shown in FIG. 5 under the above-described preconditions. It is possible to simultaneously estimate the sensitivity to the expected value of the profit of the entire bank for each planned value of the transaction amount, that is, a total of about 100,000 planned values.
  • the bank can take measures such as, for example, accumulating an appropriate allowance in consideration of the effect of the difference between the planned value and the actual value in advance, or restructuring the business plan itself. Become.
  • the Monte Carlo simulation system of the present embodiment can estimate the sensitivity of the planned value to the profit in the financial field, particularly the bank's profit prediction simulation system. Can be made more efficient.
  • the Monte Carlo simulation system of the present embodiment calculates the boundary residual R i in Equation 9 of the first embodiment.
  • the boundary remainder R i the boundary surface of the distribution range ⁇ random vector w (surface) in ⁇ ( ⁇ w / ⁇ z i)
  • R i 0 in many cases.
  • FIG. 9 shows a two-dimensional random number vector w having two elements, a uniformly distributed random number vector w 1 on the interval [a, b] and a uniformly distributed random number w 2 on the interval [c, d].
  • An example of a concept for calculating the boundary residual R i in the case of a vector is shown.
  • the range indicated by the rectangle in the figure is ⁇ .
  • the probability density function g (w) is expressed by the following equation 19.
  • storage part 300 in a present Example further contains the sample path fixed differential memory
  • the sample path fixed differential storage unit 330 holds D i, j calculated by the sample path fixed differential calculation unit 231.
  • the post-processing calculation unit 400 in the present embodiment further includes a boundary residual calculation unit 430 (not shown) that is a program.
  • the boundary residual calculation unit 430 calculates the boundary residual R i using the random number vector w, the parameter vector z, and D i, j received from the sample path fixed differential storage unit 330 as input, and calculates the calculated boundary residual R i. Is transmitted to the expected value calculation unit 410.
  • Equation 19 The first method by which the boundary residual calculation unit 430 calculates the boundary residual R i is based on the first line of Equation 9. In the condition of FIG. 9, the first line of Equation 9 is expanded as shown in Equation 20 below.
  • the boundary residual calculation unit 430 fixes w 1 to b which is a value on the boundary of ⁇ , moves w 2 , and calculates ( ⁇ w 2 / ⁇ z i )
  • x const respectively. Then, after w 1 is fixed to the boundary value a, w 2 is moved to calculate a value obtained by subtracting the integral value obtained by calculating ( ⁇ w 2 / ⁇ z i )
  • x const . Further, the boundary residual calculation unit 430 moves the w 1 after fixing w 2 to the value d on the boundary of the distribution, and calculates the integral value obtained by calculating ( ⁇ w 1 / ⁇ z i )
  • x const , respectively. from each move w 1 in terms of fixing the w 2 to the boundary value c ( ⁇ w 1 / ⁇ z i)
  • and minus the integral values calculated for x const, a value is calculated by adding the .
  • the boundary residual calculation unit 430 calculates the boundary residual R i by adding the two calculated values.
  • the boundary residual calculation unit 430 may calculate each of the integration values described above using, for example, a known numerical integration method.
  • the second method in which the boundary residual calculation unit 430 calculates the boundary residual R i is based on the third line of Equation 9. Since the third line of Equation 9 is integration over the entire distribution range ⁇ of the random number vector, the boundary residual calculation unit 430 calculates the boundary residual R i as an expected value using Monte Carlo simulation from the following Equation 21. be able to.
  • the random number generator 210 generates a large number of random numbers according to the given probability density function g (w).
  • Boundary remainder calculation unit 430 for each generated random vectors, the value of the second row of square brackets in Equation 21 is calculated, the average value of the calculated value may be set as the estimated value of R i.
  • the boundary residual calculation unit 430 may calculate div in Equation 21 by numerical differentiation. When div is given as a function by a user or the like, the boundary residual calculating unit 430 may calculate div by substituting the parameter vector z and the random number vector w into the function. Note that the boundary residual calculation unit 430 can similarly calculate the boundary residual R i when the probability density function g (w) is other than the two-dimensional uniform distribution.
  • the sensitivity calculation unit 420 calculates a difference obtained by subtracting the boundary residual R i calculated by the above method from the calculated average value. The calculated difference is an estimated value of the sensitivity.
  • the Monte Carlo simulation system of the present embodiment can calculate the value when the boundary residual R i is not zero.
  • the Monte Carlo simulation system calculates the boundary residual R i, thereby calculating x, whatever the distribution of the random number vector applied to the model for calculating the actual value x of the system behavior.
  • the sensitivity calculation method of this embodiment can be applied to all models.
  • this invention is not limited to the above-mentioned Example, Various modifications are included.
  • the above-described embodiments have been described in detail for easy understanding of the present invention, and are not necessarily limited to those having all the configurations described.
  • a part of the configuration of a certain embodiment can be replaced with the configuration of another embodiment, and the configuration of another embodiment can be added to the configuration of a certain embodiment.
  • each of the above-described configurations, functions, processing units, processing means, and the like may be realized by hardware by designing a part or all of them with, for example, an integrated circuit.
  • Each of the above-described configurations, functions, and the like may be realized by software by interpreting and executing a program that realizes each function by the processor.
  • Information such as programs, tables, and files that realize each function can be stored in a memory, a hard disk, a recording device such as an SSD (Solid State Drive), or a recording medium such as an IC card, an SD card, or a DVD.
  • control lines and information lines indicate what is considered necessary for the explanation, and not all the control lines and information lines on the product are necessarily shown. Actually, it may be considered that almost all the components are connected to each other.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Operations Research (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

　シミュレーションシステムは、複数の乱数ベクトルと、パラメータベクトルと、を取得し、取得した複数の乱数ベクトルそれぞれに対応する確率的システムの挙動の実現値を算出し、取得した複数の乱数ベクトルそれぞれと、取得したパラメータベクトルと、重み関数と、に基づき、取得したパラメータベクトルのパラメータそれぞれに対する、取得した乱数ベクトルそれぞれの重みを算出し、取得した複数の乱数ベクトルそれぞれに対応する確率的システムの挙動の評価値を算出し、取得した複数の乱数ベクトルそれぞれに対応する算出した挙動の評価値と、取得したパラメータベクトルから選択したパラメータに対する取得した乱数ベクトルそれぞれの重みと、に基づいて、選択したパラメータに対する、確率的システムの挙動の評価値の期待値の感度、を算出する。

Description

シミュレーションシステム、及びシミュレーション方法
　本発明は、シミュレーションシステム、及びシミュレーション方法に関する。
　興味の対象であるシステム（系）が与えられたとき、当該システムの振る舞いをモデル化し、計算機が当該モデルに基づいて当該システムの振る舞いをシミュレーションし、将来の振る舞いを予測することが行われる。例えば、ある種の物理現象や、多数の要素が絡み合う社会現象などは、単純な決定論的なモデルでは記述されることができず、複雑かつ確率的なモデルとして記述される。
　決定論的なシステムは、与えられた初期条件のもとで、システムの将来の振る舞いが一つに定まる。これに対して、確率的なシステムは、システムの将来の振る舞いはある確率分布にしたがって分布しており、計算機は、例えば、注目する指標の期待値等の統計量を求める場合、シミュレーションを行う必要がある。
　確率的なシステムの統計的な振る舞いを計算機が予測、シミュレーションする方法として、モンテカルロ・シミュレーション法が知られている。モンテカルロ・シミュレーション法は、対象システムに存在する確率的な振る舞いを乱数によって模擬したシミュレーションを多数回繰り返し行うことで、注目する指標の期待値を推定する方法である。このとき、シミュレーションの回数を増やすことで、大数の法則により、注目する指標の期待値が真の値に近づくと考えられる。
　本技術分野の背景技術として、特開２００４－１１７２２８号公報（特許文献１）がある。特許文献１には、「気象物理量を対象とした気象庁発表の長期予報を取得するデータ取得ステップ０２０１、気象物理量の従う時系列モデルを作成する気象時系列モデル作成ステップ０２０２、推定期間における推定地点の気象物理量をシミュレーション回数だけシミュレートする気象物理量シミュレーションステップ０２０３、気象物理量のシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力ステップ０２０４を有する。さらに、気象時系列モデル作成ステップ０２０２を、気象時系列モデルのパラメータ推定ステップ０２０２１および気象時系列モデルのパラメータ補正ステップ０２０２２から構成する。」と記載されている（要約参照）。
特開２００４－１１７２２８号公報
　例えば、特許文献１に記載の技術のシミュレーション対象である気象物理量のように、モンテカルロ・シミュレーションの対象である、確率的なモデルには、一般的に多数のパラメータが存在することが多い。このとき、モンテカルロ・シミュレーション実行にあたっては、あるパラメータ値の下での注目する指標の期待値だけではなく、注目指標の期待値のパラメータ値に対する感度を算出したい場合がある。注目指標の期待値のパラメータ値に対する感度とは、与えられたパラメータ値が微小に変化したときに、注目指標の期待値がどれくらい変化するかを表す微分係数である。
　例えば、システムのパラメータの最適化を行う、すなわち、許容されるパラメータ値の範囲の中で、ある注目指標の期待値を最大化するパラメータ値を探す場合、与えられたパラメータ値の下での感度を得られれば、勾配法やニュートン法といった公知の最適化手法によって最適化を行うことができる。
　ところが、モンテカルロ・シミュレーションにおいて、注目指標の期待値のパラメータ値に対する感度を求めることは、とくにパラメータが多数ある場合には、容易ではない。注目指標の期待値のパラメータに対する感度は、単純には、与えられたパラメータ値の下での注目指標の期待値と、パラメータをそこから微小に変化させたときの注目指標の期待値との差を、パラメータの微小変化量で割ることで算出できる。
　特許文献１に記載の技術を、多数のパラメータについて注目指標の期待値の感度の推定に適用する場合、例えば特許文献１の数５、数１２のように、多数のパラメータそれぞれを一つずつ微小変化させた注目指標の期待値それぞれを算出する必要がある。すなわち、特許文献１に記載の技術を適用した感度推定手法は、当該多数パラメータの個数と同数セットのモンテカルロ・シミュレーションを行う必要がある。
　そもそも、１セットのモンテカルロ・シミュレーションを行って注目指標の期待値を求めること自体が、大数の法則を満足するような多数回のシミュレーションをする必要がある。したがって、当該多数パラメータの個数と同数セットのモンテカルロ・シミュレーションを行うことは、計算機リソースや計算時間の制限により困難である。
　上記課題を解決するために、本発明は、例えば、以下のような構成を採用する。複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムであって、プロセッサと、メモリと、を含み、前記メモリは、１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、前記プロセッサは、複数の乱数ベクトルと、パラメータベクトルと、を取得し、前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、シミュレーションシステム。
　本発明の一態様によれば、パラメータが存在する確率的なシステムにおいて、任意に決めた注目指標の期待値のパラメータに対する感度を、１セットのモンテカルロ・シミュレーションによって得られた結果から求めることができる。これにより、注目指標の期待値のパラメータに対する感度を求めるために必要な、計算機リソースや計算時間を削減することができる。
実施例１におけるモンテカルロ・シミュレーション・システムの構成例を表すブロック図である。 実施例１における乱数重要度算出部の構成例を表すブロック図である。 実施例１のモンテカルロ・シミュレーションの全体処理の例を表すフローチャートである。 実施例１におけるシミュレーション実行処理の一例を表すフローチャートである。 実施例１におけるシミュレーション実行処理の別の一例を表すフローチャートである。 実施例１におけるモンテカルロ・シミュレーション・システムのユーザインタフェースの一例である。 実施例２における６個のトランジスタで構成されるＳＲＡＭセルの回路の例である。 実施例３における、銀行の収支計画表の例である。 実施例４において、乱数ベクトルが２次元一様分布に従う場合における、境界残余を算出する概念の例を示す図である。
　以下、本発明の実施形態について述べる。本実施形態は本発明を実現するための一例に過ぎず、本発明の技術的範囲を限定するものではないことに注意すべきである。本実施形態においては、便宜上その必要があるときは、複数のセクションまたは実施形態に分割して説明するが、特に明示した場合を除き、それらは互いに無関係なものではなく、一方は他方の一部または全部の変形例、詳細、補足説明等の関係にある。また、以下の実施形態において、要素の数等（個数、数値、量、範囲等を含む）に言及する場合、特に明示した場合および原理的に明らかに特定の数に限定される場合等を除き、その特定の数に限定されるものではなく、特定の数以上でも以下でもよい。
　以下、本発明の実施形態を、添付図面を参照して説明する。なお、各図において、共通の構成には原則として同一の参照符号を付し、繰り返しの説明は省略する。
　本実施例はモンテカルロ・シミュレーション・システムについて説明する。まず、本実施例のモンテカルロ・シミュレーション・システムがシミュレーションを実行し、感度を算出するための理論的側面について説明する。
　確率変数Ｘは、本実施例のモンテカルロ・シミュレーション・システムによる、シミュレーション対象である確率的なシステムの挙動を表す。つまり、当該シミュレーション対象である確率的なシステムの挙動は確率変数Ｘでモデル化される。以下、シミュレーション対象である確率的なシステムを、単にシステムとも表記する。また、ｘはシステムの挙動Ｘの実現値を表す。
　ここで、システムの挙動Ｘは、１次元又は多次元のベクトルのいずれであってもよく、また、「時刻」と呼ばれる変数ｔで区別される無限個の確率変数の組（離散時間確率過程と呼ばれる）であってもよい。したがって、システムの挙動Ｘが離散時間確率過程である場合、Ｘは本来Ｘ（ｔ）と表記されるべきであるが、特に必要がある場合を除いて、単にＸとも表記する。また、ある時刻ｔ１におけるシステムの挙動を示す値ｘ（ｔ１）を、時刻ｔ１におけるシステムの状態と呼ぶ。
　システムの挙動Ｘの分布は、システムのパラメータｚに依存し、さらにＸの確率密度関数が、ｆ（ｘ，ｚ）という形で表記されるとする。以下、システムのパラメータｚを単にパラメータｚと表記する。パラメータｚは、Ｎ次元のベクトル（すなわち、システムの挙動はＮ個のパラメータに依存している）であり、ｚの成分は、ｚ＝（ｚｉ）（ただし、ｉ＝１，…，Ｎ）であるとする。また、１つのシステムの挙動の実現値Ｘ＝ｘに対して１つの実数値を定める、関数ａ（ｘ）が与えられているとする。ａ（ｘ）をシステムの挙動の評価関数と呼ぶ。システムの挙動の評価関数は、注目指標の一例である。さらに、ａ（ｘ）の期待値Ａ（ｘ）＝Ｅ［ａ（ｘ）］をシステム評価関数と呼ぶ。
　実際に興味あるシステムについて検討する場合、しばしば、Ｘの確率密度関数ｆ（ｘ，ｚ）が未知であることが問題となる。しかしながら、多くの場合、Ｘの確率密度関数ｆ（ｘ，ｚ）自体は未知であっても、Ｘの分布を生み出す挙動自体は既知（又はモデル化可能）である。
　この場合、モンテカルロ・シミュレーション・システムは、モンテカルロ・シミュレーションによって、システムの挙動の実現値ｘを多数生成して、Ｘの分布、すなわちｆ（ｘ，ｚ）を模擬することで、システム評価関数Ａ（ｘ）の値を近似的に算出できる。なお、生成されたシステムの挙動の実現値ｘそれぞれをパスとも呼ぶ。モンテカルロ・シミュレーションは、分布が既知である（複数の）乱数から、システムの挙動の実現値ｘを算出する計算を含む。以下、これを定式化する。
　Ｗ＝（Ｗｊ）をＭ次元の確率変数とし（ただし、ｊ＝１，…，Ｍ）、確率変数Ｗの実現値をｗ＝（ｗｊ）と書く。また、Ｗの（同時）確率密度関数をｇ（ｗ）と書く。ここで、ｇ（ｗ）は既知であり、また計算機上でこの分布に従う乱数を容易に生成できるものとする。ｗは、例えば、［０，１］上の一様分布、又は、標準正規分布Ｎ（０，１）等に従うＭ個の独立な乱数からなるベクトルである。このとき、パラメータｚ＝（ｚｊ）のもとで、乱数ベクトルｗから、対応するＸの実現値ｘを計算する方法、すなわち、関数ｘ＝ｘ（ｗ,ｚ）が与えられているとする。このとき、以下の数１が成立する。
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
　数１８は、パラメータの微小変化Δｚｉを打ち消す乱数の微小変化Δｗｊが、初期時刻ｔ＝０からシミュレーションの終了時刻ｔ＝Ｔまで順に定まることを示す。
　以下、本実施例のモンテカルロ・シミュレーション・システムがシミュレーションを実行し、感度を算出する実践的側面について説明する。図１は、本実施例のモンテカルロ・シミュレーション・システム１００の構成例を示す。モンテカルロ・シミュレーション・システム１００は、例えば、ＣＰＵ１１０、メモリ１２０、入出力インタフェース１３０、二次記憶装置１４０を含む計算機上に構成される。また、図１において、モンテカルロ・シミュレーション・システムは１つの計算機から構成されているが、複数の計算機から構成されてもよい。
　ＣＰＵ１１０は、プログラムに従って動作するプロセッサ及び／又は論理回路を含み、データの入力／出力、読み込み／書き込みを行い、さらに、後述する各プログラムを実行する。メモリ１２０は、ＣＰＵ１１０が実行するプログラム及びデータを一時的にロードして記憶し、さらに各プログラム及び各データを保持する。入出力インタフェース１３０は、外部装置等からデータ等の入力、及び外部装置等に対してデータ等の出力を行うインタフェースである。二次記憶装置１４０は、不揮発性の記憶装置であり、例えばＨＤＤ等である。プログラムやデータ等は二次記憶装置１４０に格納されてもよい。
　プログラムはＣＰＵ１１０によって実行されることで、定められた処理をメモリ１２０及び入出力インタフェース１３０を用いながら行う。従って、本実施例及び他の実施例においてプログラムを主語とする説明は、ＣＰＵ１１０を主語とした説明でもよい。若しくは、プログラムが実行する処理は、そのプログラムが動作する計算機及び計算機システムが行う処理である。
　ＣＰＵ１１０は、プログラムに従って動作することによって、所定の機能を実現する機能部として動作する。例えば、ＣＰＵ１１０は、後述する乱数発生部２１０に従って動作することで乱数発生部として機能し、後述するモデル挙動更新部２２０に従って動作することでモデル挙動更新部として機能する。さらに、ＣＰＵ１１０は、各プログラムが実行する複数の処理のそれぞれを実現する機能部としても動作する。計算機及び計算機システムは、これらの機能部を含む装置及びシステムである。プログラムは二次記憶装置１４０に格納されていてもよい。
　メモリ１２０は、１又は複数のシナリオ・シミュレータ２００と、記憶部３００と、後処理計算部４００と、モンテカルロ・シミュレーション全体制御部５００と、を含む。シナリオ・シミュレータ２００は、異なる乱数ベクトルを用いてＬ回のモンテカルロ・シミュレーションを行う。記憶部３００は、シナリオ・シミュレータ２００から出力された数値等を記憶する。後処理計算部４００は、シミュレーション終了後に期待値および感度の計算を行う。
　モンテカルロ・シミュレーション全体制御部５００は、モンテカルロ・シミュレーション・システム１００全体の制御を行う。また、モンテカルロ・シミュレーション全体制御部５００は、ユーザ等が指定したシミュレーション回数によって算出されるシステム評価関数及び感度の推定値の計算精度を算出し、出力してもよい。また、モンテカルロ・シミュレーション全体制御部５００は、ユーザ等が指定した、システム評価関数及び感度の推定値の計算精度を保証するために必要なシミュレーション回数等を算出、出力してもよい。モンテカルロ・シミュレーション全体制御部５００は、例えば、指定された又は算出したシミュレーション回数、例えばＬ回、シナリオ・シミュレータ２００を起動し、Ｌ回のシミュレーションが終了すると、後処理計算部４００を起動する。
　シナリオ・シミュレータ２００は、ユーザ等から与えられるシミュレーションモデルの１次元又は多次元のパラメータのベクトルｚを入力として、シミュレーションに必要な乱数ベクトルを内部で生成しながらモデルのシミュレーションを実行する。また、シナリオ・シミュレータ２００は、システムの挙動の評価関数ａ（ｋ）と、感度計算に必要な等価微分重要度ｈｉ （ｋ）を出力する。
　シナリオ・シミュレータ２００は、それぞれプログラムである、乱数発生部２１０、モデル挙動更新部２２０、乱数重要度算出部２３０、指標算出部２６０、等価微分重要度算出部２７０、及びシミュレーション制御部２８０を含む。また、シナリオ・シミュレータ２００は、それぞれデータを格納する領域である、モデル挙動記憶部２４０、及び乱数重要度記憶部２５０を含む。
　乱数発生部２１０は、予め定められた、又はユーザ等によって指定されたシミュレーションモデルに対応する確率密度関数に従う乱数ベクトルを生成する。モデル挙動更新部２２０は、ユーザ等によって与えられたパラメータと乱数発生部２１０が生成した乱数ベクトルとを入力として、入力されたパラメータと乱数ベクトルとを予め定められた、又はユーザ等によって指定されたシミュレーションモデル（関数）に適用し、システムの挙動の実現値をシミュレートする。
　乱数重要度算出部２３０は、パラメータと乱数発生部２１０が生成した乱数ベクトルとを入力として、入力された各パラメータについて、ｗｊ （ｋ）の各パラメータに関する乱数重要度、すなわち数８の総和記号の中の項を算出する。モデル挙動記憶部２４０は、モデル挙動更新部２２０が算出したシステムの挙動の実現値を格納する。乱数重要度記憶部２５０は、乱数重要度算出部２３０が算出した乱数重要度を格納する。指標算出部２６０は、例えば、モデル挙動記憶部２４０に格納されたシステムの挙動の実現値の平均を算出することにより、システムの挙動の評価関数を算出する。
　等価微分重要度算出部２７０は、例えば、乱数重要度記憶部２５０に格納された乱数重要度の和を算出することにより等価微分重要度を算出する。シミュレーション制御部２８０は、シナリオ・シミュレータ２００全体の動作を制御する。
　記憶部３００は、指標算出部２６０が算出したシステムの挙動の評価関数ａ（ｋ）を格納する指標値記憶部３１０と、等価微分重要度算出部２７０が算出した等価微分重要度を格納する等価微分重要度記憶部３２０と、を含む。
　後処理計算部４００は、それぞれプログラムである、期待値計算部４１０、及び感度計算部４２０を含む。期待値計算部４１０は、指標値記憶部３１０に格納されたＬ個のシステムの挙動の評価関数ａ（ｋ）の平均値を計算し、当該平均値を期待値、すなわちシステム評価関数Ａ（ｚ）として出力する。このとき、期待値計算部４１０は、さらにＬ個のシステムの挙動の評価関数ａ（ｋ）の分散を計算し、Ａ（ｚ）の推定値の精度（信頼区間）を同時に出力してもよい。
　感度計算部４２０は、指標値記憶部３１０に格納されたシステムの挙動の評価関数ａ（ｋ）と、等価微分重要度記憶部３２０に格納されたＮ個のパラメータそれぞれに関する等価微分重要度ｈｉ （ｋ）の組について、それぞれ対応する値同士の積を計算する。さらに、感度計算部４２０は、算出した積の平均値を計算して、Ｎ個のパラメータそれぞれに関するシステム評価関数Ａ（ｚ）の感度∂Ａ／∂ｚｉとして出力する。このとき、感度計算部４２０は、さらにａ（ｋ）とｈｉ （ｋ）の積の分散を計算し、感度∂Ａ／∂ｚｉの推定値の精度（信頼区間）を同時に出力してもよい。
　図２は、乱数重要度算出部２３０の構成例を示す。乱数重要度算出部２３０は、いずれもプログラムである、サンプルパス固定微分計算部２３１、サンプルパス固定２階微分計算部２３２、及び乱数重要度計算部２３３を含む。
　サンプルパス固定微分計算部２３１は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出する。サンプルパス固定２階微分計算部２３２は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として∂Ｄｉ，ｊ／∂ｗｊを算出する（Ｄｉ，ｊ （ｋ）＝（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ）。サンプルパス固定２階微分計算部２３１は、サンプルパス固定微分計算部２３１が算出した（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔと、パラメータと、を入力として、∂Ｄｉ，ｊ／∂ｗｊを算出してもよい。また、サンプルパス固定２階微分計算部２３２が（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出してもよく、このとき、サンプルパス固定微分計算部２３１は（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出しなくてもよい。
　乱数重要度計算部２３３は、乱数発生部２１０から入力された乱数ベクトルと、サンプルパス固定微分計算部２３１又はサンプルパス固定２階微分計算部２３２が算出した（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔと、サンプルパス固定２階微分計算部２３２が算出した∂Ｄｉ，ｊ／∂ｗｊと、を入力として、各パラメータについて、ｗｊ （ｋ）の各パラメータに関する乱数重要度、すなわち数８の総和記号の中の項を算出する。
　図３は、本実施例におけるモンテカルロ・シミュレーションの全体処理の一例を示す。モンテカルロ・シミュレーションの全体処理は、シナリオ・シミュレータ２００が異なる乱数ベクトルを用いてＬ回のシミュレーションを行うシミュレーション実行処理と、シミュレーション終了後に後処理計算部４００が行う期待値・感度計算処理と、を含む。
　シナリオ・シミュレータ２００は、シミュレーション実行処理を行う（Ｓ３０１）。シミュレーション実行処理は、１セットのモンテカルロ・シミュレーションを構成するＬ回のシミュレーションを実行する処理である。シナリオ・シミュレータ２００は、異なる乱数ベクトルを用いて個々のシミュレーションを実行する。
　なお、本実施例のモンテカルロ・シミュレーション・システムは、ｋ回目のシミュレーションにより得られたシステムの挙動の実現値ｘ（ｋ）から、システムの挙動の評価関数ａ（ｋ）に加えて、Ｎ個のパラメータそれぞれについて等価微分重要度ｈｉ （ｋ）を計算する。なお、シミュレーション実行処理におけるＬ回のシミュレーションそれぞれは、互いに独立であるため、Ｌ回のシミュレーションを逐次的に行うのではなく、例えば、多数のシナリオ・シミュレータ２００又は多数の計算機を用いて同時に並列実行してもよい。等価微分重要度ｈｉ （ｋ）を算出する処理を定める一連の関数を、重み関数と呼ぶ。
　続いて、後処理計算部４００は、期待値・感度計算処理を行う（Ｓ３０２）。期待値・感度計算処理は、１セットのモンテカルロ・シミュレーションを構成するＬ回のシミュレーションの実行後に、後処理計算部４００が、注目指標の期待値の推定値と、注目指標の期待値のパラメータに対する感度の推定値と、を計算する処理である。
　注目指標の期待値、すなわちシステム評価関数Ａ（ｚ）の推定値は、期待値計算部４１０が、指標値記憶部に格納されたＬ個のシステムの挙動の評価関数ａ（ｋ）の平均を計算することで得られる。本実施例におけるモンテカルロ・シミュレーション・システムは、これに加えて、注目指標の期待値のＮ個のパラメータそれぞれに対する感度の推定値を計算する。
　具体的には、感度計算部４２０は、システム評価関数Ａ（ｚ）のｉ番目のパラメータそれぞれに対する感度の推定値を、指標値記憶部３１０に格納された対応するシステムの挙動の評価関数ａ（ｋ）に、等価微分重要度記憶部３２０に格納されたｉ番目のパラメータについての等価微分重要度ｈｉ （ｋ）を掛けることで重みづけをしたものの平均を計算することで得る。なお、モンテカルロ・シミュレーション・システムは、例えば、ユーザ等が指定したパラメータについてのみ、ステップＳ３０１における等価微分重要度の算出処理、及びステップＳ３０２における感度の推定値の算出処理を行ってもよい。
　なお、前述のように、例えば、システムの挙動を表す時系列ｘ（ｔ）が、数１５又は数１７の形で表される場合、本実施例における注目指標の期待値のパラメータに対する感度の推定手法が適用できる。図３、図４は、それぞれ、システムの挙動を表す時系列ｘ（ｔ）が、数１５、数１７の形で表されたモデルで算出される場合における、ｋ回目のシミュレーション実行処理の詳細な処理内容を示す。
　図４は、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、すなわち、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）を定めると、時刻ｔでのシステムの状態ｘ（ｔ）が，時刻ｔとｙのみで決定論的に定まる場合について、シミュレーション実行処理の詳細な処理内容の一例を示す。
　まず、シナリオ・シミュレータ２００は、乱数ベクトルｗ（ｋ）を生成し、ｗ（ｋ）とパラメータｚとからｙ（ｋ）を計算する（Ｓ４０１）。具体的には、乱数発生部２１０は、予め定められた確率密度関数ｇ（ｗ）にしたがうＭ次元の乱数ベクトルｗ（ｋ）を生成し、生成した乱数をモデル挙動更新部２２０及び乱数重要度算出部２３０に送信する。なお、乱数発生部２１０は、例えば、既知の疑似乱数発生アルゴリズム、又は物理乱数発生装置といった公知の手段を用いて乱数ベクトルを生成すればよい。モデル挙動更新部２２０は、受信したｗ（ｋ）と、与えられたＮ次元のパラメータ値ｚを、予め与えられた関係式ｙ（ｗ，ｚ）に代入してｙ（ｋ）を計算する。
　続いて、シナリオ・シミュレータ２００は、Ｎ個のパラメータそれぞれについて等価微分重要度ｈｉ （ｋ）を計算する（Ｓ４０２）。これは、注目指標の期待値のパラメータに対する感度の推定のために必要な処理であり、重み関数が定める処理は当該処理を含む。なお、シナリオ・シミュレータ２００は、ユーザ等が指定したパラメータについてのみ、ステップＳ４０２における処理を行ってもよい。ｉ番目のパラメータについての等価微分重要度ｈｉ （ｋ）は、具体的には以下の４つのステップによって算出される、すなわちステップＳ４０２は以下の４つのステップを含む。
　第１のステップは、サンプルパス固定微分計算部２３１が、システムの挙動ｘを一定に保つような、乱数ベクトルｗ（ｋ）のｊ番目の要素とｉ番目のパラメータの微小変化の関係Ｄｉ，ｊ （ｋ）＝（∂ｗｊ （ｋ）／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出するステップである。ここで、添え字ｊは、Ｍ次元の乱数ベクトルｗ（ｋ）のｊ番目の要素についての計算を表す。
　システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、数１６が成り立つから、サンプルパス固定微分計算部２３１は、（∂ｗｊ／∂ｚｉ）｜（ｙ＾（ｋ）＝ｃｏｎｓｔ）を求めればよい。つまり、サンプルパス固定微分計算部２３１は、例えば、関数－（∂ｙ（ｋ）／∂ｚｉ）／（∂ｙ（ｋ）／∂ｗｊ （ｋ））がユーザ等によって与えられた場合、入力されたパラメータｚｉ及び入力された乱数ベクトルｗｊ （ｋ）を当該関数に代入することにより、Ｄｉ，ｊ （ｋ）を算出すればよい。サンプルパス固定微分計算部２３１は、当該関数から高速にＤｉ，ｊ （ｋ）を算出することができる。
　また、サンプルパス固定微分計算部２３１は、例えば、パラメータの微小変化Δｚｉに対して、ｙを一定に保つようなΔｗを数値的に探索してΔｗ／Δｚｉを計算することにより、（∂ｗｊ／∂ｚｉ）｜（ｙ＾（ｋ）＝ｃｏｎｓｔ）、すなわちＤｉ，ｊ （ｋ）を近似的に算出してもよい。例えば、ユーザ等が関数ｘ及びｙの具体的な形がわからない場合であっても、サンプルパス固定微分計算部２３１は、上述の方法により、Ｄｉ，ｊ （ｋ）を算出することができる。サンプルパス固定微分計算部２３１は算出したＤｉ，ｊ （ｋ）を、サンプルパス固定２階微分計算部２３２及び乱数重要度計算部２３３に送信する。
　第２のステップは、サンプルパス固定２階微分計算部２３２が、受信したＤｉ，ｊ （ｋ）を、ｗｊ （ｋ）で微分するステップである。Ｄｉ，ｊ （ｋ）を、ｗｊ （ｋ）で微分した結果をＤＤｉ，ｊ （ｋ）と表記する。サンプルパス固定２階微分計算部２３２は、例えば、受信したＤｉ，ｊ （ｋ）を数値的に微分することによって、ＤＤｉ，ｊ （ｋ）を算出すればよい。また、サンプルパス固定２階微分計算部２３２は、例えば、ＤＤｉ，ｊ （ｋ）を算出する関数が与えられている場合、入力されたパラメータｚｉ及び入力された乱数ベクトルｗｊ （ｋ）当該関数に代入することによって、ＤＤｉ，ｊ （ｋ）を算出してもよい。サンプルパス固定２階微分計算部２３２は、算出したＤＤｉ，ｊ （ｋ）を、乱数重要度計算部２３３に送信する。
　第３のステップは、乱数重要度計算部２３３が、乱数ベクトルｗ（ｋ）のｊ番目の要素ｗｊ （ｋ）に関する等価微分重要度Ｈｉ，ｊ （ｋ）を計算するステップである。第３のステップは、数８における総和記号の中の項を計算することに相当する。具体的には、例えば、乱数重要度計算部２３３は、受信したＤｉ，ｊ （ｋ）と、受信したＤＤｉ，ｊ （ｋ）と、予め定められた確率密度関数ｇ（ｗ）から一意に定まる関数ＬＧに乱数ベクトルｗｊ （ｋ）を代入した値ＬＧj（ｗｊ （ｋ））＝（１／ｇ（ｗｊ （ｋ）））×（∂ｇ／∂ｗｊ （ｋ））と、を図中に示した関係式に代入して乱数重要度Ｈｉ，ｊ （ｋ）を計算すればよい。乱数重要度計算部２３３は、算出したＨｉ，ｊ （ｋ）を乱数重要度記憶部２５０に格納する。
　なお、前述のように、第３のステップにおける計算は、Ｍ次元のｗが一様分布の乱数ベクトルである場合、又は、ｗがＭ次元独立標準正規分布の乱数ベクトルである場合には、それぞれ数１１又は数１２のように表される。従って、例えば、乱数ｗｊ （ｋ）が一様分布に従う場合、数１１より、乱数重要度はＤＤｉ，ｊ （ｋ）で表されるため、乱数重要度算出部２３０は、第１のステップのみを行うことにより乱数重要度を算出することができる。シナリオ・シミュレータ２００は、第１のステップから第３のステップにおける処理を、全てのｊ（ただし、ｊ＝１…Ｍ）について行う。
　第４のステップは、等価微分重要度算出部２７０がｉ番目のパラメータに関する等価微分重要度ｈｉ （ｋ）を計算するステップである。第４のステップは、前述の数８における総和計算に相当する。等価微分重要度算出部２７０は、乱数重要度記憶部２５０に格納された乱数重要度Ｈｉ，ｊ （ｋ）を取得し、ｊについての総和を計算する。なお、各パラメータｚｉについての等価微分重要度ｈｉ （ｋ）の計算は互いに独立であるため、１つのシナリオ・シミュレータ２００によってＮ個のパラメータｚ１，…，ｚNについて等価微分重要度の計算を逐次的に行うのではなく、例えば、多数の計算機を用いて同時に並列に当該計算を実行してもよい。
　続いて、モデル挙動更新部２２０は、は、実際にシステムの挙動ｘ（ｋ）（ｔ）をシミュレートし計算する（Ｓ４０３）。モデル挙動更新部２２０は、時刻ｔを初期時刻から終了時刻に向かって増加させながら、時刻ｔでのシステムの状態を逐次的に計算することで当該処理を行うことができるが、これによらずその他の公知の効率的なシミュレーションの計算手法を用いてもよい。また、モデル挙動更新部２２０は算出したｘ（ｋ）（ｔ）をモデル挙動記憶部２４０に格納する。なお、ステップＳ４０３における処理は、例えば、ステップＳ４０２における処理と並行して行われてもよいし、ステップＳ４０２における処理の前に行われてもよい。
　続いて、指標算出部２６０は、モデル挙動記憶部２４０から、システムの挙動ｘ（ｋ）（ｔ）を取得し、注目指標の評価値、すなわち、システムの挙動の評価関数ａ（ｋ）を計算する（Ｓ４０４）。指標算出部２６０は、算出したａ（ｋ）を指標値記憶部３１０に格納する。
　以上のように、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、図４に示した処理により、シナリオ・シミュレータ２００はシステム挙動の評価関数ａ（ｋ）、およびＮ個のパラメータそれぞれについての等価微分重要度ｈｉ （ｋ）を計算することができる。
　図５は、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される場合、すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの状態と、時刻ｔで新たに生成する（時刻毎に異なる）乱数ｗｊ（複数個の乱数を含む乱数ベクトルもよい）と、パラメータｚと、によって定まる場合の、シミュレーション実行処理の詳細な処理内容を示す。
　シナリオ・シミュレータ２００は、時刻ｔを、シミュレーション開始時刻ｔ＿ｓｔａｒｔから、シミュレーション終了時刻ｔ＿ｅｎｄまで、シミュレーション時間ステップΔｔずつ増加させながら、時刻ｔにおけるシステムの状態ｘ（ｋ）（ｔ）と、Ｎ個のパラメータ全てについてｊ番目の乱数ｗｊ （ｋ）に関する乱数重要度Ｈｉ，ｊ （ｋ）と、を計算する処理を繰り返し行うことでシミュレーションを行う。さらに、シナリオ・シミュレータ２００は、シミュレーションの結果得られたシステムの挙動ｘ（ｋ）（ｔ）の時系列を用いて、システムの挙動の評価関数ａ（ｋ）、およびＮ個のパラメータそれぞれについての等価微分重要度ｈｉ （ｋ）を計算する。
　まず、シミュレーション制御部２８０は、時刻ｔをシミュレーション開始時刻ｔ＿ｓｔａｒｔにセットし、時刻ｔに関するループ処理の開始をシナリオ・シミュレータ２００の各部に通知する（Ｓ５０１）。時刻ｔに関するループ処理は、モデル挙動更新部２２０が時刻ｔでのシステムの状態ｘ（ｋ）（ｔ）を計算する処理（Ｓ５０２）と、乱数重要度算出部２３０が時刻ｔで生成した乱数のＮ個のパラメータそれぞれに関して乱数重要度Ｈｉ，ｊ （ｋ）を計算する処理（Ｓ５０３）と、を含む。
　ステップＳ５０２における具体的な処理について説明する。まず、乱数発生部２１０は、時刻ｔでのシステムの状態ｘ（ｋ）（ｔ）の計算に必要な乱数ｗｊ （ｋ）を、予め定められた確率密度関数ｇ（ｗ）にしたがって生成し、モデル挙動更新部２２０、及び乱数重要度算出部２３０に生成した乱数を送信する。モデル挙動更新部２２０は、受信した乱数ｗｊ （ｋ）と過去（時刻ｔ以前）のシステムの挙動とパラメータｚと時刻ｔとを与えられた漸化式（数１７）に代入して、時刻ｔでのシステムの状態ｘ（ｋ）（ｔ）を計算する。
　続いて、乱数重要度算出部２３０は、Ｎ個のパラメータｚ１，…，ｚNそれぞれについて、時刻ｔで生成した乱数ｗｊ （ｋ）のパラメータに関する乱数重要度Ｈｉ，ｊ （ｋ）を計算する（Ｓ５０３）。なお、シナリオ・シミュレータ２００は、ユーザ等が指定したパラメータについてのみ、ステップＳ５０３における処理を行ってもよい。
　ｉ番目のパラメータについての乱数重要度Ｈｉ、ｊ （ｋ）は、具体的には以下の３つのステップにより算出される、すなわちステップＳ５０３は以下の３つのステップを含む。なお、各パラメータについての等価微分重要度ｈｉ （ｋ）の計算は互いに独立であるため、１つのシナリオ・シミュレータ２００によってＮ個のパラメータｚ１，…，ｚNについて等価微分重要度の計算を逐次的に行うのではなく、例えば多数の計算機を用いて当該計算を同時に並列実行してもよい。
　第１のステップは、サンプルパス固定微分計算部２３１が、システムの挙動ｘを一定に保つような、乱数ｗｊ （ｋ）とｉ番目のパラメータの微小変化の関係Ｄｉ，ｊ （ｋ）＝（∂ｗｊ （ｋ）／∂ｚｉ）｜ｘ＿ｊ＋１＾（ｋ）＝ｃｏｎｓｔを求めるステップである。
　サンプルパス固定微分計算部２３１は、例えば、関数－（∂ｆ／∂ｚｉ）／（∂ｆ／∂ｗｊ （ｋ））が与えられた場合、入力されたパラメータｚｉ及び入力された乱数ｗｊ （ｋ）を当該関数に代入することにより、Ｄｉ，ｊ （ｋ）を算出すればよい。また、サンプルパス固定微分計算部２３１は、例えば、パラメータの微小変化Δｚｉに対して、ｆを一定に保つようなΔｗを数値的に探索してΔｗ／Δｚｉを計算することにより、（∂ｗｊ／∂ｚｉ）｜（x(t)＝ｃｏｎｓｔ）、すなわちＤｉ，ｊ （ｋ）を近似的に算出してもよい。サンプルパス固定微分計算部２３１は、算出したＤｉ，ｊ （ｋ）を、サンプルパス固定２階微分計算部２３２、及び乱数重要度計算部２３３に送信する。
　第２のステップは、サンプルパス固定２階微分計算部２３２が、受信したＤｉ，ｊ （ｋ）を、ｗｊ （ｋ）で微分し、ＤＤｉ，ｊ （ｋ）を算出するステップである。第２のステップは、ステップＳ４０２における第２のステップと同様である。
　第３のステップは、乱数重要度計算部２３３が、時刻ｔにおける乱数ｗｊ （ｋ）に関する乱数重要度Ｈｉ，ｊ （ｋ）を計算するステップである。第３のステップは、ステップＳ４０２における第３のステップと同様である。
　続いて、シミュレーション制御部２８０は、時刻をｔ＋Δｔにセットする（Ｓ５０４）。時刻ｔ＋Δｔがシミュレーション終了時刻ｔ＿ｅｎｄ以前である場合（Ｓ５０５：Ｙｅｓ）、モンテカルロ・シミュレーション全体制御部５００は、シナリオ・シミュレータ２００を再度起動し、シナリオ・シミュレータ２００は、ステップＳ５０２～ステップＳ５０４の処理を行う。
　時刻ｔ＋Δｔがシミュレーション終了時刻ｔ＿ｅｎｄより後である場合（Ｓ５０５：Ｎｏ）、指標算出部２６０は、モデル挙動記憶部２４０から、システムの挙動ｘ（ｋ）（ｔ）の時系列を取得し、注目指標の評価値、すなわちシステムの挙動の評価関数ａ（ｋ）を計算する（Ｓ５０６）。指標算出部２６０は、算出したａ（ｋ）を指標値記憶部３１０に格納する。
　続いて、等価微分重要度算出部２７０は、乱数重要度記憶部２５０から、全てのｊについてＨｉ，ｊ （ｋ）を取得し、ｊについての総和を計算することで、ｉ番目のパラメータに関する等価微分重要度ｈｉ （ｋ）を算出する（Ｓ５０７）。なお、等価微分重要度算出部２７０は、は、ユーザ等が指定したパラメータについてのみ、ステップＳ５０７における処理を行ってもよい。なお、ステップＳ５０６における処理とステップＳ５０７における処理は並行して行われてもよいし、順序を入れ替えて行われてもよい。重み関数が定める処理は、Ｓ５０３及びＳ５０７における処理を含む。
　以上のように、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される場合に、図５に示した処理により、シナリオ・シミュレータ２００は、システム挙動の評価関数ａ（ｋ）およびＮ個のパラメータそれぞれについて等価微分重要度ｈｉ （ｋ）を計算することができる。
　図６は、実施例１におけるモンテカルロ・シミュレーション・システムのユーザインタフェースの一例である。ユーザインタフェース６００は、シミュレーション条件設定セクション６１０、期待値計算セクション６２０、及び感度計算セクション６３０を含む。
　シミュレーション条件設定セクション６１０は、入力受付セクション６１１、６１２、及び表示セクション６１３を含む。入力受付セクション６１１は、シミュレーションモデル、すなわちシステムの挙動の実現値ｘを定める関数ｘ（ｗ，ｔ）等の入力を受け付ける。入力受付セクション６１２は、当該シミュレーションモデルによるシミュレーション回数の入力を受け付ける。表示セクション６１３は、入力受付セクション６１１～６１２が入力を受け付けた内容のシミュレーションを実行するために必要な計算時間を表示する。
　期待値計算セクション６２０は、入力受付セクション６２１、及び表示セクション６２２を含む。入力受付セクション６２１は、注目指標、すなわちシステムの挙動の評価関数ａ（ｘ）の入力を受け付ける。表示セクション６２２は、システム評価関数Ａ（ｘ）の計算精度を表示する。
　感度計算セクション６３０は、入力受付セクション６３２、６３４、チェックボックス６３１、６３３、及び表示セクション６３５を含む。入力受付セクション６３２は、感度計算の対象パラメータ、すなわち感度∂Ａ／∂ｚｉを算出する対象となるパラメータｚｉの入力を受け付ける。チェックボックス６３１は、全てのパラメータｚｉを感度計算の対象パラメータとして選択するためのチェックボックスである。入力受付セクション６３４は、Ｄｉ，ｊを計算する方法の入力を受け付ける。チェックボックス６３３は、Ｄｉ，ｊを計算する方法を自動で選択するためのチェックボックスである。表示セクション６３５は、感度∂Ａ／∂ｚｉの計算精度を表示する。
　なお、ユーザインタフェース６００における、入力セクション、表示セクションは、それぞれ表示セクション、入力セクションであってもよいし、それぞれ表示セクション、入力セクションを兼ねてもよい。例えば、表示セクション６２２がシステム評価関数Ａ（ｘ）の計算精度の入力を受け付け、表示セクション６３５が感度∂Ａ／∂ｚｉの計算精度の入力を受け付けてもよい。このとき、入力受付セクション６１１は、例えば、入力されたシステム評価関数Ａ（ｘ）の計算精度と、入力された感度∂Ａ／∂ｚｉの計算精度と、を保証するために必要なシミュレーション回数を表示してもよい。
　以上、本実施例のモンテカルロ・シミュレーション・システムは、１セットのモンテカルロ・シミュレーションの結果を用いて、システム評価関数の各パラメータに関する感度を算出することができ、ひいては、計算時間及び計算リソースを削減することができる。
　また、本実施例のモンテカルロ・シミュレーション・システムが、算出した感度を用いて、ユーザ等は、パラメータの最適化等、すなわちパラメータ値の許容される範囲で、システム評価関数を最大化するパラメータ値の探索等を行うことができる。
　本実施例は、多数のパラメータが存在する確率的なモデルにおいて、注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって同時に求める手段を、半導体集積回路の電子回路における製造ばらつきを考慮した設計に応用する例を説明する。
　本実施例は、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される、すなわち、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）が定まると、時刻ｔでのシステムの状態ｘ（ｔ）が、時刻ｔとｙのみで決定論的に定まる（図３に相当）一例である。
　近年、半導体集積回路の微細化が進むにしたがって、集積回路を構成する個々のトランジスタの特性のばらつきが大きくなっている。集積回路を構成するトランジスタの特性ばらつきは、例えば、集積回路の製造時において、形成される回路要素の物理的なサイズ、基板に注入される不純物イオンの数、又は格子欠陥の数等の製造時に不可避的に生じる確率的な現象に起因する。つまり、当該回路は、個々のトランジスタの特性が確率的にばらつくことを想定して設計される必要がある。
　したがって、トランジスタの特性のばらつきを確率分布でモデル化した上で、例えば、回路全体の特性の期待値やばらつき（分散）の大きさ、又は集積回路の歩留まり、すなわち集積回路を構成する回路全体の特性が望ましい範囲におさまる確率、等をモンテカルロ・シミュレーションによって推定することが行われる。回路設計者は、例えば、設計した回路のばらつきを考慮した特性が与えられた仕様を満たすまで、設計パラメータをチューニングしては、モンテカルロ・シミュレーションを行ってばらつきを考慮した特性を確認する、というプロセスを繰り返し行う。
　個々のトランジスタにおける設計パラメータは、例えば、ゲート幅Ｗ、ゲート長Ｌ、拡散層幅ＬＯＤ、トランジスタ間距離ＰＤＸといった、ばらつき・電力・レイアウト面積等に関係するものを１０程度含む。したがって、一般に、回路ブロックの設計パラメータの数は、回路を構成するトランジスタ数の１０倍程度になる。このため、多数のトランジスタで構成される複雑な回路ブロックにおける設計パラメータは膨大な数となりえるため、上記のように個々の設計パラメータのチューニングを人手で行うことは困難である。
　これに対して、本実施例における手法では、任意の注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって求めることができる。本実施例の手法を用いることで、回路設計者は、単に、現在のパラメータ設定のもとでの回路ブロック全体の特性（性能期待値や歩留まり）や電力の値のみならず、全てのトランジスタの全てのパラメータについて、それらが回路ブロック全体の特性（性能期待値や歩留まり）や電力にどのように影響するのか、その感度を同時に得ることができる。
　したがって、回路設計者は、例えば、どのトランジスタのどのパラメータを変えることで、消費電力を増やさずに特性を向上させることができるか等を、やみくもにパラメータをチューニングすることなく知ることができ、設計の効率を大幅に向上できる。
　図７は、６個のトランジスタで構成されるＳＲＡＭセルの回路の例を示す。ＳＲＡＭセルは、通常時にはＳＲＡＭセルに記憶された値が失われないことと、記憶の書き換え時にはＳＲＡＭセルの記憶値が意図した値にきちんと書き変わること、という互いに相反する特性を満たす必要がある。
　ＳＲＡＭに含まれる各トランジスタのパラメータは、集積回路の製造時のばらつきによって、製造された集積回路における個々のトランジスタの実際のパラメータ値は設計時に指定した値にはならず、例えば、ある確率分布にしたがって分布する。その結果、例えば、通常にＳＲＡＭセルに記憶していた値が失われる、又は、記憶の書き換え時にはＳＲＡＭセルの記憶値が意図した値にきちんと書き変わらない、等の不良が起こりえる。
　そこで、回路設計者は、設計パラメータ値をいろいろな値に変えてモンテカルロ・シミュレーションを行い、実際に製造されるＳＲＡＭセル回路の不良率が与えられた閾値以下になるように設計パラメータを最適化する必要がある。この例に、実施例１のモンテカルロ・シミュレーション・システムのシミュレーション手法を適用する。
　本実施例におけるパラメータベクトルｚは、トランジスタの設計パラメータからなる。例えば、個々のトランジスタに１０の設計パラメータがある場合、６個のトランジスタで構成されるＳＲＡＭセル回路におけるパラメータベクトルｚは、６０次元のベクトルである。
　ところで、前述のように、ＳＲＡＭ回路の特性ばらつきは、トランジスタの実際のパラメータ値が製造時にばらつくことに起因する。すなわち、トランジスタの特性を決める６０次元の設計パラメータｚに対して、実際には製造時のばらつきにより、例えば、６０次元のパラメータｙが生成される。
　乱数発生部２１０は、シミュレーション毎に確率密度関数ｇ（ｗ）に従う多次元の乱数ベクトルｗを生成する。モデル挙動更新部２２０は、乱数ベクトルｗと６０次元の設計パラメータベクトルｚと、から製造後の実際のパラメータ値ｙを関数ｙ＝ｙ（ｗ，ｚ）によって定める。このとき、製造後の実際のパラメータ値ｙの確率分布が、設計パラメータベクトルがｚであるときの、製造ばらつきを表すように、乱数ベクトルｗが従う確率密度関数ｇ（ｗ）と、関数ｙ（ｗ，ｚ）が定められているものとする。
　モデル挙動更新部２２０は、製造後の実際のパラメータｙを定めれば、ＳＲＡＭセル回路の挙動ｘ（ｔ）を、通常の回路シミュレーションによって決定論的に定めることができる。したがって、本実施例におけるシステムの挙動を表す時系列ｘ（ｔ）は、数１５の形で表される、すなわち、乱数ベクトルｗとパラメータベクトルｚで定まる値ｙ（ｗ,ｚ）を定めると、時刻ｔでのシステムの状態ｘ（ｔ）が，時刻ｔとｙのみで決定論的に定まる。また、システムの挙動の評価関数ａ（ｘ）を、そのセルの挙動が正常動作の範囲内であれば０、不良であれば１と定義すれば、システム評価関数Ａ（ｘ）、すなわち、ａ（ｘ）の期待値Ｅ［ａ（ｘ）］は、ＳＲＡＭセルの不良率を示す。
　したがって、モンテカルロ・シミュレーション・システムは、図４で説明した処理を行うことにより、ＳＲＡＭセルの不良率の６０個の設計パラメータ全てに対する感度を算出することができる。また、モンテカルロ・シミュレーション・システムは、は、算出した感度から、勾配法やニュートン法といった公知の最適化手法を用いて最急降下法ＳＲＡＭセル不良率を最小にするような６０個の設計パラメータ値の組み合わせを、算出してもよい。
　以上のように、本実施例のモンテカルロ・シミュレーション・システムは、多数のトランジスタで構成される半導体集積回路の電子回路における、製造ばらつきを考慮した設計パラメータ値の最適値を効率的に求めることができる。
　本実施例は、多数のパラメータが存在する確率的なモデルにおいて、注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって同時に求める手段を、金融分野、とくに、銀行の収益予測シミュレーションシステムに適用した例を説明する。本実施例における、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される、すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成する（時刻毎に異なる）乱数ｗｊと、パラメータｚと、で定まる。
　金融機関においては、将来の事業計画の立案時に、取り扱っている金融商品の将来の取扱高やその収益について計画を立てて、当該事業計画の妥当性を経営判断する。
　図８は、銀行における収支計画表の例を示す。収支計画表の最左列は、当該銀行における取扱金融商品の科目を示す。取扱金融商品は、例えば、企業向けの事業融資、個人向けの住宅ローン貸付、カードローン貸付、および銀行が自ら保有している有価証券配当について、それぞれＡからＣまでの３区分の科目を含む。
　実際の銀行においては、数千程度の金融商品の科目を保有することがある。表の各行は、最上行に示される期間における、金融商品の科目それぞれについての取扱高を示す。例えば、現時点が２０１４年２Ｑ終了時とすれば、‘１４／１Ｑおよび‘１４／２Ｑの取扱高は実績値を示すが、‘１４／３Ｑ以降の取扱高は現時点（２０１４年２Ｑ終了時点）における計画値を示す。図８の例では四半期単位で計６期間についての計画を示しているが、実際の銀行においては、例えば、月単位又は週単位で５年程度の計画を立てることもあり、このとき、収支計画表は数百程度の期間を含む。
　各科目における取扱高と収益との関係は、例えば、主に銀行の資金調達金利、すなわちスポットレートに依存する。スポットレートと収益の関係については、科目毎にある程度の確度をもったモデル化が可能であるのに対して、将来のスポットレートは不確定であり確率的な予測しかできない。
　そこで、銀行の経営においては、モンテカルロ・シミュレーションによって収益の期待値やばらつきをシミュレーションし、事業計画の妥当性を判断する。具体的には、例えば、乱数を用いて将来のスポットレートの実現パスを生成し、各科目の計画取扱高をもとに収益値を計算し、全科目の収益を合計することで銀行全体の収益値を計算する、という処理を多数回行う。当該多数回の処理によって、与えられた事業計画の下での銀行全体の収益の期待値やばらつきを推定する。
　ところで銀行において、事業計画の妥当性を判断するにあたっては、将来の各科目の取扱高の計画値が完全に達成された場合の銀行全体の期待収益やばらつきを求めるだけでは不足である。各科目の将来の実際の取扱高は、計画値からずれることが当然に想定されるためである。このとき、計画値と将来の実績がずれたときの影響について、事業計画立案時に予め考慮しておく必要がある。
　具体的には、例えば、数千程度の科目それぞれについての数百程度の未来の各期間における取扱高の計画値、すなわち、合計十万個程度の計画値それぞれについて、将来の実績値が計画値からずれたときの銀行全体の収益に対する影響を計算する必要がある。もし実績値が計画値からほんの少しだけずれるだけで、銀行全体の収益に大きな影響を及ぼすような計画となっていることが分かれば、適切な引当金を積んでおく、又は、事業計画そのものを立て直す、といった対策が必要になるであろう。
　以下、実施例１のモンテカルロ・シミュレーション・システムによるシミュレーション手法を銀行の収益予測シミュレーションに適用する例を説明する。
　本実施例におけるパラメータベクトルｚは、例えば、全ての科目それぞれについて未来の各期間における取扱高の計画値である。例えば、科目数は数千程度、期間は数百期間である場合、パラメータベクトルｚは十万次元程度のベクトルである。システムの挙動の時系列ｘ（ｔ）は、例えば、将来の各期間における各科目の収益値の集合（例えば、数千次元のベクトルになりえる）と、金利のスポットレートと、を合わせた多次元のベクトルである。
　ｔ番目の期間における各科目の収益値は、例えば、スポットレートと、その時刻における各科目の取扱高の計画値であるパラメータｚと、に依存する。将来のスポットレートは、例えば、Ｈｕｌｌ－Ｗｈｉｔｅモデル、ＣＩＲモデル、又はＢＤＴモデルなどといった確率的なモデルで表現される。これらのモデルは、時刻ｔ＋１でのスポットレートを、時刻ｔでのスポットレートと、（時刻ｔで新たに生成する）乱数と、に基づいて定める。
　したがって、本実施例におけるシステムの挙動ｘ（ｔ）は、全体として、数１７の形で表される。すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）は、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成される（時刻毎に異なる）乱数ｗｊと、パラメータｚと、からで定まる。
　システムの挙動の評価関数ａ（ｘ）は、例えば、全科目の全期間の収益の合計値である。このとき、システム評価関数Ａ（ｘ）、すなわちａ（ｘ）の期待値Ｅ［ａ（ｘ）］は、与えられた計画のもとでの銀行全体の将来の収益の期待値を示す。
　本実施例のモンテカルロ・シミュレーション・システムは、上述した前提条件のもとで、図５に示した処理を行うことで、例えば、数千程度ある全ての科目それぞれについて数百程度ある未来の各期間における取扱高の計画値、すなわち、合計十万個程度の計画値それぞれについて、銀行全体の収益の期待値に対する感度、を同時に推定することが可能になる。
　銀行は、推定した感度に基づいて、例えば、計画値と実績値のずれの影響を予め考慮して適切な引当金を積んでおく、又は事業計画そのものを立て直す、といった対策を立てることが可能となる。
　以上のように、本実施例のモンテカルロ・シミュレーション・システムは、金融分野、とくに銀行の収益予測シミュレーションシステムにおいて計画値の収益に対する感度を推定することが可能となるため、銀行は事業計画立案のプロセスを効率化することができる。
　本実施例のモンテカルロ・シミュレーション・システムは、実施例１の数９における境界残余Ｒｉを算出する。前述のように、境界残余Ｒｉは、乱数ベクトルｗの分布範囲Ωの境界面（表面）∂Ωで（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ＝０（ベクトル）であれば、Ｒｉ＝０となる。現実的にモンテカルロ・シミュレーションの考察対象となるシステムでは、この条件が成立して、Ｒｉ＝０となる場合が多いと考えられる。
　しかしながら、例えば、乱数ベクトルｗが一様分布乱数である場合などに、この条件が成立せず、具体的に境界残余Ｒｉを算出する必要にせまられる場合がある。境界残余Ｒｉの算出には、数９の１行目に基づく方法と、数９の３行目に基づく方法と、の２通りの方法が考えられる。
　図９は、乱数ベクトルｗが、区間［ａ，ｂ］上の一様分布乱数ベクトルｗ１と、区間［ｃ，ｄ］上の一様分布乱数ｗ２と、の２つの要素をもつ２次元のベクトルである場合の境界残余Ｒｉを算出する概念の例を示す。図中に矩形で示された範囲がΩである。確率密度関数ｇ（ｗ）は、以下の数１９で表される。
Figure JPOXMLDOC01-appb-M000021
　なお、本実施例における記憶部３００は、サンプルパス固定微分記憶部３３０（不図示）をさらに含む。サンプルパス固定微分記憶部３３０は、サンプルパス固定微分計算部２３１が算出したＤｉ，ｊを保持する。また、本実施例における後処理計算部４００は、プログラムである境界残余算出部４３０（不図示）をさらに含む。境界残余算出部４３０は、乱数ベクトルｗと、パラメータベクトルｚと、サンプルパス固定微分記憶部３３０から受信するＤｉ，ｊと、を入力として境界残余Ｒｉを算出し、算出した境界残余Ｒｉを期待値計算部４１０に送信する。
　以下、乱数ベクトルｗが従う確率密度関数ｇが数１９で表される例を用いて、境界残余Ｒｉの２通りの算出方法を説明する。境界残余算出部４３０が、境界残余Ｒｉを算出する第１の方法は、数９の１行目に基づくものである。図９の条件において、数９の１行目は、以下の数２０のように展開される。
Figure JPOXMLDOC01-appb-M000022
　したがって、境界残余算出部４３０は、ｗ１をΩの境界上の値であるｂに固定したうえでｗ２を動かしてそれぞれ（∂ｗ２／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値から、ｗ１を境界値ａに固定したうえでｗ２を動かしてそれぞれ（∂ｗ２／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値を減算した値を算出する。また、境界残余算出部４３０は、ｗ２を分布の境界上の値であるｄに固定したうえでｗ１を動かしてそれぞれ（∂ｗ１／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値から、ｗ２を境界値ｃに固定したうえでｗ１を動かしてそれぞれ（∂ｗ１／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値を減算したものと、を加算した値を算出する。
　境界残余算出部４３０は、算出した２つの値を加算することで、境界残余Ｒｉを算出する。なお、境界残余算出部４３０は、例えば、既知の数値積分の手法を用いて、上述した積分値それぞれを算出すればよい。
　境界残余算出部４３０が、境界残余Ｒｉを算出する第２の方法は、数９の３行目に基づくものである。数９の３行目は、乱数ベクトルの分布範囲Ω全域における積分であるから、境界残余算出部４３０は、以下の数２１より、モンテカルロ・シミュレーションを用いた期待値として境界残余Ｒｉを算出することができる。
Figure JPOXMLDOC01-appb-M000023
　すなわち、乱数発生部２１０が、与えられた確率密度関数ｇ（ｗ）に従う乱数を多数発生させる。境界残余算出部４３０は、発生した乱数ベクトルそれぞれに対して、数２１の２行目の角括弧内の値を計算し、計算した値の平均値をＲｉの推定値とすればよい。
　なお、境界残余算出部４３０は、数値的な微分により、数２１中のｄｉｖを算出すればよい。また、ユーザ等によって関数としてｄｉｖが与えられた場合、境界残余算出部４３０は、当該関数にパラメータベクトルｚと乱数ベクトルｗとを代入することにより、ｄｉｖを算出してもよい。なお、境界残余算出部４３０は、確率密度関数ｇ（ｗ）が、２次元一様分布以外である場合も同様に境界残余Ｒｉを算出することができる。感度計算部４２０は、上述の方法で算出された境界残余Ｒｉを、算出した平均値から引いた差を算出する。当該算出した差が当該感度の推定値である。
　以上のように、本実施例のモンテカルロ・シミュレーション・システムは、境界残余Ｒｉがゼロでな場合に、その値を算出することが可能となる。また、モンテカルロ・シミュレーション・システムは、境界残余Ｒｉを算出することにより、システムの挙動の実現値ｘを算出するモデルに適用する乱数ベクトルがどのような分布であっても、即ちｘを算出する全てのモデルに対して、本実施例の感度計算手法を適用することができる。
　なお、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。例えば、上記した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることも可能であり、また、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。
　また、上記の各構成、機能、処理部、処理手段等は、それらの一部又は全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリや、ハードディスク、ＳＳＤ（Ｓｏｌｉｄ　Ｓｔａｔｅ　Ｄｒｉｖｅ）等の記録装置、または、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に置くことができる。
　また、制御線や情報線は説明上必要と考えられるものを示しており、製品上必ずしも全ての制御線や情報線を示しているとは限らない。実際には殆ど全ての構成が相互に接続されていると考えてもよい。

Claims (15)

  1. 　複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムであって、
    　プロセッサと、
    　メモリと、を含み、
    　前記メモリは、
    　１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、
    　前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、
    　乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、
    　前記プロセッサは、
    　複数の乱数ベクトルと、パラメータベクトルと、を取得し、
    　前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、
    　前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、
    　前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、
    　前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、シミュレーションシステム。
  2. 　請求項１に記載のシミュレーションシステムであって、
    　前記重み関数は、
    　前記第１の微分係数それぞれの、前記第１の微分係数それぞれに対応する乱数による、第２の微分係数を決定し、
    　前記第２の微分係数に基づき前記パラメータに対する前記乱数ベクトルの重みを決定する、シミュレーションシステム。
  3. 　請求項２に記載のシミュレーションシステムであって、
    　前記重み関数における重みは、下記数式で表わされ、
    Figure JPOXMLDOC01-appb-M000001
     　上記数式における、ｗは前記乱数ベクトル、ｚは前記パラメータベクトル、Ｍは前記乱数ベクトルの次元、ｗｊは前記乱数ベクトルｗのｊ番目の要素、ｚｉは前記パラメータベクトルのｉ番目の要素、ｘは前記乱数ベクトルｗに対応する前記確率的システムの挙動値、ｇは前記乱数ベクトルｗが従う確率密度関数である、シミュレーションシステム。
  4. 　請求項２に記載のシミュレーションシステムであって、
    　前記取得した複数の乱数ベクトルは、一様分布に従い、
    　前記重み関数において、前記重みは、前記乱数ベクトルの全ての乱数の前記第２の微分係数の和である、シミュレーションシステム。
  5. 　請求項２に記載のシミュレーションシステムであって、
    　前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度は、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、の積の平均値で表わされる、シミュレーションシステム。
  6. 　請求項１に記載のシミュレーションシステムであって、
    　前記挙動関数の独立変数は、時刻及び時刻非依存の内部関数であり、
    　前記内部関数は、確率分布に従う乱数ベクトルと、パラメータベクトルと、に基づいて、第１ベクトルを決定し、
    　前記挙動関数は、一つの乱数ベクトルから、前記確率的システムの挙動の実現値の一つの時系列を決定し、
    　前記挙動評価関数は、前記一つの時系列に基づいて、前記確率的システムの挙動の一つの評価値を決定し、
    　前記重み関数は、前記内部関数による第１ベクトルを一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定し、
    　前記プロセッサは、
    　前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値の時系列を算出し、
    　前記算出した実現値の時系列それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出する、シミュレーションシステム。
  7. 　請求項１に記載のシミュレーションシステムであって、
    　前記挙動関数は、前記確率的システムの挙動の実現値の時系列を決定し、
    　前記挙動関数は、前記時系列における各時刻の実現値を、各時刻に対応付けられている乱数と各時刻より前の時刻の実現値とに基づき決定し、
    　前記挙動評価関数は、乱数ベクトルに対応する前記確率的システムの挙動の実現値の時系列に基づいて、前記乱数ベクトルに対応する前記確率的システムの挙動の評価値を決定し、
    　前記重み関数は、前記乱数ベクトルの乱数それぞれに対応する時刻における挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定し、
    　前記プロセッサは、
    　前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値の時系列を算出し、
    　前記算出した実現値の時系列それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出する、シミュレーションシステム。
  8. 　請求項５に記載のシミュレーションシステムであって、
    　前記プロセッサは、
    　前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値、前記取得した複数の乱数ベクトルが発生する確率、及び前記取得した複数の乱数ベクトルの乱数それぞれの第１の微分係数からなるベクトルと前記取得した複数の乱数ベクトルが従う確率密度関数の分布範囲の境界の外向き単位法線ベクトルとの内積、の積の前記境界上における積分値を算出し、
    　前記平均値と前記積分値とに基づいて、前記確率的システムの挙動の評価値の期待値の感度を算出する、シミュレーションシステム。
  9. 　請求項５に記載のシミュレーションシステムであって、
    　前記プロセッサは、
    　前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値、前記取得した複数の乱数ベクトルが発生する確率、及び前記取得した複数の乱数ベクトルの乱数それぞれ第１の微分係数からなるベクトル、の積の発散と、前記取得した複数の乱数ベクトルが発生する確率の逆数と、の積の第１の期待値を算出し、
    　前記平均値と前記第１の期待値とに基づいて、前記確率的システムの挙動の評価値の期待値の感度を算出する、シミュレーションシステム。
  10. 　請求項１に記載のシミュレーションシステムであって、
    　前記取得した複数の乱数ベクトルが従う確率密度関数は、前記取得したパラメータベクトルのパラメータそれぞれが微小変化する条件において、前記乱数ベクトルに対応する挙動関数の実現値を一定に保つ、乱数ベクトルの微小変化に対して、定義域が一定である、シミュレーションシステム。
  11. 　複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムにおけるシミュレーション方法であって、
    　前記シミュレーションシステムは、プロセッサと、メモリと、を含み、
    　前記メモリは、
    　１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、
    　前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、
    　乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、
    　前記方法は、
    　前記プロセッサが、
    　複数の乱数ベクトルと、パラメータベクトルと、を取得し、
    　前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、
    　前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、
    　前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、
    　前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、ことを含む方法。
  12. 　請求項１１に記載のシミュレーション方法であって、
    　前記重み関数は、
    　前記第１の微分係数それぞれの、前記第１の微分係数それぞれに対応する乱数による、第２の微分係数を決定し、
    　前記第２の微分係数に基づき前記パラメータに対する前記乱数ベクトルの重みを決定する、方法。
  13. 　請求項１１に記載のシミュレーション方法であって、
    　前記重み関数における重みは、下記数式で表わされ、
    Figure JPOXMLDOC01-appb-M000002
     　上記数式における、ｗは前記乱数ベクトル、ｚは前記パラメータベクトル、Ｍは前記乱数ベクトルの次元、ｗｊは前記乱数ベクトルｗのｊ番目の要素、ｚｉは前記パラメータベクトルのｉ番目の要素、ｘは前記乱数ベクトルｗに対応する前記確率的システムの挙動値、ｇは前記乱数ベクトルｗが従う確率密度関数である、方法。
  14. 　請求項１２に記載のシミュレーション方法であって、
    　前記取得した複数の乱数ベクトルは、一様分布に従い、
    　前記重み関数において、前記重みは、前記乱数ベクトルの全ての乱数の前記第２の微分係数の和である、方法。
  15. 　請求項１２に記載のシミュレーション方法であって、
    　前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度は、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、の積の平均値で表わされる、方法。
PCT/JP2014/072132 2014-08-25 2014-08-25 シミュレーションシステム、及びシミュレーション方法 WO2016030938A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
PCT/JP2014/072132 WO2016030938A1 (ja) 2014-08-25 2014-08-25 シミュレーションシステム、及びシミュレーション方法
US15/502,303 US11003810B2 (en) 2014-08-25 2014-08-25 Simulation system and simulation method
JP2016545097A JP6219528B2 (ja) 2014-08-25 2014-08-25 シミュレーションシステム、及びシミュレーション方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2014/072132 WO2016030938A1 (ja) 2014-08-25 2014-08-25 シミュレーションシステム、及びシミュレーション方法

Publications (1)

Publication Number Publication Date
WO2016030938A1 true WO2016030938A1 (ja) 2016-03-03

Family

ID=55398884

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/072132 WO2016030938A1 (ja) 2014-08-25 2014-08-25 シミュレーションシステム、及びシミュレーション方法

Country Status (3)

Country Link
US (1) US11003810B2 (ja)
JP (1) JP6219528B2 (ja)
WO (1) WO2016030938A1 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10445442B2 (en) * 2016-09-01 2019-10-15 Energid Technologies Corporation System and method for game theory-based design of robotic systems
US10579754B1 (en) * 2018-09-14 2020-03-03 Hewlett Packard Enterprise Development Lp Systems and methods for performing a fast simulation
WO2020101971A1 (en) * 2018-11-13 2020-05-22 Seddi, Inc. Procedural model of fiber and yarn deformation
US11398072B1 (en) * 2019-12-16 2022-07-26 Siemens Healthcare Gmbh Method of obtaining a set of values for a respective set of parameters for use in a physically based path tracing process and a method of rendering using a physically based path tracing process

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007213172A (ja) * 2006-02-07 2007-08-23 Japan Aerospace Exploration Agency モンテカルロ評価における危険不確定パラメータの検出方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3966139B2 (ja) 2002-09-27 2007-08-29 株式会社日立製作所 気象物理量の推定方法
US7653522B2 (en) * 2005-12-07 2010-01-26 Utah State University Robustness optimization system
US9201993B2 (en) * 2011-05-11 2015-12-01 Apple Inc. Goal-driven search of a stochastic process using reduced sets of simulation points

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007213172A (ja) * 2006-02-07 2007-08-23 Japan Aerospace Exploration Agency モンテカルロ評価における危険不確定パラメータの検出方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
NIKOLAIDIS, E. ET AL.: "Probabilistic Reanalysis Using Monte Carlo Simulation", SAE TECHNICAL PAPER SERIES, vol. 15, pages 1 - 14 *

Also Published As

Publication number Publication date
US11003810B2 (en) 2021-05-11
JPWO2016030938A1 (ja) 2017-06-15
JP6219528B2 (ja) 2017-10-25
US20170235862A1 (en) 2017-08-17

Similar Documents

Publication Publication Date Title
Mirhoseini et al. Chip placement with deep reinforcement learning
Mirhoseini et al. A graph placement methodology for fast chip design
Singaravel et al. Deep-learning neural-network architectures and methods: Using component-based models in building-design energy prediction
JP7520995B2 (ja) ニューラルネットワークを使った集積回路配置の生成
Couckuyt et al. Fast calculation of multiobjective probability of improvement and expected improvement criteria for Pareto optimization
Srivastava et al. Statistical analysis and optimization for VLSI: Timing and power
TWI559156B (zh) 識別稀有事件故障率的方法與系統
Xu et al. Learning viscoelasticity models from indirect data using deep neural networks
US9524365B1 (en) Efficient monte carlo flow via failure probability modeling
CN109426698A (zh) 预测半导体集成电路良率的装置和半导体器件的制造方法
Yondo et al. A review of surrogate modeling techniques for aerodynamic analysis and optimization: current limitations and future challenges in industry
JP6219528B2 (ja) シミュレーションシステム、及びシミュレーション方法
Savari et al. NN-SSTA: A deep neural network approach for statistical static timing analysis
Cook et al. Robust airfoil optimization and the importance of appropriately representing uncertainty
EP2051176A1 (en) Parametric yield improvement flow incorporating sigma to target distance
US8813009B1 (en) Computing device mismatch variation contributions
US10803218B1 (en) Processor-implemented systems using neural networks for simulating high quantile behaviors in physical systems
WO2016194051A1 (ja) 確率的システムの注目指標の統計量を最小化するパラメータセットを探索するシステム
Zhang et al. Multi-objective optimization for design under uncertainty problems through surrogate modeling in augmented input space
Burnicki et al. Propagating error in land-cover-change analyses: impact of temporal dependence under increased thematic complexity
Kong et al. New machine learning techniques for simulation-based inference: InferoStatic nets, kernel score estimation, and kernel likelihood ratio estimation
Sheta et al. Software effort and function points estimation models based radial basis function and feedforward artificial neural networks
El-Morshedy et al. A New Probability Heavy‐Tail Model for Stochastic Modeling under Engineering Data
Brusamarello et al. Fast and accurate statistical characterization of standard cell libraries
CN114997060A (zh) 一种声子晶体时变可靠性测试方法、计算设备及存储介质

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2016545097

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14900907

Country of ref document: EP

Kind code of ref document: A1