(57)【要約】【課題】精度よく、対象物の生存確率を推定できる。【解決手段】初期化処理部３０が、対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、イベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化する。発生時刻推定部３２が、イベントの発生が観測されたｉ番目のデータの各々に関して、イベント未発生観測経過時間Ｔｉ，Ｎと、最終観測経過時間Ｔｉと、初期化又は更新されたパラメータとに基づいて、イベント発生経過時間Ｔｉ，Ｒを推定する。パラメータ推定部３４が、ｉ番目のデータの各々についてのイベント発生経過時間Ｔｉ，Ｒ、最終観測経過時間Ｔｉ、及びイベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、パラメータを推定する。【選択図】図４
2018005580.tif 000002
【特許請求の範囲】
【請求項１】
イベントが発生するまで生存する対象物の生存確率を推定する生存確率推定装置であって、
前記対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、観測が開始されてからイベントが実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化する初期化処理部と、
前記イベントの発生が観測されたｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎと、観測が開始されてからイベントが発生したか否かが最後に観測された時刻までの最終観測経過時間Ｔｉと、初期化又は更新された前記パラメータとに基づいて、前記イベント発生経過時間Ｔｉ，Ｒを推定する発生時刻推定部と、
ｉ番目のデータの各々についての前記イベント発生経過時間Ｔｉ，Ｒ、前記最終観測経過時間Ｔｉ、及び前記イベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、前記パラメータを推定するパラメータ推定部と、
前記発生時刻推定部による推定、及び前記パラメータ推定部による推定を、予め定めた条件を満たすまで繰り返す収束判定処理部と、
を含む生存確率推定装置。
【請求項２】
前記パラメータ推定部は、前記モデルとして、指数分布のモデルを用いる請求項１に記載の生存確率推定装置。
【請求項３】
前記発生時刻推定部により推定された前記イベント発生経過時間Ｔｉ，Ｒを出力する出力部を更に含む請求項１又は請求項２に記載の生存確率推定装置。
【請求項４】
前記発生時刻推定部は、ｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までの経過時間Ｔｉ，Ｎにおける、イベントが発生していない確率Ｓ（Ｔｉ，Ｎ）と、前記最終観測経過時間Ｔｉにおける、イベントが発生していない確率Ｓ（Ｔｉ）とを、初期化又は更新された前記パラメータを用いた前記モデルに基づいて計算し、前記確率Ｓ（Ｔｉ，Ｎ）から前記確率Ｓ（Ｔｉ）までの範囲から定まる確率に基づいて、前記イベント発生経過時間Ｔｉ，Ｒを推定する請求項１〜請求項３のいずれか１項に記載の生存確率推定装置。
【請求項５】
前記パラメータ推定部は、前記イベントの発生が観測されたｉ番目のデータの各々に関する、前記イベント発生経過時間Ｔｉ，Ｒにおける、イベントが発生する確率ｆ（Ｔｉ，Ｒ）と、前記イベントの発生が観測されていないｉ番目のデータの各々に関する前記確率Ｓ（Ｔｉ）とを用いて表される前記尤度式を最適化するように、前記パラメータを推定する請求項１〜請求項４のいずれか１項に記載の生存確率推定装置。
【請求項６】
イベントが発生するまで生存する対象物の生存確率を推定する生存確率推定装置における生存確率推定方法であって、
初期化処理部が、前記対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、観測が開始されてからイベントが実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化するステップと、
発生時刻推定部が、前記イベントの発生が観測されたｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎと、観測が開始されてからイベントが発生したか否かが最後に観測された時刻までの最終観測経過時間Ｔｉと、初期化又は更新された前記パラメータとに基づいて、前記イベント発生経過時間Ｔｉ，Ｒを推定するステップと、
パラメータ推定部が、ｉ番目のデータの各々についての前記イベント発生経過時間Ｔｉ，Ｒ、前記最終観測経過時間Ｔｉ、及び前記イベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、前記パラメータを推定するステップと、
収束判定処理部が、前記発生時刻推定部による推定、及び前記パラメータ推定部による推定を、予め定めた条件を満たすまで繰り返すステップと、
を含む生存確率推定方法。
【請求項７】
出力部が、前記発生時刻推定部により推定された前記イベント発生経過時間Ｔｉ，Ｒを出力するステップを更に含む請求項６に記載の生存確率推定方法。
【請求項８】
コンピュータを、請求項１〜請求項５のいずれか１項に記載の生存確率推定装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、生存確率推定装置、方法、及びプログラムに係り、特に、イベントが発生するまで生存する対象物の生存確率を推定するための生存確率推定装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
機械の耐久年数や、生物の生存期間を推定する代表的な技術として生存時間分析というものが広く利用されている。生存時間分析では経過時間ｔ（例えば生物が誕生してからの生存期間（年齢））において、該当データにイベントが発生する（例えば生物の死亡）確率を算出する。生存時間分析の既存技術としては、１．カプランマイヤー、２．Ｃｏｘ比例ハザードモデル、３．パラメトリックモデルが存在する。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】David M Diez.:Survival Analysis in R,インターネット<URL:http://www.openintro.org/stat/down/Survival-Analysis-in-R.pdf>,(検索日時:2016/02/01),2013.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
既存技術のうち、１．カプランマイヤーは入力されたデータを集計処理することで生存確率を算出し、生存確率を表す関数のパラメータ推定を行わないのに対し、２．Ｃｏｘ比例ハザードモデル、３．パラメトリックモデルではパラメータ推定を行う違いがある。また、いずれの技術も入力データとしては「イベントの発生有無」（例えば生物の死亡が確認されたかどうか）、および「該当データの観測開始から、イベントの発生が確認されるまでの経過時間（イベントが発生していない場合は該当データの観測期間）」の情報を元に生存確率の推定を行う。しかしながら既存技術は共通してある問題を抱えている。それは既存技術の入力データとして用いている「経過時間」というのは、「イベントが実際に発生するまでの経過時間」ではなく「イベントの発生が確認されるまでの経過時間」であるため、両者が大きくずれていた場合、生存確率の推定が不正確になる（図１参照）。
【０００５】
このようなケースは往々にして発生する。例えば非常用電源は通常通電していないため、劣化故障などが見つけにくく、メンテナンス時に故障が発見されるケースがよくある。またビルのメンテナンスでも、非破壊検査が行われる頻度は非常に低いため、実際にビルに破損が発生した時期を推定することは困難である。そのため非常用電源を例にとると、潜在的な故障が発見されにくく、対象の本質的な耐用年数の推定が難しいため、適切なメンテナンス時期がわからず、実際に利用しようという段階において故障している事態が発生しやすい。
【０００６】
本発明は、上記課題を鑑みて成されたものであり、精度よく、対象物の生存確率を推定できる生存確率推定装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記目的を達成するために、第１の発明に係る生存確率推定装置は、イベントが発生するまで生存する対象物の生存確率を推定する生存確率推定装置であって、前記対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、観測が開始されてからイベントが実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化する初期化処理部と、前記イベントの発生が観測されたｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎと、観測が開始されてからイベントが発生したか否かが最後に観測された時刻までの最終観測経過時間Ｔｉと、初期化又は更新された前記パラメータとに基づいて、前記イベント発生経過時間Ｔｉ，Ｒを推定する発生時刻推定部と、ｉ番目のデータの各々についての前記イベント発生経過時間Ｔｉ，Ｒ、前記最終観測経過時間Ｔｉ、及び前記イベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、前記パラメータを推定するパラメータ推定部と、前記発生時刻推定部による推定、及び前記パラメータ推定部による推定を、予め定めた条件を満たすまで繰り返す収束判定処理部と、を含んで構成されている。
【０００８】
また、第１の発明に係る生存確率推定装置において、前記パラメータ推定部は、前記モデルとして、指数分布のモデルを用いるようにしてもよい。
【０００９】
また、第１の発明に係る生存確率推定装置において、前記発生時刻推定部により推定された前記イベント発生経過時間Ｔｉ，Ｒを出力する出力部を更に含むようにしてもよい。
【００１０】
また、第１の発明に係る生存確率推定装置において、前記発生時刻推定部は、ｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までの経過時間Ｔｉ，Ｎにおける、イベントが発生していない確率Ｓ（Ｔｉ，Ｎ）と、前記最終観測経過時間Ｔｉにおける、イベントが発生していない確率Ｓ（Ｔｉ）とを、初期化又は更新された前記パラメータを用いた前記モデルに基づいて計算し、前記確率Ｓ（Ｔｉ，Ｎ）から前記確率Ｓ（Ｔｉ）までの範囲から定まる確率に基づいて、前記イベント発生経過時間Ｔｉ，Ｒを推定するようにしてもよい。
【００１１】
また、第１の発明に係る生存確率推定装置において、前記パラメータ推定部は、前記イベントの発生が観測されたｉ番目のデータの各々に関する、前記イベント発生経過時間Ｔｉ，Ｒにおける、イベントが発生する確率ｆ（Ｔｉ，Ｒ）と、前記イベントの発生が観測されていないｉ番目のデータの各々に関する前記確率Ｓ（Ｔｉ）とを用いて表される前記尤度式を最適化するように、前記パラメータを推定するようにしてもよい。
【００１２】
第２の発明に係る生存確率推定方法は、イベントが発生するまで生存する対象物の生存確率を推定する生存確率推定装置における生存確率推定方法であって、初期化処理部が、前記対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、観測が開始されてからイベントが実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化するステップと、発生時刻推定部が、前記イベントの発生が観測されたｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎと、観測が開始されてからイベントが発生したか否かが最後に観測された時刻までの最終観測経過時間Ｔｉと、初期化又は更新された前記パラメータとに基づいて、前記イベント発生経過時間Ｔｉ，Ｒを推定するステップと、パラメータ推定部が、ｉ番目のデータの各々についての前記イベント発生経過時間Ｔｉ，Ｒ、前記最終観測経過時間Ｔｉ、及び前記イベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、前記パラメータを推定するステップと、収束判定処理部が、前記発生時刻推定部による推定、及び前記パラメータ推定部による推定を、予め定めた条件を満たすまで繰り返すステップと、を含んで実行することを特徴とする。
【００１３】
また、第２の発明に係る生存確率推定方法において、出力部が、前記発生時刻推定部により推定された前記イベント発生経過時間Ｔｉ，Ｒを出力するステップを更に含むようにしてもよい。
【００１４】
第３の発明に係るプログラムは、コンピュータを、第１の発明に係る生存確率推定装置の各部として機能させるためのプログラムである。
【発明の効果】
【００１５】
本発明の生存確率推定装置、方法、及びプログラムによれば、対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、観測が開始されてからイベントが実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化し、イベントの発生が観測されたｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎと、観測が開始されてからイベントが発生したか否かが最後に観測された時刻までの最終観測経過時間Ｔｉと、初期化又は更新されたパラメータとに基づいて、イベント発生経過時間Ｔｉ，Ｒを推定し、ｉ番目のデータの各々についてのイベント発生経過時間Ｔｉ，Ｒ、最終観測経過時間Ｔｉ、及びイベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、パラメータを推定することを、予め定めた条件を満たすまで繰り返すことにより、精度よく、対象物の生存確率を推定できる、という効果が得られる。
【図面の簡単な説明】
【００１６】
【図１】既存技術と本発明の実施の形態の開発技術による生存確率の推定の一例を示す図である。
【図２】本発明の実施の形態に係る生存確率推定装置の機械のデータの入力例を示す図である。
【図３】本発明の実施の形態に係る生存確率推定装置による出力例を示す図である。
【図４】本発明の実施の形態に係る生存確率推定装置の構成を示すブロック図である。
【図５】故障有と故障無における経過時間Ｔｉ，Ｎ、Ｔｉ，Ｒ及び、Ｔｉの関係を表す図である。
【図６】経過時間と生存確率の関係を表す図である。
【図７】本発明の実施の形態に係る生存確率推定装置における生存確率推定処理ルーチンを示すフローチャートである。
【図８】本発明の実施の形態のシミュレーション結果の一例を示す図である。
【発明を実施するための形態】
【００１７】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１８】
＜本発明の実施の形態に係る概要＞
【００１９】
まず、本発明の実施の形態における概要を説明する。
【００２０】
本発明の実施の形態に係る開発技術はパラメトリックモデルのみに適用可能であり、カプランマイヤー、Ｃｏｘ比例ハザードモデルには適用できない。この手法は故障発生時刻が不明なデータについての、故障発生時刻の予測も行う。
【００２１】
以下、対象物を「機械」、イベントを「機械の故障」、イベントの発生が「指数分布」に従う場合を例として説明する（他には対象データを「生物」、イベントを「生物の死亡」、イベントの発生が「ワイブル分布」に従うなどの例が考えられる）。
【００２２】
本実施の形態の開発技術における、生存確率推定装置の入力と出力のイメージは図２、及び図３のようになる。図２の入力において着目してほしい点としては、従来技術では利用していなかった「故障無」が確認された時刻の情報を用いている点である。例えば非常用電源の場合、頻度は少ないが「試運転」を行うことや、周辺に設置したセンサ情報を確認することによって故障が発生していない時期の情報が得られることがある。開発技術ではこれらの「故障無」情報を利用することで、実際に故障が発生した時刻を絞り込み、より正確な故障発生時期の予測に役立てる。以下、本実施の形態の詳細を説明する。
【００２３】
＜本発明の実施の形態に係る生存確率推定装置の構成＞
【００２４】
次に、本発明の実施の形態に係る生存確率推定装置の構成について説明する。図４に示すように、本発明の実施の形態に係る生存確率推定装置１００は、ＣＰＵと、ＲＡＭと、後述する生存確率推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この生存確率推定装置１００は、機能的には図４に示すように入力部１０と、演算部２０と、出力部５０とを備えている。
【００２５】
入力部１０は、Ｎ個の機械の各々についての運転開始以降のデータに関して、以下の情報を受け付ける。
【００２６】
入力部１０は、ｉ番目の機械のデータに関して、故障したか否かを表すｃｉ（ｉ＝１，２，．．．，Ｍ）を受け付ける。ｃｉ＝１は故障発生観測、ｃｉ＝０は故障発生未観測を表す。
【００２７】
また、ｉ番目の機械のデータに関して、機械が設置（製造）されてから故障か否かの情報が最後に観測された時刻までの最終観測経過時間Ｔｉを受け付ける。
【００２８】
また、ｉ番目の機械のデータに関して、故障発生が観測された時（ｃｉ＝１）、機械が設置（製造）されてから故障していないことが確認された最後の時刻までのイベント未発生経過時間Ｔｉ，Ｎを受け付ける。
【００２９】
演算部２０は、初期化処理部３０と、発生時刻推定部３２と、パラメータ推定部３４と、収束判定部３６とを含んで構成されている。
【００３０】
初期化処理部３０は、以下に説明するように、ｉ番目の機械のデータに関して、最終観測経過時間Ｔｉによって、故障発生が観測された時（ｃｉ＝１）、機械が設置（製造）されてから故障が実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化する。また、初期化処理部３０は、経過時間ｔにおいて故障が発生する確率ｆ（ｔ）及び経過時間ｔまでに故障が発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化する。
【００３１】
ここで、パラメータの種類は、選択するモデルによるが、本実施の形態では、モデルとして、指数分布のモデルを用いる。指数分布のモデルの場合は、パラメータはλとする。また、最終観測経過時間Ｔｉ、イベント未発生観測経過時間Ｔｉ，Ｎ、イベント発生経過時間Ｔｉ，Ｒは、それぞれｉ番目の機械のデータが観測されてからの経過時間を表す（Ｔｉ，Ｎ＜Ｔｉ，Ｒ≦Ｔｉ）。図５は故障有と故障無におけるＴｉ，Ｎ、Ｔｉ，Ｒ及び、Ｔｉの関係を表したものである。Ｔｉ、Ｔｉ，Ｎは入力データから一意に定まるのに対して、Ｔｉ，Ｒは本実施の形態の開発技術によって推定を行う。初期化処理部３０においては、以下（１）式に従ってイベント発生経過時間Ｔｉ，Ｒを初期化する。
【００３２】
【数１】
2018005580.tif 000003
【００３３】
関数である、経過時間ｔにおいて故障が発生する確率ｆ（ｔ）と経過時間ｔまでに故障が発生していない確率Ｓ（ｔ）は、選択するモデルとモデルのパラメータによって定まり、以下（２）式の関係を持つ。
【００３４】
【数２】
2018005580.tif 000004
【００３５】
モデルによってｆ（ｔ）、及びＳ（ｔ）は変わりうるが、片方が定まれば他方は一意に定まる関係を持つ。例えば指数分布の場合、パラメータλによって以下（３）式、及び（４）式のように表される。
【００３６】
【数３】
2018005580.tif 000005
【００３７】
初期化処理部３０ではパラメータλの初期値を設定する。初期値に関する制約は基本的にはないが、例えば指数分布のモデルを選択した場合、パラメータλは以下（５）式で算出するなどが考えられる。
【００３８】
【数４】
2018005580.tif 000006
【００３９】
発生時刻推定部３２は、故障の発生が観測されたｉ番目の機械のデータの各々に関して、イベント未発生観測経過時間Ｔｉ，Ｎと、最終観測経過時間Ｔｉと、初期化又は更新されたパラメータλとに基づいて、イベント発生経過時間Ｔｉ，Ｒを推定する。
【００４０】
発生時刻推定部３２では、すべてのデータについて推定を行うのではなく、故障発生が観測されたデータ（ｃｉ＝１）のみに対しＴｉ，Ｒの推定を行う。故障発生が観測されていないデータ（ｃｉ＝０）については初期化処理に設定した値を保持する。まず、ｉ番目の機械のデータの各々に関して、観測が開始されてから故障が発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎにおける、故障が発生していない確率Ｓ（Ｔｉ，Ｎ）と、最終観測経過時間Ｔｉにおける、故障が発生していない確率Ｓ（Ｔｉ）とを、初期化又は更新されたパラメータλを用いたモデルに基づいて計算する。次に、確率Ｓ（Ｔｉ，Ｎ）から確率Ｓ（Ｔｉ）までの範囲から、一様乱数を用いて確率Ｓ（Ｔｉ，Ｒ）を求め、逆関数を用いることにより、イベント発生経過時間Ｔｉ，Ｒを推定する。図６は経過時間のＴｉ、Ｔｉ，Ｎ及び、Ｔｉ，Ｒと、確率のＳ（Ｔｉ）、Ｓ（Ｔｉ，Ｎ）、及びＳ（Ｔｉ，Ｒ）との関係を表したものである。
【００４１】
パラメータ推定部３４は、故障の発生が観測されたｉ番目の機械のデータの各々に関する、イベント発生経過時間Ｔｉ，Ｒにおける、故障が発生する確率ｆ（Ｔｉ，Ｒ）と、故障の発生が観測されていないｉ番目の機械のデータの各々に関する確率Ｓ（Ｔｉ）とを用いて表される以下（６）式の尤度式を最適化するように、パラメータλを推定する。
【００４２】
【数５】
2018005580.tif 000007
【００４３】
例えば指数分布のモデルを選択した場合は初期化処理部３０において説明したように、パラメータλとなり、以下（７）式によって、上記（６）式を最適化するパラメータλを算出することができる。
【００４４】
【数６】
2018005580.tif 000008
【００４５】
ただし、￣ｃｉは、ｃｉの反転を表す。
【００４６】
【数７】
2018005580.tif 000009
【００４７】
収束判定部３６は、発生時刻推定部３２による推定、及びパラメータ推定部３４による推定を、予め定めた条件を満たすまで繰り返す。収束判定処理としては様々なやり方が考えられる。例えば、条件を、一定回数（例えば１００回）、発生時刻推定部３２による推定、及びパラメータ推定部３４による推定を繰り返すこととし、繰り返し終了後に、出力部５０に移行する。推定対象であるパラメータの変化が一定値以下になった時、出力に移行する。
【００４８】
出力部５０は、発生時刻推定部３２により推定された経過時間ｔのイベント発生経過時間Ｔｉ，Ｒを出力する。また、経過時間ｔにおける生存確率Ｓ（ｔ）の情報（グラフや式）が主な出力対象となる。また、パラメータ（指数分布の場合はλ）等も出力可能である。
【００４９】
＜本発明の実施の形態に係る生存確率推定装置の作用＞
【００５０】
次に、本発明の実施の形態に係る生存確率推定装置１００の作用について説明する。入力部１０において、Ｎ個の機械の各々について、故障したか否かを表すｃｉ、機械が設置（製造）されてから故障か否かの情報が最後に観測された時刻までの最終観測経過時間Ｔｉ、及び故障発生が観測された時（ｃｉ＝１）、機械が設置（製造）されてから故障していないことが確認された最後の時刻までのイベント未発生経過時間Ｔｉ，Ｎを受け付けると、生存確率推定装置１００は、図７に示す生存確率推定処理ルーチンを実行する。
【００５１】
まず、ステップＳ１００では、ｉ番目の機械のデータ（ｉ＝１，・・・，Ｎ）の各々に関して、最終観測経過時間Ｔｉによって、故障発生が観測された時（ｃｉ＝１）、機械が設置（製造）されてから故障が実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを上記（１）式に従って初期化する。また、経過時間ｔにおいて故障が発生する確率ｆ（ｔ）及び経過時間ｔまでに故障が発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータλを上記（５）式に従って初期化する。
【００５２】
次に、ステップＳ１０２では、ｉ番目の機械のデータの各々に関して、観測が開始されてから故障が発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎにおける、故障が発生していない確率Ｓ（Ｔｉ，Ｎ）と、最終観測経過時間Ｔｉにおける、故障が発生していない確率Ｓ（Ｔｉ）とを、初期化又は更新されたパラメータλを用いたモデルに基づいて上記（３）式、及び（４）式に従って計算する。
【００５３】
ステップＳ１０４では、ｉ番目の機械のデータの各々に関して、ステップＳ１０２で計算した確率Ｓ（Ｔｉ，Ｎ）から確率Ｓ（Ｔｉ）までの範囲から、一様乱数を用いて確率Ｓ（Ｔｉ，Ｒ）を求め、逆関数を用いることにより、イベント発生経過時間Ｔｉ，Ｒを推定する。
【００５４】
ステップＳ１０６では、故障の発生が観測されたｉ番目の機械のデータの各々に関する、イベント発生経過時間Ｔｉ，Ｒにおける、故障が発生する確率ｆ（Ｔｉ，Ｒ）と、故障の発生が観測されていないｉ番目の機械のデータの各々に関する確率Ｓ（Ｔｉ）とを用いて表される上記（６）式の尤度式を最適化するように、上記（７）式に従って、パラメータλを推定する。
【００５５】
ステップＳ１０８では、収束判定に関する条件を満たしているか否かを判定し、条件を満たしていると判定された場合にはステップＳ１１０へ移行し、条件を満たしていないと判定された場合にはステップＳ１０４に戻って処理を繰り返す。
【００５６】
ステップＳ１１０では、ステップＳ１０６で推定された経過時間ｔのイベント発生経過時間Ｔｉ，Ｒを出力し、処理を終了する。
【００５７】
＜シミュレーション結果＞
【００５８】
本実施の形態の効果を検証するため以下のシミュレーションを行った。本実施の形態に係る開発技術、及び既存技術を用いて以下の人工データの生存確率推定を行う。
【００５９】
データは以下の（ｃｉ，Ｔｉ，Ｎ，Ｔｉ，Ｒ，Ｔｉ）で構成されている。
・故障発生が確認されているか否か（ｃｉ）
・機械が設置（製造）されてから故障していないことが確認された最後の時刻までのイベント未発生経過時間（Ｔｉ，Ｎ）
・機械が設置（製造）されてから故障が実際に発生した時刻までのイベント発生経過時間（Ｔｉ，Ｒ）
・機械が設置（製造）されてから故障か否かの情報が最後に観測された時刻までの最終観測経過時間（Ｔｉ）
【００６０】
ただしＴｉ，Ｒは正解データであって、通常の分析の際はわからない。
【００６１】
本シミュレーションでは、人工データとして以下を用いる。
【００６２】
・（１，７，９，１０）
・（１，３，５，７）
・（１，６，７，８）
・（１，２，３，５）
・（１，１０，１２，１３）
・（０，１３，−，１５）
・（０，４，−，８）
【００６３】
パラメトリックモデルの一つである指数分布（分布パラメータ）をベースに、既存技術と本実施の形態の開発技術を適用するものとする。開発技術の収束判定処理における繰返し回数は１０回とする。また、各技術の精度の検証のため正解データであるＴｉ，Ｒを用いた場合の推定結果と比較を行う。正解データがすべて判明している場合、既存技術と開発技術は同じ結果となる。
【００６４】
この時、既存技術を用いて推定された分布パラメータは
【００６５】
【数８】
2018005580.tif 000010
【００６６】
となった。
【００６７】
また、本実施の形態の開発技術を用いて推定された分布パラメータは
【００６８】
【数９】
2018005580.tif 000011
【００６９】
となった。また、正解データを用いた場合における推定された分布パラメータは
【００７０】
【数１０】
2018005580.tif 000012
【００７１】
となった。
【００７２】
分布パラメータλの各推定結果を元に、横軸：経過時間、縦軸：生存確率としたのが図８である。図８より本実施の形態の開発技術の方が、より正解データに近い推定結果を出せていることがわかる。
【００７３】
以上説明したように、本発明の実施の形態に係る生存確率推定装置によれば、対象物の各々について観測されたイベントに関するデータに基づいて、ｉ番目のデータの各々に関して、観測が開始されてからイベントが実際に発生した時刻までのイベント発生経過時間Ｔｉ，Ｒを初期化し、経過時間ｔにおいてイベントが発生する確率ｆ（ｔ）及び経過時間ｔまでにイベントが発生していない確率Ｓ（ｔ）を推定するためのモデルに関するパラメータを初期化し、イベントの発生が観測されたｉ番目のデータの各々に関して、観測が開始されてからイベントが発生していないことが最後に確認された時刻までのイベント未発生観測経過時間Ｔｉ，Ｎと、観測が開始されてからイベントが発生したか否かが最後に観測された時刻までの最終観測経過時間Ｔｉと、初期化又は更新されたパラメータとに基づいて、イベント発生経過時間Ｔｉ，Ｒを推定し、ｉ番目のデータの各々についてのイベント発生経過時間Ｔｉ，Ｒ、最終観測経過時間Ｔｉ、及びイベント発生経過時間Ｔｉ，Ｒにより表される尤度式を最適化するように、パラメータを推定することを、予め定めた条件を満たすまで繰り返すことにより、精度よく、対象物の生存確率を推定できる。
【００７４】
また、潜在的なイベント（例：機械の故障など）をモデル化し、イベント発生時期（例：故障機械の実際の故障時期）を推定することによって、対象の本質的な生存期間（例：機械の耐用年数）を測ることが可能になる。例えば機械の故障を例にとると、メンテナンスコストの増大を抑えつつ、機器が故障して利用できないリスクを削減でき、メンテナンススケジュールの最適化が見込まれる。また、機器のより本質的な耐用年数を推定することによって、機器の故障原因の推定、例えば共通部品の経年劣化なのか、機器固有の原因か、人為的な原因か、などをより精度良く推定することも可能となる。
【００７５】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００７６】
例えば上述した実施の形態では、対象物を機械とし、イベントを故障の発生として生存確率を推定する場合について説明したが、これに限定されるものではなく、イベントが発生するまで生存する対象物であればどのようなものでもよい。
【符号の説明】
【００７７】
１０   入力部
２０   演算部
３０   初期化処理部
３２   発生時刻推定部
３４   パラメータ推定部
３６   収束判定部
５０   出力部
１００ 生存確率推定装置
【図１】
2018005580.tif 000013
【図２】
2018005580.tif 000014
【図３】
2018005580.tif 000015
【図４】
2018005580.tif 000016
【図５】
2018005580.tif 000017
【図６】
2018005580.tif 000018
【図７】
2018005580.tif 000019
【図８】
2018005580.tif 000020