WO2020170867A1 - 推定装置、推定方法、及びプログラム - Google Patents

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

Info

Publication number
WO2020170867A1
WO2020170867A1 PCT/JP2020/004909 JP2020004909W WO2020170867A1 WO 2020170867 A1 WO2020170867 A1 WO 2020170867A1 JP 2020004909 W JP2020004909 W JP 2020004909W WO 2020170867 A1 WO2020170867 A1 WO 2020170867A1
Authority
WO
WIPO (PCT)
Prior art keywords
data
observed
sample
probability distribution
model
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
PCT/JP2020/004909
Other languages
English (en)
French (fr)
Inventor
匡宏 幸島
達史 松林
浩之 戸田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
NTT Inc
Original Assignee
Nippon Telegraph and Telephone Corp
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 Nippon Telegraph and Telephone Corp filed Critical Nippon Telegraph and Telephone Corp
Priority to US17/431,791 priority Critical patent/US20220138375A1/en
Publication of WO2020170867A1 publication Critical patent/WO2020170867A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; 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 OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N99/00—Subject matter not provided for in other groups of this subclass
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present invention relates to an estimation device, an estimation method, and a program, and more particularly to an estimation device, an estimation method, and a program for estimating parameters of a mixed model from censored data.
  • -Censored data refers to data for which the value of a sample whose observed value is greater than or equal to a certain threshold value (or less than a certain threshold value) is not observed and only the information above the threshold value is obtained.
  • a lot of data such as clinical data describing the onset of illness and death of people, contract history data of Internet line users, service usage history data of e-commerce sites, etc. are expressed as censored data.
  • a typical example of a problem using censored data is a survival time analysis that estimates a distribution of time required for a certain device to fail, for example.
  • the survival time distribution is often multimodal due to the presence of early failures and deterioration failures, so distribution estimation using a mixed model is widely used.
  • the first is the existence of local optimal solutions. That is, although the EMCM guarantees a monotonic decrease of the objective function, the convergence destination changes depending on the initial value, so that there is a problem that repeated execution from different initial values is required.
  • the second is the need to calculate the statistics of the truncated distribution of the probability distribution used in the model (the distribution that has been modified so that it takes values only within a certain range). That is, except for exceptional one-dimensional cut normal distribution which is a cut distribution of one-dimensional normal distribution, this statistic cannot be calculated analytically, and numerical calculation such as Monte Carlo method can be used. There was a problem that it was necessary. Since EMCM is an algorithm that repeats parameter updating and statistic calculation many times, it is considered desirable to avoid repeating numerical calculation at each iteration.
  • the present invention has been made in view of the above points, and provides an estimation device, an estimation method, and a program capable of accurately estimating a parameter of a model representing a probability distribution of censored data while suppressing a calculation time.
  • the purpose is to
  • the estimation device includes observation data of a sample in which an observation value is observed, observation data of a sample in which an observation value is not observed, and a variable indicating whether or not the observation value is observed in each sample.
  • An estimation device that estimates a parameter of a model that represents a probability distribution of censored data, and an input unit that receives an input of the censored data, and a distribution of observed values corresponding to each sample of the censored data that the input unit has received.
  • the estimation method for estimating a parameter of a model representing a probability distribution of censored data, wherein an input unit receives an input of the censored data, and a parameter estimation unit, each sample of the censored data received by the input unit.
  • a model representing the probability distribution of the censored data and the input unit accepts
  • the parameters of the model are estimated by optimizing the objective function, which is the divergence with the probability distribution of the censored data, obtained from the censored data.
  • the program according to the present invention stores the observation data of the sample in which the observation value is observed, the observation data of the sample in which the observation value is not observed, and the variable representing whether or not the observation value is observed in each sample.
  • each sample of data represented by a mixed model of each component representing the distribution of observed values, represented using a probability density function of the observed data, a model representing the probability distribution of the censored data, and
  • a program that causes a computer to execute a process including estimating parameters of the model by optimizing an objective function that is obtained from the censored data received by the input unit and is a divergence with the probability distribution of the censored data. is there.
  • the input unit has the observation data of the sample in which the observation value is observed, the observation data of the sample in which the observation value is not observed, and the observation value for each sample.
  • Input of censored data including a variable indicating whether or not observed is accepted.
  • the parameter estimation unit which corresponds to each sample of the censored data accepted by the input unit, is represented by a mixed model of each component that represents the distribution of observed values, and is represented using a probability density function of observed data,
  • the parameters of the model are estimated by optimizing the objective function which is the divergence between the model representing the probability distribution of the censored data and the probability distribution of the censored data obtained from the censored data received by the input unit.
  • each of the censored data including the observed data of the sample in which the observed value was observed, the observed data of the sample in which the observed value was not observed, and the variable indicating whether or not the observed value was observed for each sample
  • the model that represents the probability distribution of the censored data corresponds to the sample, represented by the mixture model of each component that represents the distribution of the observed values, represented by using the probability density function of the observation data, the model that represents the probability distribution of the censored data, and obtained from the censored data
  • the model representing the probability distribution of the censored data of the estimation device is represented by using a probability density function of the observation data and a length of time until the end of the observation given in advance for each sample.
  • the probability distribution of the variable, the probability density function of the observation data, and the length of time to the end of observation given in advance for each sample And the probability distribution of the observation data.
  • the objective function of the estimation device may be Kullback-Leibler divergence or L 2 divergence between the model representing the probability distribution of the censored data and the probability distribution of the censored data.
  • the estimation device the estimation method, and the program of the present invention, it is possible to accurately estimate the parameters of the model representing the probability distribution of the censored data, while suppressing the calculation time.
  • the first is an estimation method based on Kullback-Leibler (KL: Kullback-Leibler) divergence minimization between probability distributions
  • the second is an estimation method based on L 2 divergence minimization between probability distributions.
  • the first estimation method allows parameter estimation with a very simple iterative calculation
  • the second estimation method allows parameter estimation analytically without even needing iteration.
  • a key point in constructing the method according to the embodiment of the present invention is the use of an approach (reference 1) called an exemplar based model (hereinafter, emb).
  • an approach hereinafter, emb.
  • the parameters of each component of the mixed model are not treated explicitly as parameters, but each component is placed at the point where data points exist.
  • an algorithm that converges to a global optimal solution can be constructed by formulating only the mixture ratio (weight parameter of each component).
  • the above approach can be viewed as an evolution of this approach to use censored data as input.
  • the method of the reference proposes a method of performing emb by KL divergence minimization using normal data, but it cannot handle censored data.
  • FIG. 1 shows an example of one-dimensional cutoff data. As shown in the example of FIG. 1, description will be given using device failure data, which is a typical example of termination data.
  • the time when each device was installed and the time of the failure when the failure occurred are recorded.
  • the time is recorded because a failure occurs during the observation period.
  • the devices 3 and 4 did not fail during the observation period, and the observation was terminated, so that the time of the failure is not recorded.
  • both the device 3 and the device 4 always fail sometime, and the time of the failure is after the observation end time.
  • data consisting of a set of data such as the devices 1 and 2 whose observed values (failure times) are known and devices 3 and 4 whose observed values (failure time) are equal to or greater than a certain value. Is called censored data.
  • the device failure data in Fig. 1 is one-dimensional censored data, but this method can handle censored data of two or more dimensions, so it will be described here.
  • 2 is two-dimensional censored data that represents the usage period of users of a certain two services.
  • User 1 started using both services at the same time and canceled at the same time during the observation period, so the cancellation times of both services are recorded.
  • User 2 starts to use both services at the same time and cancels only service 2 during the observation period.
  • User 3 starts using both services at the same time and cancels only service 1 during the observation period.
  • the user 4 first starts using the service 2, and additionally contracts the service 1 during the observation period. Therefore, due to the termination of the observation, the cancellation time of the service 1 of the user 2, the cancellation time of the service 2 of the user 3, and the cancellation times of the services 1 and 2 of the user 4 are not recorded.
  • FIG. 3 shows the observation data of FIG. 1 expressed by survival time
  • FIG. 4 shows the observation data of FIG. 2 expressed by survival time. Since the data is defined as multidimensional censored data, FIG. 4 will be used as an example for description below.
  • Is. x ij represents the usage time of the service j of the user i
  • v ij represents the length from the service j use start time of the user i to the observation end time.
  • K is the number of mixtures
  • Probability distribution of components For example, a Gaussian distribution represented by the following formula can be used.
  • ⁇ represent the mean and standard deviation of the Gaussian distribution.
  • ⁇ represents a value randomly generated from a probability distribution (eg, exponential distribution) having a value of 0 or more.
  • a probability distribution e.g, exponential distribution
  • the standard deviation ⁇ can be determined by a cross-validation method or the like.
  • f tr is a truncated distribution of a distribution in which f is marginalized with respect to an unobserved element, and is represented by the following formula.
  • f tr complies with the following equations (5) to (6).
  • L 2 divergence is defined as the squared error of two probability density functions. Which divergence to use depends on the problem. Therefore, in the embodiment of the present invention, two types of methods are constructed so that they can be applied regardless of which divergence is used.
  • the constraint condition (the parameter element is 0 or more and the sum is 1) is for the mixed model f to have a probability distribution.
  • optimization can be performed by repeating the update of.
  • ⁇ k is smaller than a certain threshold value (for example, 10 ⁇ 3 /n)
  • ⁇ k 0.
  • a renormalization operation may be performed in which the sum is adjusted to be 1.
  • the objective function can be designed as in the following equation by replacing the place where the KL divergence/log likelihood ratio is used in each of these two terms with L 2 divergence.
  • Is a value that can be analytically calculated when a Gaussian distribution is used as the component distribution, and can be expressed as in the following equation (21).
  • the numerical problem solver can directly solve the optimization problem of the above formula (22). In that case, you may decide to solve the following optimization problem which added the regularization term.
  • ⁇ represents a hyper parameter.
  • the following approximate method may be used.
  • the optimum solution of the problem in which the constraint condition is removed from the optimization problem of the above formula (22) is obtained as in the following formulas (23) and (24).
  • Is a regularization term and ⁇ is a hyperparameter, which has an effect of preventing the divergence of the parameter.
  • the parameter of the model is estimated by one of the above-described two methods, thereby suppressing the calculation time and accurately estimating the parameter of the model representing the probability distribution of the censored data.
  • FIG. 5 is a block diagram showing a schematic configuration of a computer that functions as the estimation device 1 according to the embodiment of the present invention.
  • FIG. 6 is a block diagram showing the configuration of the estimation device 1 according to the embodiment of the present invention.
  • the estimation device 1 includes a CPU 110, a memory 120 such as a RAM, a communication interface (IF) unit 130, an input unit 140 such as a keyboard, a display unit 150 such as a display, and an estimation described later.
  • the computer is provided with a storage unit 160 such as a ROM that stores a program 170 for executing the processing routine.
  • the CPU 110, the memory 120, the communication IF unit 130, the input unit 140, the display unit 150, and the storage unit 160 are connected via the bus 100.
  • the communication IF unit 130 is connected to the external device 2 by a communication line such as a LAN cable. Note that the external device 2 may be connected via a network (not shown).
  • the estimation device 1 includes a data processing unit 10, a parameter estimation unit 20, a parameter output unit 30, a storage unit 40, an input unit 50, and an output unit 60. Prepared for.
  • the data processing unit 10 is the observation data of the sample in which the observation value received by the input unit 50 is observed. And the observation data of the sample for which the observation value was not observed And a variable that indicates whether the observed value was observed for each sample Censored data including and Are stored in the data storage unit 41.
  • the sample is, for example, each device in the examples of FIGS. 1 and 3 and each user in the examples of FIGS.
  • the parameter estimation unit 20 uses the censored data received by the input unit 50.
  • the censored data received by the input unit 50 Censoring data obtained from True probability distribution of By optimizing the objective function, which is the divergence with To estimate.
  • the parameter estimation unit 20 first determines the censored data True probability distribution of Ask for.
  • the parameter estimation unit 20 sets the objective function to the censored data.
  • Model representing the probability distribution of And the censored data True probability distribution of The parameter is estimated as KL divergence with or L 2 divergence.
  • a model that represents the probability distribution of censored data Is the probability distribution function of the observation data, which is expressed using the probability density function of the observation data and the length of time until the end of the observation, which is given in advance for each sample, as shown in the above equation (8). And a probability density function of the observation data, and a probability distribution of the variables under the given variables, which is represented by the length of time to the end of the observation given in advance for each sample. And are represented using.
  • the parameter estimation unit 20 estimates the parameter by repeating the parameter update of the equation (15). Further, the parameter estimator 20, when using a L 2 divergence, the equation (24), to estimate parameters using (25).
  • the parameter estimation unit 20 calculates the estimated parameter Is stored in the parameter storage unit 42.
  • the parameter output unit 30 is the parameter of the parameter storage unit 42. Is obtained, and the parameter is output to the parameter output unit 30. give.
  • the storage unit 40 includes a data storage unit 41 and a parameter storage unit 42.
  • the data storage unit 41 stores censored data. Is stored.
  • the parameter storage unit 42 stores the model parameters. Is stored.
  • the input unit 50 receives the censored data from the external device 2. Accepts input. Then, the input unit 50 receives the censored data To the data processing unit 10.
  • the output unit 60 is the model parameter received from the parameter output unit 30. Is output to the external device 2.
  • FIG. 7 is a flowchart showing an estimation processing routine according to the embodiment of the present invention.
  • step S100 the input unit 50 sets the observation data of the sample in which the observation value is observed. And the observation data of the sample for which the observation value was not observed And a variable that indicates whether the observed value was observed for each sample Censored data including and Accepts input.
  • the input unit 50 also receives an input of the length of time until the observation end time for each sample i.
  • step S110 the parameter estimation unit 20 sets the objective function to the censored data.
  • Model representing the probability distribution of And the censored data Probability distribution of As KL divergence or L 2 divergence, the parameters are estimated by repeating the parameter update of the above equation (15), or the parameters are estimated using the above equations (24) and (25).
  • step S120 the output unit 60 uses the parameters estimated in step S110. Is output.
  • the observed data of the sample in which the observed value is observed, the observed data of the sample in which the observed value is not observed, and the observed value for each sample are Corresponding to each sample of censored data including a variable representing whether or not observed, represented by a mixed model of each component representing the distribution of observed values, represented using the probability density function of the observed data,
  • the parameters of the model are estimated, so that the calculation time is suppressed. , It is possible to accurately estimate the parameters of the model that represents the probability distribution of censored data.
  • the present invention is not limited to this and the present invention can be applied to any censored data that is not time-series data. It can be applied.
  • estimation device 1 has been described as a configuration in which the processing of each unit is constructed as a program and installed and executed in a computer used as the estimation device, but as a configuration to be distributed via a network. Good.
  • the program can be stored in a computer-readable recording medium and provided.

Landscapes

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

Abstract

計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができるようにする。　パラメタ推定部２０が、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定する。

Description

推定装置、推定方法、及びプログラム
　本発明は、推定装置、推定方法、及びプログラムに係り、特に、打ち切りデータから混合モデルのパラメタを推定するための推定装置、推定方法、及びプログラムに関する。
　打ち切りデータとは、観測値がある閾値以上（またはある閾値以下）であるサンプルについては、値が観測されず、閾値以上である、という情報しか得られないデータのことを指す。病気の発症や人の死亡などを記述する臨床データや、インターネット回線利用者の契約履歴データ、Ｅコマースサイトのサービス利用履歴データなど多くのデータが打ち切りデータとして表現される。打ち切りデータを用いる問題の代表例は、例えばある機器が故障するまでに要する時間の分布を推定する、生存時間分析である。生存時間分布は、初期故障や劣化故障などの存在によりしばしば多峰性をもつため、混合モデルを用いた分布推定が広く用いられる。
　打ち切りデータから混合モデルのパラメタを推定するには、文献（非特許文献１）で提案されているＥＭＣＭ（Ｅｘｐｅｃｔａｔｉｏｎ　Ｍａｘｉｍｉｚａｔｉｏｎ　ｆｏｒ　Ｃｅｎｓｏｒｅｄ　Ｍｉｘｔｕｒｅｍｏｄｅｌｓ）アルゴリズムが適用できる。
Didier Chauveau. A stochastic em algorithm for mixtures with censored data. Journal of statistical planning and inference, Vol.46, No.1, 1995, pp.1-25.
　しかし、このアプローチには以下の２つの問題があった。
　１つ目は、局所最適解の存在である。すなわち、ＥＭＣＭは目的関数の単調減少は保証されるものの収束先は初期値に依存して変わるため、異なる初期値からの繰り返し実行が必要となってしまう、という問題があった。
　２つ目は、モデルで利用する確率分布の切断分布（ある一定の範囲でしか値をとらないように変更した分布）の統計量の計算の必要性である。すなわち、１次元正規分布の切断分布である１次元切断正規分布などの例外的なものを除き、この統計量は解析的に計算することができず、モンテカルロ法などの数値計算を利用することが必要となる、という問題があった。ＥＭＣＭは、パラメタの更新と統計量の計算を何度も繰り返すアルゴリズムであるから、各反復で数値計算を繰り返すことは避けることが望ましいと考えられる。
　本発明は上記の点に鑑みてなされたものであり、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる推定装置、推定方法、及びプログラムを提供することを目的とする。
　本発明に係る推定装置は、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定装置であって、前記打ち切りデータの入力を受け付ける入力部と、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定するパラメタ推定部と、を備えて構成される。
　また、本発明に係る推定方法は、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定方法であって、入力部が、前記打ち切りデータの入力を受け付け、パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する。
　また、本発明に係るプログラムは、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する処理をコンピュータに実行させるプログラムであって、入力部が、前記打ち切りデータの入力を受け付け、パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定することを含む処理をコンピュータに実行させるプログラムである。
　本発明に係る推定装置、推定方法及びプログラムによれば、入力部が、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの入力を受け付ける。
　そして、パラメタ推定部が、入力部が受け付けた打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、打ち切りデータの確率分布を表すモデルと、入力部が受け付けた打ち切りデータから得られる、打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定する。
　このように、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定することにより、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
　また、本発明に係る推定装置の前記打ち切りデータの確率分布を表すモデルは、前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数の確率分布と、前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数が与えられた下での、前記観測データの確率分布と、を用いて表されることができる。
　また、本発明に係る推定装置の前記目的関数は、前記打ち切りデータの確率分布を表すモデルと、前記打ち切りデータの確率分布とのカルバックライブラーダイバージェンス、又はＬ２ダイバージェンスであることができる。
　本発明の推定装置、推定方法、及びプログラムによれば、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
１次元の打ち切りデータの例を表すイメージ図である。 ２次元の打ち切りデータの例を表すイメージ図である。 １次元の打ち切りデータの生存時間表現の例を表すイメージ図である。 ２次元の打ち切りデータの生存時間表現の例を表すイメージ図である。 本発明の実施の形態に係る推定装置として機能するコンピュータの概略構成を示すブロック図である。 本発明の実施の形態に係る推定装置の構成を示すブロック図である。 本発明の実施の形態に係る推定装置の推定処理ルーチンを示すフローチャートである。
　以下、本発明の実施の形態について図面を用いて説明する。
＜本発明の実施の形態に係る推定装置の原理＞
　まず、本発明の実施形態の原理について説明する。
　本発明の実施の形態では、次に示す目的関数の異なる２つの手法を構築した。１つ目は確率分布間のＫｕｌｌｂａｃｋ－Ｌｅｉｂｌｅｒ（ＫＬ：カルバック－ライブラー）ダイバージェンス最小化に基づく推定手法、２つ目は確率分布間のＬ２ダイバージェンス最小化に基づく推定手法である。１つ目の推定手法では非常にシンプルな繰り返し計算でパラメタの推定が可能であり、２つ目の推定手法では、繰り返しの必要すらなく解析的にパラメタの推定が可能である。
　本発明の実施の形態に係る手法を構築する上でポイントとなるのは、ｅｘａｍｐｌａｒ　ｂａｓｅｄ　ｍｏｄｅｌ（以下、ｅｍｂ）と呼ばれるアプローチ（参考文献１）の利用である。
［参考文献１］Danial Lashkari and Polina Golland, "Convex clustering with exemplar-based models", In Advances in neural information processing systems, 2008, pp.825-832.
　ｅｍｂアプローチでは、混合モデルの各コンポーネントのパラメタを陽にパラメタとしては扱わず、データ点の存在する点に各コンポーネントを配置する。これによって、混合比（各コンポーネントの重みパラメタ）のみの推定を行う定式化をすることで大域的最適解へ収束するアルゴリズムが構築できる。上記の手法は、このアプローチに対して打ち切りデータを入力として利用できるよう発展させたものとみなすことができる。
　更に、上述したＥＭＣＭで必要となる切断分布の統計量はコンポーネントのパラメタの推定に必要なものであったため、このアプローチの採用により数値計算の繰り返し実行の必要のない手法を構築することができる。
　なお、文献（参考文献１）の方法では、通常のデータを用いてＫＬダイバージェンス最小化によるｅｍｂを行う手法を提案しているが、打ち切りデータを扱うことができない。
　確率分布の誤差の関数にＬ２ダイバージェンスを用いるアプローチは文献（参考文献２）でも考えられているが、打ち切りデータではなく通常のデータを入力する状況を考えている。
［参考文献２］DavidW Scott, "Parametric statistical modeling by minimum integrated square error", Technometrics, Vol.43, No.3, 2001, pp.274-285.
　また、打ち切りデータを利用し、かつＬ２ダイバージェンスを用いる研究には文献（参考文献３）が存在するが、指数分布、ワイブル分布のように単峰の単純な分布をモデルとして用いる場合のみが考えられている。
［参考文献３］Srabashi Basu, Ayanendranath Basu, and MCJones, "Robust and efficient parametric estimation for censored survival data", Annals of the Institute of Statistical Mathematics, Vol.58, No.2, 2006, pp.341-355.
＜＜準備＞＞
＜＜＜打ち切りデータ＞＞＞
　まず、打ち切りデータについて説明する。図１に、１次元の打ち切りデータの例を示す。図１の例に示すように、打ち切りデータの代表例である機器故障データを用いて説明する。
　機器故障データでは、各機器の設置された時刻と故障が起こった際の故障の時刻が記録されている。機器１及び２に関しては観察期間中に故障が起きているため、その時刻が記録されている。一方、機器３及び４に関しては観察期間中には故障せずに、観察が打ち切られたため、故障の時刻が記録されてはいない。
　しかしながら、機器３及び機器４の何れも、いつかは必ず故障するとして、故障の時刻は観察終了時刻以降にあるということは読み取れる。このように機器１と２のように観測値（故障の時刻）が分かるデータと機器３と４のように観測値（故障の時刻）がある値以上であることが分かるデータの組からなるデータを、打ち切りデータと呼ぶ。
　図１の機器故障データは、１次元の打ち切りデータであるが、本手法では２次元以上の打ち切りデータも扱うことができるためここで説明する。
　図２は、ある２つのサービスの利用者の利用期間を表す２次元の打ち切りデータである。図２のケースでは、利用者の（少なくとも１つの）サービスを利用開始した時刻と各サービスを追加契約または解約した時刻が記録されている。
　利用者１は両サービスを同時に利用開始し、観察期間中に同時に解約しているため両方のサービスの解約時刻が記録されている。利用者２は両サービスを同時に利用開始し、観察期間中にサービス２のみを解約している。利用者３は両サービスを同時に利用開始し、観察期間中にサービス１のみを解約している。利用者４はサービス２をまず利用開始し、観察期間中にサービス１を追加契約している。よって、観測の打ち切りによって、利用者２のサービス１の解約時刻、利用者３のサービス２の解約時刻、利用者４のサービス１と２の解約時刻は記録されていない。
　このように２次元の打ち切りデータにおいては、打ち切られている値の次元が異なる、３通りの打ち切りの種類が存在する。一般にｎ次元の打ち切りデータでは２ｎ－１種類の打ち切りが存在する。なお、以下、各次元に観測が打ち切られるかが決まる例を用いて説明するが、どれか一つの要素が打ち切られたときに全ての要素が観測されなくなるという状況でも同様のアプローチで対応することができる。
　ここで、打ち切りデータの定義を与える。扱いを簡単にするため、観測データは図１のようなカレンダー時刻ではなく、生存時間（機器の設置から故障までの時間、サービスの契約から解約までの時間）を用いて表現する。図１の観測データを生存時間で表現したものを図３、図２の観測データを生存時間で表現したものを図４に示す。多次元の打ち切りデータとして定義するため、以下では図４を説明の例に用いる。
　打ち切りデータを
Figure JPOXMLDOC01-appb-I000001

と書く。ここで、
Figure JPOXMLDOC01-appb-I000002

と、
Figure JPOXMLDOC01-appb-I000003

とは共にｄｘ次元のベクトルであり、
Figure JPOXMLDOC01-appb-I000004

、
Figure JPOXMLDOC01-appb-I000005

である。ｘｉｊが利用者ｉのサービスｊの利用時間、ｗｉｊが利用者ｉのサービスｊの解約時刻が記録されたか（ｗｉｊ＝１）、打ち切りにより記録されなかったか（ｗｉｊ＝０）を表す。同様に、ｉ番目の利用者の、観測終了時刻までの時間の長さを
Figure JPOXMLDOC01-appb-I000006

と書く。ｖｉｊが利用者ｉのサービスｊ利用開始時刻から観測終了時刻までの長さを表す。打ち切りにより観測値が観測されなかったとき（ｗｉｊ＝０）には、ｘｉｊ＝ｖｉｊと設定されているとする。
＜＜＜混合モデル＞＞＞
　次に、本発明の実施の形態で用いるモデルについて説明する。混合モデルで表される、観測値の確率密度関数は一般に下記式（１）で定義される。
Figure JPOXMLDOC01-appb-M000007
　ここで、Ｋは混合数、
Figure JPOXMLDOC01-appb-I000008

がｋ番目コンポーネントの確率分布を表す。コンポーネントの確率分布
Figure JPOXMLDOC01-appb-I000009

には、例えば下記式で表されるガウス分布が利用できる。
Figure JPOXMLDOC01-appb-M000010
　ここで、
Figure JPOXMLDOC01-appb-I000011

とσとは、ガウス分布の平均と標準偏差を表す。ただし、ｅｂｍ（参考文献１）のアプローチに従い、コンポーネントの確率分布のパラメタ
Figure JPOXMLDOC01-appb-I000012

はＫ＝ｎとして、それぞれが観測データ点に対応するよう
Figure JPOXMLDOC01-appb-I000013

と設定されたもの等であるとする。この場合には、観測値の確率密度関数は、観測データの各々に対する、各コンポーネントの混合モデルで表され、各コンポーネントに含まれる、ガウス分布の平均を、対応する観測値とする。本手法は打ち切りデータを扱うため、Ｋ＝ｎとし、値が観測されていれば（ｗｉｊ＝１）、μｉｊ＝ｘｉｊ、そうでなければ（ｗｉｊ＝０）、μｉｊ＝ｘｉｊ＋εと設定してもよい。εは０以上の値を取る確率分布（例えば指数分布）からランダムに生成した値を表す。また、データ数が多い場合には、例えばランダムに選んだ１００個のデータのみを用いてもよいし、事前知識に基づいて設定したコンポーネントを用いてもよい。標準偏差σは交差検証法等により決定することが可能である。
　上記モデルを用いたときの打ち切りデータ
Figure JPOXMLDOC01-appb-I000014

の生成過程は次のように記述できる。まず、各サンプルについて観測終了までの時間の長さ
Figure JPOXMLDOC01-appb-I000015

が既知のもと、打ち切りが起こるか否かを表す変数
Figure JPOXMLDOC01-appb-I000016

が、下記式（２）の確率分布に従い生成される。
Figure JPOXMLDOC01-appb-M000017
　ただし、
Figure JPOXMLDOC01-appb-I000018

と、
Figure JPOXMLDOC01-appb-I000019

とのうち、ｗｉｊ＝１である、観測値が観測された要素の集合を
Figure JPOXMLDOC01-appb-I000020

とし、ｗｉｊ＝０である、観測値が観測されなかった集合を
Figure JPOXMLDOC01-appb-I000021

とする。また、観測値が観測された要素と、観測値が観測されなかった要素とを区別しない場合に、観測データと呼ぶこととする。コンポーネントの確率分布
Figure JPOXMLDOC01-appb-I000022

に、上記式（２）のガウス分布のように標準的な確率分布を利用すると、
Figure JPOXMLDOC01-appb-I000023

を累積密度関数を用いて解析的に計算できる。
Figure JPOXMLDOC01-appb-I000024

であり、少なくとも一つ観測された要素がある場合、
Figure JPOXMLDOC01-appb-I000025

は下記式（４）で表される分布にしたがって生成される。
Figure JPOXMLDOC01-appb-M000026
　ここで、
Figure JPOXMLDOC01-appb-I000027

はデルタ関数であり、ｆｔｒはｆを観測されなかった要素に関して周辺化した分布の切断分布であり、下記式で表される。
Figure JPOXMLDOC01-appb-M000028
　ただし、ｆｔｒは下記式（５）～（６）に従う。
Figure JPOXMLDOC01-appb-M000029
　コンポーネントが上記式（２）のガウス分布である場合、
Figure JPOXMLDOC01-appb-I000030

は
Figure JPOXMLDOC01-appb-I000031

となる。ただし、
Figure JPOXMLDOC01-appb-I000032

はそれぞれ
Figure JPOXMLDOC01-appb-I000033

から
Figure JPOXMLDOC01-appb-I000034

と
Figure JPOXMLDOC01-appb-I000035

に対応する次元の要素を抜き出したベクトルである。
Figure JPOXMLDOC01-appb-I000036

であり、一つも観測された要素がない場合、下記式（７）で示すように、デルタ関数のみの表現となる。
Figure JPOXMLDOC01-appb-M000037
　よって、上記をまとめると、各打ち切りデータの生成確率は、下記式（８）で与えられる。
Figure JPOXMLDOC01-appb-M000038
＜＜ＫＬダイバージェンスとＬ２ダイバージェンス＞＞
　次に、提案手法の目的関数を定義する際に利用するダイバージェンスについて記す。良く知られるように、確率分布ｐ（ｘ）とｑ（ｘ）とに対するカルバックライブラー（ＫＬ）ダイバージェンスは下記式（９）で定義される。
Figure JPOXMLDOC01-appb-M000039
　これに加え、本発明の実施の形態では、下記で定義されるＬ２ダイバージェンス（参考文献２）も利用する場合についても説明する（下記式（１０））。
Figure JPOXMLDOC01-appb-M000040
　Ｌ２ダイバージェンスは、２つの確率密度関数の２乗誤差として定義されている。何れのダイバージェンスを用いるべきかについては、問題に応じて異なる。よって、本発明の実施の形態では、何れのダイバージェンスを利用する場合にも適用できるように、２種類の手法を構築することとした。
＜＜ＫＬダイバージェンスの最適化による推定＞＞
　まず、目的関数としてＫＬダイバージェンスを利用する場合の提案手法を示す。打ち切りデータの確率分布を表すモデル
Figure JPOXMLDOC01-appb-I000041

と、打ち切りデータから得られる真の確率分布
Figure JPOXMLDOC01-appb-I000042

とのＫＬダイバージェンスは、上記式（９）の定義に従い、下記式（１１）で与えられる。
Figure JPOXMLDOC01-appb-M000043
　これは、下記式（１２）のように式変形することができる。
Figure JPOXMLDOC01-appb-M000044
　定数項を除去し、未知である真の分布
Figure JPOXMLDOC01-appb-I000045

に関する期待値を標本平均で置き換えれば、下記の式（１３）が導かれる。
Figure JPOXMLDOC01-appb-M000046
　ここで、
Figure JPOXMLDOC01-appb-I000047

は、データから計算できる量であり、これを目的関数とすることでアルゴリズムを導く。具体的には、下記式（１４）の最適化問題を解けばよい。
Figure JPOXMLDOC01-appb-M000048
　ここで、制約条件（パラメタの要素が０以上、かつ、和が１）は、混合モデルｆが確率分布となるためのものである。ラグランジュの未定乗数法を用いると、上記の最適化問題の解は、下記式（１５）を満たすことが分かる。
Figure JPOXMLDOC01-appb-M000049
　よって、上記式（１５）に基づいて
Figure JPOXMLDOC01-appb-I000050

の更新を繰り返すことにより、最適化が可能となる。なお、参考文献１の手法と同様に、計算量の削減と収束を早めるため、パラメタ更新の際にθｋがある閾値（例えば１０－３／ｎ）より小さい場合には、θｋ＝０と設定した後に、全体を和が１になるように調整する、再正規化操作を行ってもよい。
＜＜Ｌ２ダイバージェンスの最適化による推定＞＞
　次に、Ｌ２ダイバージェンスを利用する場合の提案手法を示す。ＫＬダイバージェンスとは異なり、Ｌ２ダイバージェンスの定義から直接目的関数を定義することはせず、ＫＬダイバージェンスを用いた際の目的関数に注目することにより新たな目的関数を定義する。
　ＫＬダイバージェンスを用いた際の目的関数（式（１２））に注目すると、混合モデルｆの周辺分布
Figure JPOXMLDOC01-appb-I000051

と、変数が与えられた下での観測値の真の分布
Figure JPOXMLDOC01-appb-I000052

のＫＬダイバージェンスに対応する項と、観測値が観測されなかったことを表す変数のモデルの分布
Figure JPOXMLDOC01-appb-I000053

と、真の分布
Figure JPOXMLDOC01-appb-I000054

を用いた対数尤度比に対応する項との２つの項をそれぞれ
Figure JPOXMLDOC01-appb-I000055

で重み付き和を取ることで構成されていることが分かる。
　この洞察に基づき、これら２つの各項でＫＬダイバージェンス／対数尤度比を用いている箇所をＬ２ダイバージェンスに置き換えることで、目的関数を下記式のように設計することができる。
Figure JPOXMLDOC01-appb-M000056
　これを式変形すれば、
Figure JPOXMLDOC01-appb-I000057

となる。更に、定数項（Ｃｏｎｓｔ）を除去し、真の分布に関する平均を標本平均で置き換えた下記式の最適化を考える。
Figure JPOXMLDOC01-appb-M000058
　ただし、
Figure JPOXMLDOC01-appb-I000059

　とおいた。これを目的関数とすることで、アルゴリズムを導出する。
Figure JPOXMLDOC01-appb-I000060

は、行列・ベクトル形式で、下記式（１６）、（１７）のように表すことができる。
Figure JPOXMLDOC01-appb-M000061
　ただし、下記式（１８）～（２０）に従う。
Figure JPOXMLDOC01-appb-M000062
　また、
Figure JPOXMLDOC01-appb-I000063

　は、コンポーネントの分布にガウス分布を用いる場合、解析的に計算できる値であり、下記式（２１）のように表すことができる。
Figure JPOXMLDOC01-appb-M000064
　上記式（２１）中で、
Figure JPOXMLDOC01-appb-I000065

に関して総和をとる際、
Figure JPOXMLDOC01-appb-I000066

や
Figure JPOXMLDOC01-appb-I000067

は
Figure JPOXMLDOC01-appb-I000068

の値によって異なる値であることを明記しておく。これにより、目的関数は
Figure JPOXMLDOC01-appb-I000069

に関する２次の形式で表現されることが分かる。よって、パラメタの推定値は下記で示す制約付きの２次最適化問題を解くことで得ることができる。
Figure JPOXMLDOC01-appb-M000070
　数値ソルバーを用いて直接上記式（２２）の最適化問題を解くことができる。その際、正則化項を加えた下記の最適化問題を解くことにしてもよい。
Figure JPOXMLDOC01-appb-M000071
　ただし、βがハイパーパラメタを表す。また、例えば次に示す近似的な方法を用いてもよい。上記式（２２）の最適化問題から制約条件を外した問題の最適解は下記式（２３）、（２４）のように求められる。
Figure JPOXMLDOC01-appb-M000072
　ただし、
Figure JPOXMLDOC01-appb-I000073

は正則化項、βはハイパーパラメタであり、パラメタの発散を防ぐ効果がある。
Figure JPOXMLDOC01-appb-I000074

は、値が０以上かつ和が１であるから、下記式（２５）の処理で条件を満たす
Figure JPOXMLDOC01-appb-I000075

が得られる。
Figure JPOXMLDOC01-appb-M000076
　ただし、
Figure JPOXMLDOC01-appb-I000077

は、
Figure JPOXMLDOC01-appb-I000078

ノルムを表す
　本発明の実施の形態では、以上説明した２つの手法の何れかにより、モデルのパラメタを推定することにより、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
＜本発明の実施の形態に係る推定装置の構成＞
　次に、図５及び図６を参照して、本発明の実施の形態に係る推定装置１の構成について説明する。図５は、本発明の実施の形態に係る推定装置１として機能するコンピュータの概略構成を示すブロック図である。図６は、本発明の実施の形態に係る推定装置１の構成を示すブロック図である。
　図５に示すように、推定装置１は、ＣＰＵ１１０と、ＲＡＭ等のメモリ１２０と、通信インターフェース（ＩＦ）部１３０と、キーボード等の入力部１４０と、ディスプレイ等の表示部１５０と、後述する推定処理ルーチンを実行するためのプログラム１７０を記憶したＲＯＭ等の記憶部１６０とを備えたコンピュータで構成されている。また、ＣＰＵ１１０、メモリ１２０、通信ＩＦ部１３０、入力部１４０、表示部１５０、及び記憶部１６０は、バス１００を介して接続されている。また、通信ＩＦ部１３０は、ＬＡＮケーブル等の通信回線により外部装置２と接続されている。なお、ネットワーク（図示しない）を介して外部装置２と接続される構成としてもよい。
　図６に示すように、本実施形態に係る推定装置１は、データ処理部１０と、パラメタ推定部２０と、パラメタ出力部３０と、記憶部４０と、入力部５０と、出力部６０とを備えて構成される。
　データ処理部１０は、入力部５０が受け付けた観測値が観測されたサンプルの観測データ
Figure JPOXMLDOC01-appb-I000079

と、観測値が観測されなかったサンプルの観測データ
Figure JPOXMLDOC01-appb-I000080

と、各サンプルについて観測値が観測されたか否かを表す変数
Figure JPOXMLDOC01-appb-I000081

とを含む打ち切りデータ
Figure JPOXMLDOC01-appb-I000082

をデータ記憶部４１に格納する。サンプルは、例えば、上述の図１、３の例では各機器のこと、図２、４の例では各利用者のことである。
　パラメタ推定部２０は、入力部５０が受け付けた打ち切りデータ
Figure JPOXMLDOC01-appb-I000083

の各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、打ち切りデータ
Figure JPOXMLDOC01-appb-I000084

の確率分布を表すモデル
Figure JPOXMLDOC01-appb-I000085

と、入力部５０が受け付けた打ち切りデータ
Figure JPOXMLDOC01-appb-I000086

から得られる、当該打ち切りデータ
Figure JPOXMLDOC01-appb-I000087

の真の確率分布
Figure JPOXMLDOC01-appb-I000088

とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタ
Figure JPOXMLDOC01-appb-I000089

を推定する。
　具体的には、パラメタ推定部２０は、まず、当該打ち切りデータ
Figure JPOXMLDOC01-appb-I000090

の真の確率分布
Figure JPOXMLDOC01-appb-I000091

を求める。
　次に、パラメタ推定部２０は、目的関数を、当該打ち切りデータ
Figure JPOXMLDOC01-appb-I000092

の確率分布を表すモデル
Figure JPOXMLDOC01-appb-I000093

と、当該打ち切りデータ
Figure JPOXMLDOC01-appb-I000094

の真の確率分布
Figure JPOXMLDOC01-appb-I000095

とのＫＬダイバージェンス、又はＬ２ダイバージェンスとして、パラメタを推定する。
　ここで、打ち切りデータの確率分布を表すモデル
Figure JPOXMLDOC01-appb-I000096

は、上記式（８）で示すように、観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、観測データの確率分布
Figure JPOXMLDOC01-appb-I000097

と、観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、変数が与えられた下での、変数の確率分布
Figure JPOXMLDOC01-appb-I000098

と、を用いて表される。
　パラメタ推定部２０は、ＫＬダイバージェンスを用いる場合、上記式（１５）のパラメタ更新を繰り返すことにより、パラメタを推定する。また、パラメタ推定部２０は、Ｌ２ダイバージェンスを用いる場合、上記式（２４）、（２５）を用いてパラメタを推定する。
　そして、パラメタ推定部２０は、推定したパラメタ
Figure JPOXMLDOC01-appb-I000099

をパラメタ記憶部４２に格納する。
　パラメタ出力部３０は、パラメタ記憶部４２のパラメタ
Figure JPOXMLDOC01-appb-I000100

を取得し、パラメタ出力部３０に当該パラメタ
Figure JPOXMLDOC01-appb-I000101

を渡す。
　記憶部４０は、データ記憶部４１と、パラメタ記憶部４２とを含む。データ記憶部４１には、打ち切りデータ
Figure JPOXMLDOC01-appb-I000102

が格納される。また、パラメタ記憶部４２には、モデルのパラメタ
Figure JPOXMLDOC01-appb-I000103

が格納される。
　入力部５０は、外部装置２から、打ち切りデータ
Figure JPOXMLDOC01-appb-I000104

の入力を受け付ける。そして、入力部５０は、受け付けた打ち切りデータ
Figure JPOXMLDOC01-appb-I000105

を、データ処理部１０に渡す。
　出力部６０は、パラメタ出力部３０から受け取ったモデルのパラメタ
Figure JPOXMLDOC01-appb-I000106

を、外部装置２へ出力する。
＜本発明の実施の形態に係る推定装置の作用＞
　図７は、本発明の実施の形態に係る推定処理ルーチンを示すフローチャートである。
　入力部５０に打ち切りデータ
Figure JPOXMLDOC01-appb-I000107

が入力されると、推定装置１において、図７に示す推定処理ルーチンが実行される。
　まず、ステップＳ１００において、入力部５０は、観測値が観測されたサンプルの観測データ
Figure JPOXMLDOC01-appb-I000108

と、観測値が観測されなかったサンプルの観測データ
Figure JPOXMLDOC01-appb-I000109

と、各サンプルについて観測値が観測されたか否かを表す変数
Figure JPOXMLDOC01-appb-I000110

とを含む打ち切りデータ
Figure JPOXMLDOC01-appb-I000111

の入力を受け付ける。また、入力部５０は、各サンプルｉについての観測終了時刻までの時間の長さの入力を受け付ける。
　ステップＳ１１０において、パラメタ推定部２０は、目的関数を、当該打ち切りデータ
Figure JPOXMLDOC01-appb-I000112

の確率分布を表すモデル
Figure JPOXMLDOC01-appb-I000113

と、当該打ち切りデータ
Figure JPOXMLDOC01-appb-I000114

の確率分布
Figure JPOXMLDOC01-appb-I000115

とのＫＬダイバージェンス、又はＬ２ダイバージェンスとして、上記式（１５）のパラメタ更新を繰り返すことにより、パラメタを推定するか、あるいは、上記式（２４）、（２５）を用いてパラメタを推定する。
　ステップＳ１２０において、出力部６０は、上記ステップＳ１１０により推定されたパラメタ
Figure JPOXMLDOC01-appb-I000116

を出力する。
　以上説明したように、本発明の実施形態に係る推定装置によれば、観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、観測データの確率密度関数を用いて表される、当該打ち切りデータの確率分布を表すモデルと、当該打ち切りデータから得られる、当該打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、モデルのパラメタを推定するため、計算時間を抑えて、打ち切りデータの確率分布を表すモデルのパラメタを精度良く推定することができる。
　なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
　例えば、上述の実施の形態では、ダイバージェンスとしてＫＬダイバージェンス又はＬ２ダイバージェンスを用いる場合について説明したが、これに限定されるものではなく、他のダイバージェンスを用いることも可能である。
　また、上述の実施の形態では、時系列データである打ち切りデータを前提に記載しているが、これに限定されるものではなく、時系列データではない任意の打ち切りデータに対しても本発明の適用が可能である。
　また、上述の実施の形態に係る推定装置１は、各部の処理をプログラムとして構築し、推定装置として利用されるコンピュータにインストールして実行させる構成として説明したが、ネットワークを介して流通させる構成としてもよい。
　また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。
１     推定装置
２     外部装置
１０   データ処理部
２０   パラメタ推定部
３０   パラメタ出力部
４０   記憶部
４１   データ記憶部
４２   パラメタ記憶部
５０   入力部
６０   出力部
１００ バス
１１０ ＣＰＵ
１２０ メモリ
１３０ 通信ＩＦ部
１４０ 入力部
１５０ 表示部
１６０ 記憶部
１７０ プログラム

Claims (5)

  1. 　観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定装置であって、
    　前記打ち切りデータの入力を受け付ける入力部と、
    　前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定するパラメタ推定部と、
    　を含む推定装置。
  2. 　前記打ち切りデータの確率分布を表すモデルは、
    　前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数の確率分布と、
    　前記観測データの確率密度関数と、各サンプルについて予め与えられた、観測終了までの時間の長さとを用いて表される、前記変数が与えられた下での、前記観測データの確率分布と、
    　を用いて表される請求項１記載の推定装置。
  3. 　前記目的関数は、前記打ち切りデータの確率分布を表すモデルと、前記打ち切りデータの確率分布とのカルバックライブラーダイバージェンス、又はＬ２ダイバージェンスである
    　請求項１又は２記載の推定装置。
  4. 　観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する推定方法であって、
    　入力部が、前記打ち切りデータの入力を受け付け、
    　パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する
    　推定方法。
  5. 　観測値が観測されたサンプルの観測データと、観測値が観測されなかったサンプルの観測データと、各サンプルについて観測値が観測されたか否かを表す変数とを含む打ち切りデータの確率分布を表すモデルのパラメタを推定する処理をコンピュータに実行させるプログラムであって、
    　入力部が、前記打ち切りデータの入力を受け付け、
    　パラメタ推定部が、前記入力部が受け付けた前記打ち切りデータの各サンプルに対応する、観測値の分布を表す各コンポーネントの混合モデルで表される、前記観測データの確率密度関数を用いて表される、前記打ち切りデータの確率分布を表すモデルと、前記入力部が受け付けた前記打ち切りデータから得られる、前記打ち切りデータの確率分布とのダイバージェンスである目的関数を最適化することにより、前記モデルのパラメタを推定する
    　ことを含む処理をコンピュータに実行させるプログラム。
PCT/JP2020/004909 2019-02-21 2020-02-07 推定装置、推定方法、及びプログラム Ceased WO2020170867A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US17/431,791 US20220138375A1 (en) 2019-02-21 2020-02-07 Estimation device, estimation method, and program

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2019-029769 2019-02-21
JP2019029769A JP7107246B2 (ja) 2019-02-21 2019-02-21 推定装置、推定方法、及びプログラム

Publications (1)

Publication Number Publication Date
WO2020170867A1 true WO2020170867A1 (ja) 2020-08-27

Family

ID=72144974

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/004909 Ceased WO2020170867A1 (ja) 2019-02-21 2020-02-07 推定装置、推定方法、及びプログラム

Country Status (3)

Country Link
US (1) US20220138375A1 (ja)
JP (1) JP7107246B2 (ja)
WO (1) WO2020170867A1 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7107246B2 (ja) * 2019-02-21 2022-07-27 日本電信電話株式会社 推定装置、推定方法、及びプログラム
US11921488B2 (en) * 2020-12-15 2024-03-05 Xerox Corporation System and method for machine-learning-enabled micro-object density distribution control with the aid of a digital computer
US20230359882A1 (en) * 2022-05-06 2023-11-09 International Business Machines Corporation Training a neural network to achieve average calibration
CN119128368B (zh) * 2024-11-14 2025-02-21 中国建设银行股份有限公司 回归模型的构建方法及装置、程序产品、存储介质

Family Cites Families (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040117051A1 (en) * 2000-06-06 2004-06-17 Ford Dean M. Method of determining a cumulative distribution function confidence bound
US8438126B2 (en) * 2007-07-13 2013-05-07 The Regents Of The University Of California Targeted maximum likelihood estimation
US8204714B2 (en) * 2010-02-26 2012-06-19 International Business Machines Corporation Method and computer program product for finding statistical bounds, corresponding parameter corners, and a probability density function of a performance target for a circuit
US20140149205A1 (en) * 2012-11-29 2014-05-29 Adobe Systems Incorporated Method and Apparatus for an Online Advertising Predictive Model with Censored Data
JP6194450B2 (ja) * 2013-04-15 2017-09-13 株式会社メガチップス 状態推定装置、プログラムおよび集積回路
JP6366999B2 (ja) * 2014-05-22 2018-08-01 株式会社メガチップス 状態推定装置、プログラムおよび集積回路
PT3188060T (pt) * 2014-08-27 2023-06-27 Nec Corp Dispositivo de simulação, método de simulação e meios de memória
WO2017170086A1 (ja) * 2016-03-31 2017-10-05 日本電気株式会社 情報処理システム、情報処理装置、シミュレーション方法およびシミュレーションプログラムが記録された記録媒体
US20190057065A1 (en) * 2016-04-01 2019-02-21 Hitachi, Ltd. Modeling device
JP6686151B2 (ja) * 2017-01-27 2020-04-22 三菱日立パワーシステムズ株式会社 モデルパラメータ値推定装置及び推定方法、プログラム、プログラムを記録した記録媒体、モデルパラメータ値推定システム
US10528644B1 (en) * 2017-06-30 2020-01-07 Cadence Design Systems, Inc. Estimation and visualization of a full probability distribution for circuit performance obtained with Monte Carlo simulations over scaled sigma sampling
JP6935765B2 (ja) * 2018-02-13 2021-09-15 日本電信電話株式会社 動的分布推定装置、方法、及びプログラム
JP7064356B2 (ja) * 2018-03-14 2022-05-10 株式会社日立製作所 将来状態推定装置および将来状態推定方法
WO2019217876A1 (en) * 2018-05-10 2019-11-14 Equifax Inc. Training or using sets of explainable machine-learning modeling algorithms for predicting timing of events
JP7172616B2 (ja) * 2019-01-11 2022-11-16 日本電信電話株式会社 データ解析装置、方法、及びプログラム
US10977341B2 (en) * 2019-01-18 2021-04-13 Hide Housing Corporation Insolation probability distribution analysis method, insolation probability distribution analysis system, insolation probability distribution analysis program product, insolation normalization statistical analysis method, insolation normalization statistical analysis system, and insolation normalization statistical analysis program product
JP7107246B2 (ja) * 2019-02-21 2022-07-27 日本電信電話株式会社 推定装置、推定方法、及びプログラム
WO2020261447A1 (ja) * 2019-06-26 2020-12-30 日本電信電話株式会社 パラメタ推定装置、パラメタ推定方法、及びパラメタ推定プログラム
US11010222B2 (en) * 2019-08-29 2021-05-18 Sap Se Failure mode specific analytics using parametric models
US20220343199A1 (en) * 2019-10-02 2022-10-27 Nippon Telegraph And Telephone Corporation Parameter estimation device, parameter estimation method, and parameter estimation program
JP7374868B2 (ja) * 2020-08-28 2023-11-07 株式会社東芝 情報処理装置、情報処理方法およびプログラム

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
BASU, SRABASHI ET AL.: "Robust and efficient parametric estimation for censored survival data", ANNALS OF THE INSTITUTE OF STATISTICAL MATHEMATICS [ ONLINE, vol. 58, no. 2, 2006, pages 341 - 355, XP019390068, ISSN: 0020-3157, Retrieved from the Internet <URL:https://www.ism.ac.jp/editsec/aism/pdf/058_2_0341.pdf> [retrieved on 20200318], DOI: 10.1007/s10463-005-0004-x *
KOHJIMA, MASAHIRO; KIYOTAKE, HIROSHI; MATSUBAYASHI, TATSUSHI; SHIOHARA, HISAKO; TODA, HIROYUKI: "Introduction of Online EM Method for Mixed Models for Censored Data, Arrival Time Estimation for Large-Scale Gathering Events. DEIM Forum 2018", 16TH ANNUAL MEETING OF THE DATABASE SOCIETY OF JAPAN, DEIM FORUM 2018, March 2018 (2018-03-01), pages 8 pp., XP009523480 *

Also Published As

Publication number Publication date
JP7107246B2 (ja) 2022-07-27
JP2020135554A (ja) 2020-08-31
US20220138375A1 (en) 2022-05-05

Similar Documents

Publication Publication Date Title
WO2020170867A1 (ja) 推定装置、推定方法、及びプログラム
Taylor et al. Performance estimation toolbox (PESTO): Automated worst-case analysis of first-order optimization methods
Eduarda Da Silva et al. Difference equations for the higher‐order moments and cumulants of the INAR (1) model
Ma Adaptive shrinkage in Pólya tree type models
Pan et al. Adaptive importance sampling for extreme quantile estimation with stochastic black box computer models
CN115828414A (zh) 雷达罩结构的分布参数不确定性可靠性灵敏度分析方法
Mohammadi et al. An Introduction to the BDgraph for Bayesian Graphical Models
Antoniano‐Villalobos et al. A nonparametric model for stationary time series
Cui et al. Multilevel dimension-independent likelihood-informed MCMC for large-scale inverse problems
WO2020145252A1 (ja) データ解析装置、方法、及びプログラム
Song et al. Wasserstein generative regression
Vexler et al. Computing critical values of exact tests by incorporating Monte Carlo simulations combined with statistical tables
de la Fuente et al. An efficient nonlinear programming strategy for PCA models with incomplete data sets
JP2018067227A (ja) データ分析装置、データ分析方法、データ分析処理プログラム
Lu et al. Likelihood based confidence intervals for the tail index
US11443136B2 (en) System and method for signal pre-processing based on data driven models and data dependent model transformation
Wen et al. Laplace Error Penalty‐based Variable Selection in High Dimension
Liu et al. On Random Batch Methods (RBM) for interacting particle systems driven by Lévy processes
JP6398991B2 (ja) モデル推定装置、方法およびプログラム
Wang et al. Constrained spline regression in the presence of AR (p) errors
Alanazi A mixture of regular vines for multiple dependencies
Sköld et al. Density estimation for the Metropolis–Hastings algorithm
Chesi Stabilization and entropy reduction via SDP-based design of fixed-order output feedback controllers and tuning parameters
Delattre et al. Local distortion and μ-mass of the cells of one dimensional asymptotically optimal quantizers
Chan et al. Model predictive control of Hammerstein systems with multivariable nonlinearities

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

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20758795

Country of ref document: EP

Kind code of ref document: A1