(57)【要約】【課題】精度良く離脱確率の分布を推定することができる。【解決手段】データ取得部２４により、複数の顧客のデータを取得し、パラメータ初期値設定部２６により、定常要因ｘの各々について、定常要因ｘによる現象の離脱確率の分布を表す分布パラメータの初期値を設定すると共に、非定常要因ｙの各々について、非定常要因ｙによる現象の離脱確率の分布を表す分布パラメータの初期値を設定し、モデル推定部３０により、顧客の各々について所属する要因を決定し、定常要因ｘの分布パラメータを更新し、非定常要因ｙの分布パラメータを更新し、繰り返し更新された定常要因ｘ毎の分布パラメータに基づいて、定常要因ｘ毎の分布パラメータを決定し、非定常要因分布パラメータ更新手段により繰り返し更新された非定常要因ｙ毎の分布パラメータに基づいて、非定常要因ｙ毎の分布パラメータを決定するパラメータを推定する。【選択図】図１
2015011692.tif 000002
【特許請求の範囲】
【請求項１】
多数の個体に発生する推定対象の現象の離脱確率の分布を推定する生存分析方法であって、
データ取得手段によって、複数の個体ｉの各々について収集された、観測終了時点までの前記現象が継続した観測継続期間ｓi、観測終了時点において前記現象が観測されているか否かを示す情報ｃi、及び観測開始時点から前記現象の現象開始時点までの現象開始前期間ｔiと、絶対的な時期に非定常的に発生する前記現象の離脱要因である非定常要因の各々の発生時点とを取得し、
初期値設定手段によって、前記現象の定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記現象の離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定し、
所属要因決定手段によって、前記複数の個体ｉの各々について、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記個体ｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記個体ｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定し、
定常要因分布パラメータ更新手段によって、前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定された個体ｉの各々の前記情報ｃi及び前記観測継続期間ｓiに基づいて、前記定常要因ｘの前記分布パラメータを更新し、
非定常要因分布パラメータ更新手段によって、前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定された個体ｉの各々の前記情報ｃｉ及び前記観測継続期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新し、
反復判定手段によって、予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返し、
分布パラメータ推定手段によって、前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する、
生存分析方法。
【請求項２】
前記所属要因決定手段によって前記所属する要因を決定することは、
前記個体ｉが非定常要因ｙに所属する度合いを計算するとき、前記個体ｉの観測継続期間ｓｉ及び現象開始前期間ｔｉと、前記非定常要因ｙの発生時点とに基づいて、前記非定常要因ｙの発生時点において、前記個体ｉの前記現象が観測されないと判断された場合に、前記個体ｉが非定常要因ｙに所属する度合いを０とする請求項１記載の生存分析方法。
【請求項３】
離脱数予測手段によって、前記定常要因ｘ毎に、現象が継続している個体ｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記定常要因ｘの前記分布パラメータとに基づいて、入力された予測対象時期における前記定常要因ｘによる離脱数を予測し、
前記非定常要因ｙ毎に、前記現象が継続している個体ｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記非定常要因ｙの前記分布パラメータとに基づいて、前記入力された予測対象時期における前記定常要因ｙによる離脱数を予測することを更に含む請求項１又は２記載の生存分析方法。
【請求項４】
離脱数予測手段によって、仮想的に発生する前記非定常要因ｙについて、現象が継続している個体ｉの各々の観測継続期間ｓｉと、入力された前記非定常要因ｙの仮想的な発生時点と、前記分布パラメータ推定手段によって推定された前記非定常要因ｙの前記分布パラメータとに基づいて、入力された予測対象時期における前記非定常要因ｙによる離脱数を予測することを更に含む請求項１又は２記載の生存分析方法。
【請求項５】
契約サービスのユーザの離脱確率の分布を推定する生存分析方法であって、
データ取得手段によって、複数のユーザｉの各々について収集された、観測終了時点までの前記契約サービスを継続した観測継続期間ｓｉ、観測終了時点において前記契約サービスが契約されているか否かを示す情報ｃｉ、及び観測開始時点から前記契約サービスの初回契約時点までの初回契約前期間ｔｉと、非定常的に発生する前記契約サービスの離脱要因である非定常要因の各々の発生時点とを取得し、
初期値設定手段によって、前記契約サービスの定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記契約サービスの離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定し、
所属要因決定手段によって、前記複数のユーザｉについて、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記ユーザｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記ユーザｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定し、
定常要因分布パラメータ更新手段によって、前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記定常要因ｘの前記分布パラメータを更新し、
非定常要因分布パラメータ更新手段によって、前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新し、
反復判定手段によって、予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返し、
分布パラメータ推定手段によって、前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する
生存分析方法。
【請求項６】
前記所属要因決定手段によって前記所属する要因を決定することは、
前記ユーザｉが非定常要因ｙに所属する度合いを計算するとき、前記ユーザｉの観測継続期間ｓｉ及び初回契約前期間ｔｉと、前記非定常要因ｙの発生時点とに基づいて、前記非定常要因ｙの発生時点において、前記ユーザｉの前記契約サービスが契約されてないと判断された場合に、前記ユーザｉが非定常要因ｙに所属する度合いを０とする請求項５記載の生存分析方法。
【請求項７】
離脱数予測手段によって、前記定常要因ｘ毎に、前記契約サービスが継続しているユーザｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記定常要因ｘの前記分布パラメータとに基づいて、入力された予測対象時期における前記定常要因ｘによる離脱数を予測し、
前記非定常要因ｙ毎に、前記契約サービスが継続しているユーザｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記非定常要因ｙの前記分布パラメータとに基づいて、前記入力された予測対象時期における前記定常要因ｙによる離脱数を予測することを更に含む請求項５又は６記載の生存分析方法。
【請求項８】
離脱数予測手段によって、仮想的に発生する前記非定常要因ｙについて、前記契約サービスが継続しているユーザｉの各々の観測継続期間ｓｉと、入力された前記非定常要因ｙの仮想的な発生時点と、前記分布パラメータ推定手段によって推定された前記非定常要因ｙの前記分布パラメータとに基づいて、入力された予測対象時期における前記非定常要因ｙによる離脱数を予測することを更に含む請求項５又は６記載の生存分析方法。
【請求項９】
前記定常要因ｘによる前記離脱確率の分布は、ワイブル分布であり、前記定常要因ｘの前記分布パラメータは、ワイブル分布の形状パラメータａx及びスケールパラメータθxであり、
前記非定常要因ｙによる前記離脱確率の分布は、指数分布であり、前記非定常要因ｙの前記分布パラメータは、指数分布のパラメータλyである請求項１〜請求項８の何れか１項記載の生存分析方法。
【請求項１０】
多数の個体に発生する推定対象の現象の離脱確率の分布を推定する生存分析装置であって、
複数の個体ｉの各々について収集された、観測終了時点までの前記現象が継続した観測継続期間ｓi、観測終了時点において前記現象が観測されているか否かを示す情報ｃi、及び観測開始時点から前記現象の現象開始時点までの現象開始前期間ｔiと、絶対的な時期に非定常的に発生する前記現象の離脱要因である非定常要因の各々の発生時点とを取得するデータ取得手段と、
前記現象の定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記現象の離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定する初期値設定手段と、
前記複数の個体ｉの各々について、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記個体ｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記個体ｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定する所属要因決定手段と、
前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定された個体ｉの各々の前記情報ｃi及び前記観測継続期間ｓiに基づいて、前記定常要因ｘの前記分布パラメータを更新する定常要因分布パラメータ更新手段と、
前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定された個体ｉの各々の前記情報ｃｉ及び前記観測継続期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新する非定常要因分布パラメータ更新手段と、
予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返す反復判定手段と、
前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する分布パラメータ推定手段と、
を含む、生存分析装置。
【請求項１１】
契約サービスのユーザの離脱確率の分布を推定する生存分析装置であって、
複数のユーザｉの各々について収集された、観測終了時点までの前記契約サービスを継続した観測継続期間ｓｉ、観測終了時点において前記契約サービスが契約されているか否かを示す情報ｃｉ、及び観測開始時点から前記契約サービスの初回契約時点までの初回契約前期間ｔｉと、非定常的に発生する前記契約サービスの離脱要因である非定常要因の各々の発生時点とを取得するデータ取得手段と、
前記契約サービスの定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記契約サービスの離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定する初期値設定手段と、
前記複数のユーザｉについて、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記ユーザｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記ユーザｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定する所属要因決定手段と、
前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記定常要因ｘの前記分布パラメータを更新する定常要因分布パラメータ更新手段と、
前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新する非定常要因分布パラメータ更新手段と、
予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返す反復判定手段と、
前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する分布パラメータ推定手段と、
を含む、生存分析装置。
【請求項１２】
コンピュータに、請求項１〜９の何れか１項記載の生存分析方法を構成する各ステップを実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、多数の個体に発生する推定対象の現象の離脱確率の分布を推定する生存分析方法、装置、及びプログラムに関する。
【背景技術】
【０００２】
回線契約等の契約型サービスにおいて、ある時点に発生した競合のキャンペーンや集中営業がユーザの離脱にどの程度の影響を有しているかを測る必要がある。
【０００３】
従来、現象発生からの経過期間（生存期間）を複数のワイブル分布モデルの混合値としてモデル化することで、定常的な生存期間分布を推定している（非特許文献１、非特許文献２）。
【０００４】
また、特定の絶対時間に発生する非定常な離脱を分析し、非定常なイベントを検出する方法が知られている（非特許文献３）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】”実践編３  生存分析”、［online］［ 平成２５年６月２４日検索］、インターネット＜URL:http://www.med.nagoya-u.ac.jp/spss/files/jissen3.pdf＞
【非特許文献２】J.M. Marin, M.R. Rodriguez Bernal and M.P. Wiper, USING WEIBULL MIXTURE DISTRIBUTIONS TO MODEL HETEROGENEOUS SURVIVAL DATA, Universidad Carlos III de Madrid Working Paper 03-32, Statistics and Econometrics Series 08, 2003
【非特許文献３】山本  拓、”時系列分析とその経済分析への応用”、［online］［ 平成２５年６月２４日検索］、インターネット＜URL: http://www.mof.go.jp/pri/publication/financial_review/fr_list2/r23/r_23_048_072.pdf＞
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかし、世の中に発生する事象としては、ある特定の絶対時間に開始されるイベントがあり（例えば、他社サービスの開始等）、これらのイベントに基づいて発生する事象（以下、非定常要因とする）を、定常要因とみなした従来のモデルの上に分類することは困難であり、分離を図ってモデル化しなければならないという問題がある。
【０００７】
また、分離モデルを構築する際、既モデルのように全てのユーザに適用するモデルを構築した場合、非定常要因発生前に離脱したユーザや非定常要因発生後に開始したユーザが、非定常要因による離脱と判定される可能性があるため、非定常要因が過大に評価されるという問題がある。
【０００８】
また、非特許文献３については、契約開始期間からの経過に表出する定常的な生存期間を扱えず、契約日の集中による離脱多発を誤検出してしまうという問題がある。
【０００９】
本発明は、上記問題点を解決するために成されたものであり、精度良く離脱確率の分布を推定することができる生存分析方法、装置、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
上記目的を達成するために、第１の発明の生存分析方法は、多数の個体に発生する推定対象の現象の離脱確率の分布を推定する生存分析方法であって、データ取得手段によって、複数の個体ｉの各々について収集された、観測終了時点までの前記現象が継続した観測継続期間ｓi、観測終了時点において前記現象が観測されているか否かを示す情報ｃi、及び観測開始時点から前記現象の現象開始時点までの現象開始前期間ｔiと、絶対的な時期に非定常的に発生する前記現象の離脱要因である非定常要因の各々の発生時点とを取得し、初期値設定手段によって、前記現象の定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記現象の離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定し、所属要因決定手段によって、前記複数の個体ｉの各々について、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記個体ｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記個体ｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定し、定常要因分布パラメータ更新手段によって、前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定された個体ｉの各々の前記情報ｃi及び前記観測継続期間ｓiに基づいて、前記定常要因ｘの前記分布パラメータを更新し、非定常要因分布パラメータ更新手段によって、前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定された個体ｉの各々の前記情報ｃｉ及び前記観測継続期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新し、反復判定手段によって、予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返し、分布パラメータ推定手段によって、前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する。
【００１１】
第２の発明の生存分析装置は、多数の個体に発生する推定対象の現象の離脱確率の分布を推定する生存分析装置であって、複数の個体ｉの各々について収集された、観測終了時点までの前記現象が継続した観測継続期間ｓi、観測終了時点において前記現象が観測されているか否かを示す情報ｃi、及び観測開始時点から前記現象の現象開始時点までの現象開始前期間ｔiと、絶対的な時期に非定常的に発生する前記現象の離脱要因である非定常要因の各々の発生時点とを取得するデータ取得手段と、前記現象の定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記現象の離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定する初期値設定手段と、前記複数の個体ｉの各々について、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記個体ｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記個体ｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定する所属要因決定手段と、前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定された個体ｉの各々の前記情報ｃi及び前記観測継続期間ｓiに基づいて、前記定常要因ｘの前記分布パラメータを更新する定常要因分布パラメータ更新手段と、
前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定された個体ｉの各々の前記情報ｃｉ及び前記観測継続期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新する非定常要因分布パラメータ更新手段と、予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返す反復判定手段と、前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する分布パラメータ推定手段と、を含んで構成されている。
【００１２】
第１の発明及び第２の発明によれば、多数の個体に発生する現象の離脱確率の分布を精度よく推定することができる。
【００１３】
第１の発明の生存分析方法によれば、前記所属要因決定手段によって前記所属する要因を決定することは、前記個体ｉが非定常要因ｙに所属する度合いを計算するとき、前記個体ｉの観測継続期間ｓｉ及び現象開始前期間ｔｉと、前記非定常要因ｙの発生時点とに基づいて、前記非定常要因ｙの発生時点において、前記個体ｉの前記現象が観測されないと判断された場合に、前記個体ｉが非定常要因ｙに所属する度合いを０とすることができる。
【００１４】
第１の発明の生存分析方法によれば、離脱数予測手段によって、前記定常要因ｘ毎に、現象が継続している個体ｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記定常要因ｘの前記分布パラメータとに基づいて、入力された予測対象時期における前記定常要因ｘによる離脱数を予測し、前記非定常要因ｙ毎に、前記現象が継続している個体ｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記非定常要因ｙの前記分布パラメータとに基づいて、前記入力された予測対象時期における前記定常要因ｙによる離脱数を予測することを更に含むことができる。
【００１５】
第１の発明の生存分析方法によれば、離脱数予測手段によって、仮想的に発生する前記非定常要因ｙについて、現象が継続している個体ｉの各々の観測継続期間ｓｉと、入力された前記非定常要因ｙの仮想的な発生時点と、前記分布パラメータ推定手段によって推定された前記非定常要因ｙの前記分布パラメータとに基づいて、入力された予測対象時期における前記非定常要因ｙによる離脱数を予測することを更に含むことができる。
【００１６】
第３の発明の生存分析方法は、契約サービスのユーザの離脱確率の分布を推定する生存分析方法であって、データ取得手段によって、複数のユーザｉの各々について収集された、観測終了時点までの前記契約サービスを継続した観測継続期間ｓｉ、観測終了時点において前記契約サービスが契約されているか否かを示す情報ｃｉ、及び観測開始時点から前記契約サービスの初回契約時点までの初回契約前期間ｔｉと、非定常的に発生する前記契約サービスの離脱要因である非定常要因の各々の発生時点とを取得し、初期値設定手段によって、前記契約サービスの定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記契約サービスの離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定し、所属要因決定手段によって、前記複数のユーザｉについて、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記ユーザｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記ユーザｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定し、定常要因分布パラメータ更新手段によって、前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記定常要因ｘの前記分布パラメータを更新し、非定常要因分布パラメータ更新手段によって、前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新し、反復判定手段によって、予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返し、分布パラメータ推定手段によって、前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する。
【００１７】
第４の発明の生存分析装置は、  契約サービスのユーザの離脱確率の分布を推定する生存分析装置であって、複数のユーザｉの各々について収集された、観測終了時点までの前記契約サービスを継続した観測継続期間ｓｉ、観測終了時点において前記契約サービスが契約されているか否かを示す情報ｃｉ、及び観測開始時点から前記契約サービスの初回契約時点までの初回契約前期間ｔｉと、非定常的に発生する前記契約サービスの離脱要因である非定常要因の各々の発生時点とを取得するデータ取得手段と、前記契約サービスの定常的な離脱要因である定常要因ｘの各々について、前記定常要因ｘによる前記契約サービスの離脱確率の分布を表わす分布パラメータの初期値を設定すると共に、前記非定常要因ｙの各々について、前記非定常要因ｙによる前記現象の離脱確率の分布を表す分布パラメータの初期値を設定する初期値設定手段と、前記複数のユーザｉについて、前記定常要因ｘ毎に、前記定常要因ｘの前記分布パラメータに基づいて、前記ユーザｉが定常要因ｘに所属する度合いを計算し、前記非定常要因ｙ毎に、前記非定常要因ｙの前記分布パラメータに基づいて、前記ユーザｉが非定常要因ｙに所属する度合いを計算し、前記計算された前記度合いの各々に基づいて、前記定常要因ｘの各々及び前記非定常要因ｙの各々のうちの何れか一つを、前記所属する要因として決定する所属要因決定手段と、前記定常要因ｘ毎に、前記定常要因ｘに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記定常要因ｘの前記分布パラメータを更新する定常要因分布パラメータ更新手段と、前記非定常要因ｙ毎に、前記非定常要因ｙに所属すると決定されたユーザｉの各々の前記情報ｃｉ及び前記観測期間ｓｉに基づいて、前記非定常要因ｙの前記分布パラメータを更新する非定常要因分布パラメータ更新手段と、予め定められた反復終了条件を満足するまで、前記所属要因決定手段、前記定常要因分布パラメータ更新手段、及び前記非定常要因分布パラメータ更新手段を繰り返す反復判定手段と、前記定常要因分布パラメータ更新手段により繰り返し更新された前記定常要因ｘ毎の前記分布パラメータに基づいて、前記定常要因ｘ毎の前記分布パラメータを決定し、前記非定常要因分布パラメータ更新手段により繰り返し更新された前記非定常要因ｙ毎の前記分布パラメータに基づいて、前記非定常要因ｙ毎の前記分布パラメータを決定する分布パラメータ推定手段と、を含んで構成されている。
【００１８】
第３及び第４の発明によれば、契約サービスについて離脱確率の分布を精度よく推定することができる。
【００１９】
第３の発明の生存分析方法によれば、前記所属要因決定手段によって前記所属する要因を決定することは、前記ユーザｉが非定常要因ｙに所属する度合いを計算するとき、前記ユーザｉの観測継続期間ｓｉ及び初回契約前期間ｔｉと、前記非定常要因ｙの発生時点とに基づいて、前記非定常要因ｙの発生時点において、前記ユーザｉの前記契約サービスが契約されてないと判断された場合に、前記ユーザｉが非定常要因ｙに所属する度合いを０とすることができる。
【００２０】
第３の発明の生存分析方法は、離脱数予測手段によって、前記定常要因ｘ毎に、前記契約サービスが継続しているユーザｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記定常要因ｘの前記分布パラメータとに基づいて、入力された予測対象時期における前記定常要因ｘによる離脱数を予測し、前記非定常要因ｙ毎に、前記契約サービスが継続しているユーザｉの各々の観測継続期間ｓｉと、前記分布パラメータ推定手段によって決定された前記非定常要因ｙの前記分布パラメータとに基づいて、前記入力された予測対象時期における前記定常要因ｙによる離脱数を予測することを更に含むことができる。
【００２１】
第３の発明の生存分析方法は、離脱数予測手段によって、仮想的に発生する前記非定常要因ｙについて、前記契約サービスが継続しているユーザｉの各々の観測継続期間ｓｉと、入力された前記非定常要因ｙの仮想的な発生時点と、前記分布パラメータ推定手段によって推定された前記非定常要因ｙの前記分布パラメータとに基づいて、入力された予測対象時期における前記非定常要因ｙによる離脱数を予測することを更に含むことができる。
【００２２】
第１及び第３の発明の生存分析方法によれば、前記定常要因ｘによる前記離脱確率の分布は、ワイブル分布であり、前記定常要因ｘの前記分布パラメータは、ワイブル分布の形状パラメータａx及びスケールパラメータθxであり、前記非定常要因ｙによる前記離脱確率の分布は、指数分布であり、前記非定常要因ｙの前記分布パラメータは、指数分布のパラメータλyとすることができる。
【００２３】
また、本発明のプログラムは、コンピュータに、上記の生存分析方法を構成する各ステップを実行させるためのプログラムである。
【発明の効果】
【００２４】
以上説明したように、生存分析方法、装置、及びプログラムによれば、精度良く離脱確率の分布を推定することができる。
【図面の簡単な説明】
【００２５】
【図１】本発明の第１の実施の形態に係る生存分析装置の機能的構成を示すブロック図である。
【図２】発生情報の例を示す図である。
【図３】仮想発生情報の例を示す図である。
【図４】顧客契約情報の例を示す図である。
【図５】入力部１０及びデータ取得部２４により取得した情報を表した例を示す図である。
【図６】出力結果の例を示す図である。
【図７】本発明の第１の実施の形態に係る生存分析装置における生存分析処理ルーチンを示すフローチャートである。
【図８】本発明の第１の実施の形態に係る生存分析装置におけるパラメータ推定処理ルーチンを示すフローチャートである。
【図９】本発明の第１の実施の形態に係る生存分析装置における顧客所属パラメータ更新処理ルーチンを示すフローチャートである。
【図１０】本発明の第１の実施の形態に係る生存分析装置における重みパラメータ更新処理ルーチンを示すフローチャートである。
【図１１】本発明の第１の実施の形態に係る生存分析装置におけるワイブルスケールパラメータ更新処理ルーチンを示すフローチャートである。
【図１２】本発明の第１の実施の形態に係る生存分析装置におけるワイブル形状パラメータ更新処理ルーチンを示すフローチャートである。
【図１３】本発明の第１の実施の形態に係る生存分析装置における指数パラメータ更新処理ルーチンを示すフローチャートである。
【図１４】本発明の第１の実施の形態に係る生存分析装置における結果生成処理ルーチンを示すフローチャートである。
【図１５】本発明の第２の実施の形態に係る生存分析装置におけるパラメータ推定処理ルーチンを示すフローチャートである。
【発明を実施するための形態】
【００２６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２７】
＜第１の実施の形態に係る生存分析装置の構成＞
まず、本発明の第１の実施の形態に係る生存分析装置の構成について説明する。図１に示すように、本発明の実施の形態に係る生存分析装置１００は、ＣＰＵとＲＡＭと後述する生存分析処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することができる。この生存分析装置１００は、機能的には図１に示すように入力部１０と、演算部２０と、出力部５０とを備えている。
【００２８】
入力部１０は、キーボードなどの入力装置から分析対象となる顧客ｉ（１≦ｉ≦Ｎ）、非定常要因の各々についてのデータ取得開始日から非定常要因の各々の発生までの期間Ｔ（Ｔ１，Ｔ２，…，ＴＹ）を表すＹ個の発生情報、及び仮想非定常要因の各々についてのデータ取得開始日から仮想非定常要因の各々の発生までの期間である発生時期Ｕ（Ｕ１，Ｕ２，…，Ｕｊ）と発生させる仮想非定常要因の各々の仮想非定常要因Ｚ（Ｚ１，Ｚ２，…，Ｚｊ）とからなるｊ個の仮想発生情報を受け付ける。なお、入力部１０は、ネットワーク等を介して外部から入力されたものを受け付けるようにしてもよい。また、Ｔ及びＵは０より大きい経過週数の数値となり、Ｚは発生情報で指定した（あるいは過去に行った分析の）非定常要因の名前である。また、Ｙ及びｊは０以上の整数である。
【００２９】
ここで、発生情報は離脱を促すイベントが発生したと思われるＹ個の時期を入力し、Ｙ個の時期の各々とデータ取得開始日との間の期間を算出することで獲得する。また、仮想発生情報は、将来的に離脱を促すイベントが発生すると思われるｊ個の時期Ｕを入力し、ｊ個の時期の各々とデータ取得開始日との間の期間を算出することで獲得し、発生するイベントである仮想非定常要因Ｚを指定することで獲得する。図２に発生情報の例を示し、図３に仮想発生情報の例を示す。
【００３０】
演算部２０は、記憶部２２と、データ取得部２４と、パラメータ初期値設定部２６と、モデル推定部３０と、パラメータ追加部３２と、予測部３４と、を含んだ構成で表すことができる。
【００３１】
記憶部２２は、全ての顧客の各々についての観測継続期間ｓ、現在の状態ｃ、及びデータ取得開始日から初回契約までの初回契約前期間ｔからなる顧客契約情報が記憶されている。観測継続期間は、観測終了時点までの契約サービスを継続した期間であり、現在の情報ｃは、観測終了時点において契約サービスが観測されているか否かを示す情報である。なお、ｓ及びｔは０より大きい経過週数であり、ｃは１か０の値（ｃ＝１が離脱を表し、ｃ＝０が継続を表す）をとる。図４に顧客契約情報の例を示す。
【００３２】
データ取得部２４は、入力部１０において受け付けた分析対象となる顧客ｉの各々を取得し、当該顧客ｉの各々の顧客契約情報（ｓｉ、ｃｉ及びｔｉ）を記憶部２２から取得する。なお、入力部１０において受け付けた情報と、データ取得部２４により取得した情報を表した例を図５に示す。
【００３３】
パラメータ初期値設定部２６は、定常分布数（以下Ｘとし、定常要因数と同数であり、任意の値とする）、及び非定常分布数（以下Ｙとし、入力部１０において受け付けた非定常要因の数と同値とする）を設定する。また、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ及びワイブルスケールパラメータθｘ、Ｘ＋Ｙ個の要因ｋの重みパラメータｗk、各顧客ｉの顧客所属パラメータｈｉ、Ｙ個の非定常要因ｙの指数パラメータλｙ、ａの事前パラメータ（αａ，βａ）、θの事前パラメータ（αθ，βθ）、及びｗの事前ディリクレパラメータφの各々に、適切な初期値を設定する。第１の実施の形態においては、ｘ＝２、ｙ＝１、ａ＝（１，１）、θ＝（１，１）、ｗ＝（１／３，１／３，１／３）、λ＝１、（αａ，βａ）＝（１，１）、（αθ，βθ）＝（１，１）、φ＝［１，１，…，１］の１行Ｎ列を、初期値として設定し、ｈは１〜Ｘ＋Ｙの離散値をランダムに与えるものとする。なお、第１の実施の形態においては、定常要因による顧客の離脱が観測継続期間ｓ上のワイブル分布に従うものとし、非定常要因による顧客の離脱が、データ取得開始日からの期間ｔ上の位置パラメータ付き指数分布に従うものと仮定する。
【００３４】
モデル推定部３０は、パラメータ初期値設定部２６において設定された各パラメータの初期値に基づいて、マルコフ連鎖モンテカルロ法（ＭＣＭＣ：Markovchain Monte Carlo methods）に従って、分析対象の顧客ｉの各々の顧客所属パラメータｈｉ、定常要因又は非定常要因ｋの各々の重みパラメータｗk、定常要因ｘの各々のワイブルスケールパラメータθｘ、定常要因ｘの各々のワイブル形状パラメータａｘ、及び非定常要因ｙの各々の指数パラメータλｙを推定する。具体的には、サイクル数＝１５０００、バーンインサイクル数＝１００００と設定し、１０００１回〜１５０００回の各サイクルのパラメータの各々をサイクル毎に記憶し、パラメータ毎に、サイクル毎に記憶された当該パラメータの値の平均値を当該パラメータの値として推定する。
【００３５】
モデル推定部３０において推定された分析対象の顧客の各々の顧客所属パラメータｈｉ、定常要因又は非定常要因ｋの各々の重みパラメータｗk、定常要因ｘの各々のワイブルスケールパラメータθｘ、定常要因ｘの各々のワイブル形状パラメータａｘ、及び非定常要因ｙの各々の指数パラメータλｙを用いると、離脱確率の分布を求めることができる。定常要因ｘの定常パラメータを、param１(ａｘ，θｘ)とし、非定常要因ｙの非定常パラメータをparam２（λｙ，Ｔｙ）とすると、モデル推定部３０において推定された各種パラメータを、下記（１）式に示す基本モデル式に代入することにより、顧客の離脱確率の分布を求めることができる。ここで、下記（１）式の１項目が定常要因に関わる部分であり、２項目が非定常要因に関わる部分である。
【００３６】
【数１】
2015011692.tif 000003
【００３７】
パラメータ追加部３２は、入力部１０において受け付けた仮想発生情報と、モデル推定部３０により推定された各パラメータの各々の値とに基づいて、仮想非定常要因のパラメータを追加する。具体的には、仮想発生情報に含まれる仮想非定常要因Ｚの各々について、当該仮想非定常要因Ｚに対応する非定常要因を特定し、モデル推定部３０によるパラメータの推定結果から、当該非定常要因の非定常要因パラメータλと、当該非定常要因パラメータλに対応する重みパラメータｗを取得する。そして、仮想非定常要因Ｚの各々について、取得した非定常要因パラメータλのＴを、当該仮想非定常要因Ｚの発生時期Ｕに置き換え、新たな非定常パラメータλ（ｊ個の仮想非定常パラメータとする）とその非定常パラメータに対応する重みパラメータｗを作成する。
【００３８】
パラメータ追加部３２は、作成したｊ個の仮想非定常パラメータをＹ個の非定常パラメータの末尾に追加し、作成したｊ個の重みパラメータをＸ＋Ｙ個の重みパラメータの末尾に追加する。これにより、非定常パラメータの数がＹ＋ｊ個となり、重みパラメータはＸ＋Ｙ＋ｊ個となる。なお、重みパラメータがｊ個追加されたため、例えば、ｗ１の値はｗ１／Σｗとするように、Ｘ＋Ｙ＋ｊ個の重みパラメータｗの合計が１となるように重みパラメータｗを正規化する。
【００３９】
予測部３４は、入力部１０において受け付けた分析対象となる顧客のうち、観測終了時である現在において契約を継続している顧客ｉ´の数Ｎ´（ｃ＝０の顧客数）をカウントする。ここで、発生情報又は仮想発生情報で指定されたＬ番目の非定常要因の総離脱数はｗX＋Ｌ×Ｎ´となる。
【００４０】
予測部３４は、顧客ｉ´（０≦ｉ´≦Ｎ´）の各々について、各週、各要因の離脱確率を算出し、要因ごとの各週期待離脱者数を集計し、出力部５０に出力する。なお、ある週の期待離脱者数は、各顧客ｉ´の当該週の離脱確率を加算することにより求められる。具体的には、顧客ｉ´がＷ週目にｘ番目の定常要因で離脱する確率を下記（２）式により求め、顧客ｉ´がＷ週目にｙ´番目（０≦ｙ´≦Ｙ´）の非定常要因で離脱する確率を下記（３）式により求める。また、ｙ´番目の非定常要因が仮想非定常要因の場合には、下記（３）式のＴはＵの値を用いる。図６に出力結果の例を示す。図６の出力結果については、小数点以下の値を切り捨てているものとする。ここで、Ｙ´＝Ｙ＋ｊである。
【００４１】
【数２】
2015011692.tif 000004
【００４２】
＜第１の実施の形態に係る生存分析装置の作用＞
次に、本発明の実施の形態に係る生存分析装置１００の作用について説明する。図７に示す生存分析処理ルーチンを実行する前に、入力部１０により全顧客の各々の顧客契約情報を受け付け記憶部２２に記憶する。そして、入力部１０により分析対象となる顧客、発生情報、及び仮想発生情報を受け付けると、生存分析装置１００は、図７に示す生存分析処理ルーチンを実行する。
【００４３】
まず、ステップＳ１００では、入力部１０により受け付けた分析対象となる顧客ｉ、発生情報、及び仮想発生情報の各々を受け付ける。
【００４４】
次に、ステップＳ１０２では、記憶部２２に記憶されているステップＳ１００において受け付けた分析対象となる顧客ｉの各々の顧客契約情報を読み込む。
【００４５】
次に、ステップＳ１０４では、ステップＳ１００において取得した発生情報と、ステップＳ１０２において取得した分析対象となる顧客ｉの各々の顧客契約情報とに基づいて、生存分析をするためのパラメータの各々を推定する。
【００４６】
次に、ステップＳ１０６では、ステップＳ１００において取得した仮想発生情報と、ステップＳ１０４において取得した生存分析をするためのパラメータの各々とに基づいて、生存分析をした結果を生成する。
【００４７】
次に、ステップＳ１０８では、ステップＳ１０６において取得した生存分析をした結果を出力部５０に出力して処理を終了する。
【００４８】
上記ステップＳ１０４は、図８に示すパラメータ推定処理ルーチンによって実現される。
【００４９】
まず、ステップＳ２００では、定常分布数（Ｘ）、非定常分布数（Ｙ）、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及び定常要因ｘの各々のワイブルスケールパラメータθｘ、Ｘ＋Ｙ個の要因ｋの重みパラメータｗｋ、各顧客ｉの顧客所属パラメータｈｉ、Ｙ個の非定常要因ｙの指数パラメータλｙ、ａの事前パラメータ（αａ，βａ）、θの事前パラメータ（αθ，βθ）、及びｗの事前ディリクレパラメータφの各々に、適切な初期値を設定する。また、ＭＣＭＣサイクルのパラメータであるサイクル数のパラメータＧを１５０００とし、バーンインサイクルburninを１００００の値に設定する。
【００５０】
次に、ステップＳ２０２では、変数ｇの値に１を設定する。
【００５１】
次に、ステップＳ２０４では、顧客を識別するための変数ｉの値に１を設定する。
【００５２】
次に、ステップＳ２０６では、ステップＳ２００において取得又は前回のサイクルの処理において更新したＸ＋Ｙ個の要因ｋの重みパラメータｗｋ、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及びワイブルスケールパラメータθｘ、及びＹ個の非定常要因ｙの指数パラメータλｙと、ステップＳ１０２において取得した顧客ｉの顧客契約情報と、ステップＳ１００において取得した発生情報とに基づいて、顧客ｉの顧客所属パラメータｈｉの値を更新する。
【００５３】
上記ステップＳ２０６は、図９に示す顧客所属パラメータ更新処理ルーチンによって実現される。
【００５４】
図９のステップＳ３００では、変数ｋの値に１を設定する。
【００５５】
次に、ステップＳ３０２では、ステップＳ２００において取得又は前回のサイクルの処理において更新したＸ＋Ｙ個の要因ｋの重みパラメータｗｋ、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及びワイブルスケールパラメータθｘ、及びＹ個の非定常要因ｙの指数パラメータλｙと、ステップＳ１０２において取得した顧客ｉの顧客契約情報と、ステップＳ１００において取得した発生情報とに基づいて、顧客ｉのｋ番目の定常要因又は非定常要因（Ｘ個の定常要因を並べた後に、Ｙ個の非定常要因を並べているものとする）の所属のしやすさであるｂａｒ（ｉ，ｋ）を下記（４）式により算出する。
【００５６】
【数３】
2015011692.tif 000005
【００５７】
上記（４）式に示すように、非定常要因の所属のしやすさを算出する場合には、当該非定常要因の発生時期に、顧客ｉの契約サービスが継続されているか否かを判断し、顧客ｉの契約サービスが継続されていない場合には、非定常要因の所属のしやすさを０とする。
【００５８】
ステップＳ３０４では、変数ｋの値がＸ＋Ｙの値と一致するか否かの判定を行う。変数ｋの値がＸ＋Ｙの値と一致する場合にはステップＳ３０８へ移行し、変数ｋの値がＸ＋Ｙの値と一致しない場合には、ステップＳ３０６へ移行する。
【００５９】
ステップＳ３０６では、変数ｋの値をｋ＋１の値に設定し、ステップＳ３０２に移行し、ステップＳ３０２〜ステップＳ３０４の処理を繰り返す。
【００６０】
ステップＳ３０８では、ステップＳ３０２において取得したｂａｒ（ｉ，ｋ）の各々の値に基づいて、重みづけルーレットにより乱数を発生させる。
【００６１】
ステップＳ３１０では、ステップＳ３０８において取得した乱数により顧客ｉの所属する定常要因又は非定常要因を決定し、顧客所属パラメータｈｉを更新し顧客所属パラメータ更新処理ルーチンを終了する。
【００６２】
次に、図８のステップＳ２０８では、変数ｉの値をｉ＋１の値に設定する。
【００６３】
次に、ステップＳ２１０では、変数ｉの値がＮ＋１の値と一致するか否か判定する。変数ｉの値がＮ＋１の値と一致する場合には、ステップＳ２１２へ移行し、変数ｉの値がＮ＋１の値と一致しない場合には、ステップＳ２０６へ移行し、ステップＳ２０６〜ステップＳ２１０の処理を繰り返す。
【００６４】
次に、ステップＳ２１２では、ステップＳ２００において取得した重みパラメータの事前ディリクレパラメータφと、ステップＳ２０６において更新した顧客ｉの各々の顧客所属パラメータｈｉとに基づいて、Ｘ＋Ｙ個の要因ｋの重みパラメータｗｋを更新する。
【００６５】
上記ステップＳ２１２は、図１０に示す重みパラメータ更新処理ルーチンによって実現される。
【００６６】
図１０のステップＳ４００では、変数ｋの値に１を設定する。
【００６７】
次に、ステップＳ４０２では、ｋ番目の定常要因又は非定常要因（Ｘ個の定常要因を並べた後に、Ｙ個の非定常要因を並べているものとする）に所属しているユーザの数を、ステップＳ３０８において取得した顧客ｉの各々の顧客所属パラメータｈｉの値から求め、求めた値をｗｅｉｇｈｔｋとする。
【００６８】
次に、ステップＳ４０４では、ステップＳ２００において取得した重みパラメータの事前ディリクレパラメータφと、ステップＳ４０２において取得したｗｅｉｇｈｔｋとに基づいて、（φ＋ｗｅｉｇｈｔｋ，１）をパラメータとするガンマ乱数を発生させ、発生した乱数を重みパラメータｗｋとする。
【００６９】
次に、ステップＳ４０６では、変数ｋの値がＸ＋Ｙの値と一致しているか否かの判定をする。変数ｋの値がＸ＋Ｙの値と一致している場合には、ステップＳ４１０へ移行し、変数ｋの値がＸ＋Ｙの値と一致していない場合には、ステップＳ４０８へ移行する。
【００７０】
ステップＳ４０８では、変数ｋの値をｋ＋１と設定し、ステップＳ４０２へ移行し、ステップＳ４０２〜ステップＳ４０６の処理を繰り返す。
【００７１】
ステップＳ４１０では、ステップＳ４０４において取得した定常要因及び非定常要因の各々についての重みパラメータｗｋの合計が１になるように正規化し、正規化後の値を、定常要因及び非定常要因の各々についての重みパラメータｗｋとし、重みパラメータ更新処理ルーチンの処理を終了する。
【００７２】
図８のステップ２１４では、ステップＳ２００において取得したワイブルスケールパラメータの事前パラメータ（αθ，βθ）と、ステップＳ２００において取得し又は前回のサイクルの処理において更新したＸ個の定常要因ｘの各々のワイブル形状パラメータａｘ、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６において更新した顧客ｉの顧客所属パラメータｈｉの各々とに基づいて、Ｘ個の定常要因ｘの各々のワイブルスケールパラメータθを更新する。
【００７３】
上記ステップＳ２１４は、図１１に示すワイブルスケールパラメータ更新処理ルーチンによって実現される。
【００７４】
図１１のステップＳ５００では、変数ｘの値に１を設定する。
【００７５】
次に、ステップＳ５０１では、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６で更新した顧客ｉの各々の顧客所属パラメータｈｉとに基づいて、Ｒｘを取得する。ここで、Ｒｘとは、ｘ番目の定常要因（Ｘ個の定常要因を並べているものとする）のｃｉ＝１の顧客数である。
【００７６】
次に、ステップＳ５０２では、ステップＳ２００において取得したワイブルスケールパラメータの事前パラメータ（αθ，βθ）と、ステップＳ２００において取得し又は前回にサイクルの処理において更新したＸ個の定常要因ｘの各々のワイブル形状パラメータａｘと、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６において更新した顧客ｉの各々の顧客所属パラメータｈｉと、ステップＳ２０６において取得したＲｘとの各々に基づいて、（Ｒｘ＋αθ），（βθ＋Σ（ｓｉ（ｈｉ＝＝ｋ）＾ａｘ））をパラメータとするガンマ乱数を発生させ、発生したガンマ乱数をワイブルスケールパラメータθｘとして更新する。
【００７７】
次に、ステップＳ５０４では、変数ｘの値がＸと一致するか否かの判定を行う。変数ｘの値がＸの値と一致する場合には、ワイブルスケールパラメータ更新処理ルーチンを終了し、変数ｘの値がＸの値と一致しない場合には、ステップＳ５０６へ移行する。
【００７８】
ステップＳ５０６では、変数ｘの値をｘ＋１と設定し、ステップＳ５０１へ移行し、ステップＳ５０１〜ステップＳ５０４の処理を繰り返す。
【００７９】
図８のステップＳ２１６では、ステップＳ２００において取得したワイブル形状パラメータの事前パラメータ（αａ，βａ）と、ステップＳ２１４において更新したＸ個の定常要因ｘの各々のワイブル形状パラメータａｘと、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６において更新した顧客ｉの各々の顧客契約パラメータｈｉの各々とに基づいて、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘを更新する。
【００８０】
上記ステップＳ２１６は、図１２に示すワイブル形状パラメータ更新処理ルーチンによって実現される。
【００８１】
図１２のステップＳ６００では、変数ｘの値を１に設定する。
【００８２】
次に、ステップＳ６０２では、ステップＳ２００において取得し又は前回のサイクルの処理において更新した定常要因ｘのワイブル形状パラメータａｘに標準正規乱数を加えた値を、ｎｅｗａとする。
【００８３】
次に、ステップＳ６０４では、ステップＳ２００において取得したワイブル形状パラメータの事前パラメータ（αａ，βａ）と、ステップＳ２００において取得し又は前回のサイクルの処理において更新した定常要因ｘのワイブル形状パラメータａｘと、ステップＳ６０２で得られたｎｅｗａと、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６において更新した顧客ｉの各々の顧客所属パラメータｈｉとに基づいて、ｎｅｗａの尤度、及びａｘの尤度の各々を、下記（５）式に従って算出し、（（ｎｅｗａの尤度）／（ａｘの尤度））であるνの値を算出する。
【００８４】
【数４】
2015011692.tif 000006
【００８５】
次に、ステップＳ６０６では、一様乱数Ａ（０＜Ａ＜１）を発生させる。
【００８６】
次に、ステップＳ６０８では、ステップＳ６０４で算出したνの値がステップＳ６０６で取得した一様乱数Ａの値よりも大きいか否かの判定を行う。νの値が一様乱数Ａよりも大きい場合には、ステップＳ６１０へ移行し、νの値が一様乱数Ａ以下の場合には、ステップＳ６１２へ移行する。
【００８７】
ステップＳ６１０では、ワイブル形状パラメータａｘの値をステップＳ６０２において取得したｎｅｗａの値に更新する。
【００８８】
次に、ステップＳ６１２では、変数ｘの値がＸの値と一致するか否かの判定を行う。変数ｘの値がＸの値と一致する場合には、ワイブル形状パラメータ更新処理ルーチンを終了し、変数ｘの値がＸの値と一致しない場合には、ステップＳ６１４へ移行する。
【００８９】
ステップＳ６１４では、変数ｘの値をｘ＋１に設定し、ステップＳ６０２へ移行し、ステップＳ６０２〜ステップＳ６１２の処理を繰り返す。
【００９０】
図８のステップＳ２１８では、ステップＳ２００において取得し又は前回のサイクルの処理において更新したＹ個の非定常要因ｙの各々の指数パラメータλｙと、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６において更新した顧客ｉの各々の顧客所属パラメータｈｉとに基づいて、Ｙ個の非定常要因ｙの各々の指数パラメータλｙの値を更新する。
【００９１】
上記ステップＳ２１８は、図１３に示す指数パラメータ更新処理ルーチンによって実現される。
【００９２】
図１３のステップ７００では、変数ｙの値を１に設定する。
【００９３】
次に、ステップＳ７０２では、ステップＳ２００において取得し又は前回のサイクルの処理において更新した非定常要因ｙの指数パラメータλｙに、標準正規乱数を加えた値を、ｎｅｗλとする。
【００９４】
次に、ステップＳ７０４では、ステップＳ２００において取得し又は前回のサイクルの処理において更新した非定常要因ｙの指数パラメータλｙと、ステップＳ７０２で得られたｎｅｗλと、ステップＳ１０２において取得した顧客ｉの各々の顧客契約情報と、ステップＳ２０６において取得した顧客ｉの各々の顧客契約パラメータｈｉとに基づいて、非定常要因ｙのｎｅｗλの尤度及びλｙの尤度を、下記（６）式により算出し、（（ｎｅｗλの尤度）／（λｙの尤度））であるν´の値を算出する。
【００９５】
【数５】
2015011692.tif 000007
【００９６】
次に、ステップＳ７０６では、一様乱数Ａ（０＜Ａ＜１）を発生させる。
【００９７】
次に、ステップＳ７０８では、ステップＳ７０４で取得したν´の値がステップＳ７０６で取得した一様乱数Ａの値よりも大きいか否かの判定を行う。ν´の値が一様乱数Ａよりも大きい場合には、ステップＳ７１０へ移行し、ν´の値が一様乱数Ａ以下の場合には、ステップＳ７１２へ移行する。
【００９８】
ステップＳ７１０では、指数パラメータλｙの値をステップＳ７０２において取得したｎｅｗλの値に更新する。
【００９９】
次に、ステップＳ７１２では、変数ｙの値がＹの値と一致するか否かの判定を行う。変数ｙの値がＹの値と一致する場合には、指数パラメータ更新処理ルーチンを終了し、変数ｙの値がＹの値と一致しない場合には、ステップＳ７１４へ移行する。
【０１００】
ステップＳ７１４では、変数ｙの値をｙ＋１に設定し、ステップＳ７０２へ移行し、ステップＳ７０２〜ステップＳ７１２の処理を繰り返す。
【０１０１】
図８のステップＳ２２０では、変数ｇの値がバーンインサイクル（burnin）の値よりも大きいか否かの判定を行う。変数ｇの値がburninの値よりも大きい場合には、ステップＳ２２４へ移行し、変数ｇの値がburninの値以下の場合には、ステップＳ２２２へ移行する。
【０１０２】
ステップＳ２２２では、変数ｇの値をｇ＋１に設定し、ステップＳ２０４へ移行する。
【０１０３】
ステップＳ２２４では、ステップＳ２０６において更新した顧客ｉの各々の顧客契約パラメータｈｉ、ステップＳ２１２において更新したＸ＋Ｙ個の要因ｋの各々の重みパラメータｗｋ、ステップＳ２１４において更新したＸ個の定常要因ｘの各々のワイブルスケールパラメータθｘ、ステップＳ２１６において更新したＸ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及びステップＳ２１８において更新したＹ個の非定常要因ｙの各々の指数パラメータλｙの各々の値を、ｇサイクル目のパラメータとして各々記憶する。
【０１０４】
次に、ステップＳ２２６では、変数ｇの値がサイクル数Ｇの値以上か否かの判定を行う。変数ｇの値がサイクル数Ｇの値以上の場合、ステップＳ２２８へ移行し、変数ｇの値がサイクル数Ｇの値未満の場合、ステップＳ２２２へ移行する。
【０１０５】
次に、ステップＳ２２８では、ステップＳ２２４において記憶された、サイクル毎に更新された顧客ｉの各々の顧客契約パラメータｈｉ、Ｘ＋Ｙ個の要因ｋの各々の重みパラメータｗｋ、Ｘ個の定常要因ｘの各々のワイブルスケールパラメータθｘ、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及びＹ個の非定常要因ｙの各々の指数パラメータλｙの各々の値に基づいて、各パラメータの平均値を各パラメータの値と推定し、パラメータ推定処理ルーチンを終了する。
【０１０６】
上記ステップＳ１０６は、図１４に示す結果生成処理ルーチンによって実現される。
【０１０７】
図１４のステップＳ８００では、ステップＳ１０４において推定したパラメータの各々を受け付ける。
【０１０８】
次に、ステップＳ８０２では、ステップＳ１００において取得した仮想発生情報が１つ以上か否かの判定を行う。仮想発生情報が１つ以上であればステップＳ８０４へ移行し、仮想発生情報を取得していない場合には、ステップＳ８１０へ移行する。
【０１０９】
次に、ステップＳ８０４では、ステップＳ１００において取得した仮想発生情報の各々について、仮想発生情報の各々に対応する非定常要因を特定し、ステップＳ８００において取得した当該非定常要因の指数パラメータλｙ、及び当該非定常要因の重みパラメータｗｋを各々取得する。
【０１１０】
次に、ステップＳ８０６では、仮想発生情報の各々について、ステップＳ８０４において取得した指数パラメータλｙ及び重みパラメータｗｋと、発生時期Ｕとを、仮想非定常要因のパラメータとして追加する。ここで、非定常要因の数は、既存の非定常要因の数Ｙと仮想非定常要因の数ｊを併せたＹ´個になる。
【０１１１】
次に、ステップＳ８０８では、重みパラメータｗＸ＋Ｙ´を正規化する。
【０１１２】
次に、ステップＳ８１０では、ステップＳ１００において取得した観測終了時である現在の顧客ｉの継続者数の総数（ｃ＝０の顧客数）をカウントし、Ｎ´とする。なおＮ´に属する顧客をｉ´とする（０≦ｉ´≦Ｎ´）。
【０１１３】
次に、ステップＳ８１２では、変数ｉ´の値を１と設定する。
【０１１４】
次に、ステップＳ８１４では、変数Ｗの値を１と設定する。
【０１１５】
次に、ステップＳ８１６では、変数ｘの値を１と設定する。
【０１１６】
次に、ステップＳ８１８では、顧客ｉ´のＷ週目にｘ番目の定常要因により離脱する離脱確率を下記（７）式に基づいて算出する。
【０１１７】
【数６】
2015011692.tif 000008
【０１１８】
次に、ステップＳ８２０では、変数ｘの値がＸの値と一致するか否かの判定を行う。変数ｘの値がＸと一致する場合には、ステップＳ８２４へ移行し、変数ｘの値がＸの値と一致しない場合には、ステップＳ８２２へ移行する。
【０１１９】
ステップＳ８２２では、変数ｘの値をｘ＋１に設定してステップＳ８１８へ移行し、ステップＳ８１８〜ステップＳ８２０の処理を繰り返す。
【０１２０】
次に、ステップＳ８２４では、変数ｙ´の値を１に設定する。
【０１２１】
次に、ステップＳ８２６では、顧客ｉ´のＷ週目にｙ´番目の非定常要因により離脱する離脱確率を下記（８）式に基づいて算出する。
【０１２２】
【数７】
2015011692.tif 000009
【０１２３】
次に、ステップＳ８２８では、変数ｙ´の値がＹ´の値と一致しているか否かを判定する。変数ｙ´の値がＹ´の値と一致している場合には、ステップＳ８３２へ移行し、変数ｙ´の値がＹ´の値と一致していない場合には、ステップＳ８３０へ移行する。
【０１２４】
ステップＳ８３０では、変数ｙ´の値をｙ´＋１に設定する。
【０１２５】
次に、ステップＳ８３２では、変数Ｗの値が、算出対象の最終週を表す値（例えば、非定常要因の発生時期Ｔ又はＵの内、最も遅い時期の値に２０を加算した値）と一致しているか否か判定する。変数Ｗが算出対象の最終週の値と一致している場合には、ステップＳ８３６へ移行し、変数Ｗが、算出対象の最終週の値と一致していない場合には、ステップＳ８３４へ移行する。
【０１２６】
ステップＳ８３４では、変数Ｗの値をＷ＋１に設定し、ステップＳ８１６へ移行する。
【０１２７】
次に、ステップＳ８３６では、変数ｉ´の値がＮ´の値と一致するか否かを判定する。変数ｉ´の値がＮ´の値と一致する場合には、ステップＳ８４０へ移行し、変数ｉ´の値がＮの値と一致しない場合には、ステップＳ８３８へ移行する。
【０１２８】
ステップＳ８３８では、変数ｉ´の値をｉ´＋１に設定し、ステップＳ８１４へ移行する。
【０１２９】
次に、ステップＳ８４０では、ステップＳ８１８で算出した各週の顧客ｉ´のＸ個の定常要因ｘの各々において離脱する確率の各々と、ステップＳ８２６で算出した各週の顧客ｉ´のＹ個の非定常要因ｙ´の各々において離脱する確率の各々とに基づいて、各週及び各要因について期待離脱者数を集計する。
【０１３０】
次に、ステップＳ８４２では、ステップＳ８４０において集計した各週及び各要因についての期待離脱者数に基づいて、結果を図表に表して処理を終了する。
【０１３１】
以上説明したように、第１の実施の形態に係る生存分析装置によれば、精度良く離脱確率の分布を推定することができる。
【０１３２】
また、生存時間上の定常要因離脱と絶対時間上の非定常離脱を分離することで、別の絶対時間上又は別の地域に非定常離脱が発生した場合の離脱者数を推定できる。
【０１３３】
また、絶対時刻に発生したイベントを分離した分析が可能になり、将来の時刻で過去のイベントを模擬し、仮想的なイベントを発生させ、その仮想イベントによって事象が停止となる確率の時間的推移を予測できる。
【０１３４】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１３５】
第１の実施の形態に係る生存分析装置においては、顧客契約情報におけるｃは現在の状態を表す場合について説明したが、これに限定されるわけではない。例えば、ｃは現在の状態ではなく、観測終了時又は別の時期の状態を表してもよい。ただし、ｃを観測終了時又は別の時期の状態を表すことにした場合、ｃの観測時に伴って観測継続期間ｓの値も再計算する。具体的には、ｃが現在（６月２０日）の状態を表し、ｃ＝１、ｓ＝３０（日単位）となる場合、ｃの観測時を６月１９日に変更すると、ｃ＝１、ｓ＝２９（日単位）となる。
【０１３６】
第１の実施の形態に係る生存分析装置においては、ｔをデータ取得日から初回契約までの経過期間とし、Ｔをデータ取得日から非定常要因発生までの期間とし、Ｕをデータ取得日から仮想非定常要因発生までの期間とする場合について説明したが、これに限定されるわけではない。例えば、ｔ、Ｔ、及びＵに係るデータ取得日を同一の任意の時点としてもよい。
【０１３７】
第１の実施の形態に係る生存分析装置においては、契約している顧客の離脱確率を推定する場合について説明したが、これに限定されるわけではない。例えば、個体に発生する現象を、製品の出荷とし、現象の離脱を、製品の故障として、製品の故障率の分布を推定してもよい。また、個体に発生する現象を生物の誕生とし、現象の離脱を生物の死亡として、生物の死亡率の分布を推定してもよい。また、個体に発生する現象を、疾患の発症とし、現象の離脱を、生物の死亡として、疾患発症からの生物の死亡率の分布を推定してもよい。このように、他の経時的現象発生に適用してもよい。
【０１３８】
第１の実施の形態に係る生存分析装置においては、ｓ、ｔ、Ｔ、及びＵの単位については週数とする場合について説明したが、これに限定されるわけではない。例えば、日数、年数等、他の期を指定してもよい。なお、他の期を指定する場合には、ｓ、ｔ、Ｔ，及びＵについては同一の単位を使用する必要がある。
【０１３９】
第１の実施の形態に係る生存分析装置においては、ワイブル形状パラメータを推定する際にランダムウォークを使用する場合について説明したが、これに限定されるわけではない。例えば、スライスサンプリング等の他の手法においてワイブル形状パラメータを推定してもよい。また、分散値は１以外を使用してもよい。
【０１４０】
第１の実施の形態に係る生存分析装置においては、指数パラメータを推定する際にランダムウォークを使用する場合について説明したが、これに限定されるわけではない。例えば、スライスサンプリング等の他の手法において指数パラメータを推定してもよい。また、分散値は１以外を使用してもよい。
【０１４１】
第１の実施の形態に係る生存分析装置においては、定常離脱をワイブル分布に従うものとし、非定常離脱を位置パラメータ付き指数分布に従うものと仮定する場合について説明したが、これに限定されるわけではない。例えば、定常離脱又は非定常離脱を別の分布に従うものと仮定してもよい。
【０１４２】
第１の実施の形態に係る生存分析装置においては、定常要因又は非定常要因の要因毎に離脱確率を算出する場合について説明したがこれに限定されるものではない。例えば、上記（１）式に基づいて、各顧客の各時期ｖの離脱確率を算出し、離脱確率の算出結果を集計して各時期の期待離脱者数を算出してもよい。
【０１４３】
次に、第２の実施の形態について説明する。なお、第１の実施の形態と同様の構成及び作用となる部分については、同一符号を付して説明を省略する。
【０１４４】
第２の実施の形態では、定常分布数を自動で決定している点が第１の実施の形態と異なっている。
【０１４５】
＜第２の実施の形態に係る生存分析装置の作用＞
次に、本発明の実施の形態に係る生存分析装置１００の作用について説明する。図７に示す生存分析処理ルーチンを実行する前に、入力部１０により全顧客の各々の顧客契約情報を受け付け記憶部２２に記憶する。そして、入力部１０により分析対象となる顧客、発生情報、及び仮想発生情報を受け付けると、生存分析装置１００は、図７に示す生存分析処理ルーチンを実行する。
【０１４６】
図７のステップＳ１０４は、図１５に示すパラメータ推定処理ルーチンによって実現される。
【０１４７】
図１５のステップＳ１０００では、定常分布数Ｘの値を１に設定し、非定常分布数の値を非定常要因数のＹに設定する。
【０１４８】
次に、ステップＳ１００２では、ｍａｘＸの値を１０に、ＭＣＭＣサイクルのパラメータであるサイクル数のパラメータＧの値を１５０００に、バーンインサイクル（burnin）の値を１００００の値に設定する。
【０１４９】
次に、ステップＳ１００４では、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及びワイブルスケールパラメータθｘ、Ｘ＋Ｙ個の要因ｋの各々の重みパラメータｗｋ、顧客ｉの各々の顧客所属パラメータｈｉ、Ｙ個の非定常要因ｙの各々の指数パラメータλｙ、ワイブル形状パラメータの事前パラメータ（αａ，βａ）、ワイブルスケールパラメータの事前パラメータ（αθ，βθ）、及び重みパラメータの事前ディリクレパラメータφの各々に適切な初期値を各々設定する。
【０１５０】
次に、ステップＳ１００６では、ステップＳ２２４において記憶された、サイクル毎に更新された顧客ｉの各々の顧客契約パラメータｈｉ、Ｘ＋Ｙ個の要因ｋの各々の重みパラメータｗｋ、Ｘ個の定常要因ｘの各々のワイブルスケールパラメータθｘ、Ｘ個の定常要因ｘの各々のワイブル形状パラメータａｘ、及びＹ個の非定常要因ｙの各々の指数パラメータλｙの各々の値に基づいて、各パラメータの平均値を定常分布数がＸの場合における各パラメータの値と推定する。
【０１５１】
次に、ステップＳ１００８では、ステップＳ１００６において取得した定常分布数がＸの場合におけるパラメータの各々の値に基づいて、ＡＩＣ値を算出し、定常分布数がＸの場合におけるＡＩＣ値として記憶する。
【０１５２】
次に、ステップＳ１０１０では、Ｘの値がｍａｘＸの値と一致するか否かの判定を行う。Ｘの値がｍａｘＸの値と一致する場合には、ステップＳ１０１４へ移行し、Ｘの値がｍａｘＸの値と一致しない場合には、ステップＳ１０１２へ移行する。
【０１５３】
ステップＳ１０１２では、Ｘの値をＸ＋１に設定し、ステップＳ１００４へ移行する。
【０１５４】
次に、ステップＳ１０１４では、ステップＳ１００８において取得した定常分布数のＸ毎のＡＩＣ値の中で最も値が大きいものを取得する。そして、取得したＡＩＣ値に対応する定常分布数Ｘを定常分布数Ｘに決定し、当該定常分布数Ｘにおける各パラメータの値を各パラメータの値として推定する。
【０１５５】
以上説明したように、第２の実施の形態に係る生存分析装置によれば、精度良く離脱確率の分布を推定することができる。
【符号の説明】
【０１５６】
１０   入力部
２０   演算部
２２   記憶部
２４   データ取得部
２６   パラメータ初期値設定部
３０   モデル推定部
３２   パラメータ追加部
３４   予測部
５０   出力部
１００ 生存分析装置
【図１】
2015011692.tif 000010
【図２】
2015011692.tif 000011
【図３】
2015011692.tif 000012
【図４】
2015011692.tif 000013
【図５】
2015011692.tif 000014
【図６】
2015011692.tif 000015
【図７】
2015011692.tif 000016
【図８】
2015011692.tif 000017
【図９】
2015011692.tif 000018
【図１０】
2015011692.tif 000019
【図１１】
2015011692.tif 000020
【図１２】
2015011692.tif 000021
【図１３】
2015011692.tif 000022
【図１４】
2015011692.tif 000023
【図１５】
2015011692.tif 000024