JP7456514B2 - パラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラム - Google Patents

Description

本発明は、秘密計算技術において、秘匿性を保ったままｃｏｘ比例ハザードモデルのパラメータ推定を実現する技術に関連するものである。

ｃｏｘの比例ハザードモデルを用いた回帰分析であるｃｏｘ比例ハザード回帰は、生存時間解析でよく用いられる分析手法である（非特許文献１）。市販の統計解析ソフトウェアやソフトウェア言語のパッケージにおいて、平文でｃｏｘ比例ハザード回帰を行うことができる。

また、暗号化された数値を復元すること無く特定の演算結果を得る方法として、秘密計算と呼ばれる方法が知られている。一例として、複数の秘密計算装置に数値の断片を分散させるという暗号化を行い、複数の秘密計算装置が協調計算を行うことにより、数値を復元すること無く、加減算、定数加算、乗算、定数倍、論理演算（否定、論理積、論理和、排他的論理和）、データ形式変換（整数と二進数）等の結果を複数の秘密計算装置に分散された状態として得ることができる。

D. R. Cox. Regression Models and Life-Tables. Journal of the Royal Statistical Society. Series B (Methodological), Vol. 34, No. 2. (1972), pp. 187-220.

ｃｏｘ比例ハザードモデルのパラメータ推定を行う際、平文では死亡や打ち切りが発生した時点毎に計算を繰り返す。しかし、秘密計算でも同様の方法を用いる場合、秘匿しておくべき時点の値を復号する必要があるという課題がある。また、ｃｏｘ比例ハザードモデルのパラメータ推定の計算では、指数、除算、ｇｒｏｕｐ－ｂｙ ｓｕｍといった秘密計算での処理コストが大きい処理が多く含まれるため、秘密計算で効率良く計算するのが難しいという課題もある。

本発明は上記の点に鑑みてなされたものであり、時点の値を復号することなく、効率良くｃｏｘ比例ハザードモデルのパラメータ推定を行うための技術を提供することを目的とする。

開示の技術によれば、ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置であって、
イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
前記演算部により推定されたパラメータを出力する出力部と
を備えるパラメータ推定装置が提供される。

開示の技術によれば、時点の値を復号することなく、効率良くｃｏｘ比例ハザードモデルのパラメータ推定を行うための技術が提供される。

本発明の実施の形態におけるパラメータ推定装置の構成図である。 装置のハードウェア構成例を示す図である。 置換表を説明するための図である。 フラグを説明するための図である。 データベースの例を示す図である。 パラメータ推定装置の処理手順の例を説明するための図である。 パラメータ推定装置の処理手順の例を説明するための図である。

以下、図面を参照して本発明の実施の形態（本実施の形態）を説明する。以下で説明する実施の形態は一例に過ぎず、本発明が適用される実施の形態は、以下の実施の形態に限られるわけではない。

（装置構成例）
図１に、本実施の形態におけるパラメータ推定装置１００の構成図を示す。図１に示すように、本実施の形態におけるパラメータ推定装置１００は、入力部１１０、演算部１２０、出力部１３０、及びデータ格納部１４０を有する。

パラメータ推定装置１００は、１つの装置（コンピュータ）で構成されてもよいし、複数のコンピュータからなるシステムとして構成されてもよい。このシステムをパラメータ推定システムと呼んでもよい。パラメータ推定システムにおいて、例えば、演算部１２０とデータ格納部１４０が別々のサーバであってもよい。

パラメータ推定装置１００の入力部１１０に、観測により得られたデータから秘匿化されたデータが入力される。特に断らなければ、パラメータ推定装置１００により扱われるデータは秘匿化されたデータであり、計算は秘密計算でなされるものとする。

入力されたデータはデータベースとしてデータ格納部１４０に格納される。演算部１２０は、データ格納部１４０のデータベースから読み出したスカラー、ベクトル、行列等のデータに対して後述する処理を行うことで、ｃｏｘ比例ハザードモデルのパラメータ推定を行う。出力部１３０は、演算部１２０で推定されたパラメータを出力する。

なお、演算部１２０で計算されたパラメータがデータ格納部１４０に格納され、外部からのアクセスに応じて、出力部１３０から出力されてもよい。演算部１２０での処理内容の詳細については後述する。

（ハードウェア構成例）
本実施の形態におけるパラメータ推定装置１００は、例えば、コンピュータに、本実施の形態で説明する処理内容を記述したプログラムを実行させることにより実現可能である。なお、この「コンピュータ」は、物理マシンであってもよいし、クラウド上の仮想マシンであってもよい。仮想マシンを使用する場合、ここで説明する「ハードウェア」は仮想的なハードウェアである。

上記プログラムは、コンピュータが読み取り可能な記録媒体（可搬メモリ等）に記録して、保存したり、配布したりすることが可能である。また、上記プログラムをインターネットや電子メール等、ネットワークを通して提供することも可能である。

図２は、上記コンピュータのハードウェア構成例を示す図である。図２のコンピュータは、それぞれバスＢで相互に接続されているドライブ装置１０００、補助記憶装置１００２、メモリ装置１００３、ＣＰＵ１００４、インタフェース装置１００５、表示装置１００６、入力装置１００７、出力装置１００８等を有する。

当該コンピュータでの処理を実現するプログラムは、例えば、ＣＤ－ＲＯＭ又はメモリカード等の記録媒体１００１によって提供される。プログラムを記憶した記録媒体１００１がドライブ装置１０００にセットされると、プログラムが記録媒体１００１からドライブ装置１０００を介して補助記憶装置１００２にインストールされる。但し、プログラムのインストールは必ずしも記録媒体１００１より行う必要はなく、ネットワークを介して他のコンピュータよりダウンロードするようにしてもよい。補助記憶装置１００２は、インストールされたプログラムを格納すると共に、必要なファイルやデータ等を格納する。

メモリ装置１００３は、プログラムの起動指示があった場合に、補助記憶装置１００２からプログラムを読み出して格納する。ＣＰＵ１００４は、メモリ装置１００３に格納されたプログラムに従って、パラメータ推定装置１００に係る機能を実現する。インタフェース装置１００５は、ネットワークに接続するためのインタフェースとして用いられる。表示装置１００６はプログラムによるＧＵＩ（Ｇｒａｐｈｉｃａｌ Ｕｓｅｒ Ｉｎｔｅｒｆａｃｅ）等を表示する。入力装置１００７はキーボード及びマウス、ボタン、又はタッチパネル等で構成され、様々な操作指示を入力させるために用いられる。出力装置１００８は演算結果を出力する。

（準備）
パラメータ推定装置１００の動作例を説明する前に、その準備として、記法、ｃｏｘ比例ハザード回帰、及び秘密計算等について説明する。

＜記法＞
ａをｂで定義することをａ：＝ｂと記載し、ベクトルを→ａ：＝（ａ_０，...，ａ_ｎ－１）と記載し、特筆しない限りＡのような大文字は行列を表し、その転置行列はＡ^Ｔと記載する。なお、本明細書のテキストにおいて、記載の便宜上、文字の頭の上に置かれるべきベクトルの記号（→）を、「→ａ」のように、文字の前に記載している。また、図面に記載したアルゴリズムにおいて、３階のテンソルは斜体の文字で記載している。本明細書のテキストにおいては、記載の便宜上、３階のテンソルの文字を、例えば^{ｔｅｎｓｏｒ}ＺＺ´のように、その文字の左上に^{ｔｅｎｓｏｒ}を付けて表す。

加減乗算において入力がベクトル→ａもしくは行列Ａと、スカラーｂの場合、→ａ、Ａの全ての要素に対してｂとの演算を行うものとする。また、特に記載が無いベクトルは列ベクトルである。横ベクトルの場合は^ｔ→ａのように左上にｔを付けることで区別する。

＜ｃｏｘ比例ハザード回帰＞
ｃｏｘ比例ハザードモデルは式（１）で表されるモデルである（非特許文献１）。

式（１）において、ｔ、β、ｚはそれぞれ時間、重み、特徴量を表し、λ_０（ｔ）、ｅｘｐ（→β^Ｔ→ｚ）はそれぞれベースラインハザード関数、相対危険度関数（ハザード）と呼ばれる。ｃｏｘ比例ハザード回帰では重みのパラメータβを推定する。重みは、式（２）に示す部分尤度を用いて最尤推定を行うことで計算する。

式（２）において、Ｄは死亡が観測された時点の数であり、→ｚ_ｉは時点ｉに死亡した患者の特徴量を表す。なお、本実施の形態では、目的となるイベントとして、死亡を用いているが、これは例である。例えば、転倒、歩行不可能、疾患発症、入院等を目的のイベントとしてもよい。打ち切りをイベントの１つと解釈してもよい。

式（２）におけるＲ_ｉは、時点ｉの直前まで打ち切りも死亡も発生していない患者の集合であり、リスクセットと呼ばれる。なお、「打ち切り」とは、観測ができなくなり、それ以降、死亡の発生の有無が不明になることである。

従って、式（２）の部分尤度関数は、時点毎に（死亡した患者のハザード）／（リスクセットのハザードの総和）を計算し、全時点分掛け合わせたものである。この部分尤度は同じ時刻に複数の打ち切りや死亡が発生していない（タイデータが無い）という仮定を置いているため、タイデータがあることの多い実データでは、式（３）に示すＢｒｅｓｌｏｗ法が良く用いられる。

式（３）は、基本的には式（２）のｃｏｘの部分尤度と同じであるが、分母がｄ_ｉ乗（ｄ_ｉは時点ｉの死亡患者数）されている点と、分子に→ｚ_ｉの代わりに、時点ｉの死亡患者の特徴量の総和ｓ_ｉを用いる点が異なる。ｄ_ｉ＝１の場合は式（２）と一致するため、式（３）は、式（２）を一般化した形と考えることができる。以下の説明では、式（３）を前提としている。

→βの最尤推定量を求める方法としては、Ｎｅｗｔｏｎ法等が一般的である。本実施の形態ではＮｅｗｔｏｎ法を用いている。Ｎｅｗｔｏｎ法では、式（３）を対数尤度関数へと変形した後、その対数尤度関数の１階微分（勾配）と２階微分（ヘシアン）を用いて計算する。対数尤度関数ｌ（→β）、その１階微分Ｕ（→β）、及び２階微分Ｉ（→β）をそれぞれ式（４）、式（５）、式（６）に示す。

Ｎｅｗｔｏｎ法では、式（５）と式（６）を用いて、以下の式（７）を反復して→βの最尤推定値を求める。式（７）は、およそ５回ほどの反復で収束する。

＜秘密計算＞
ある値ａを暗号化や秘密分散等により秘匿化した値をａの暗号文あるいは秘匿値と呼び、［ａ］と記述する。ａが秘密分散により秘匿化された場合は、［ａ］により各秘密計算装置が持つ秘密分散の断片の集合を参照する。なお、暗号文であることを示す括弧"［"、"］"について、図面及び明細書に挿入された数式での括弧の書式と少し異なるが、明細書テキストにおいては、記載の便宜上、"［"、"］"を使用している。

次に、秘密計算における各種の演算を説明する。

■四則演算
２つの暗号文［ａ］、［ｂ］の加算、減算、乗算は、それぞれ暗号文［ａ＋ｂ］、［ａ‐ｂ］、［ａ×ｂ］を計算する処理である。これらの演算をそれぞれ、［ａ］＋［ｂ］、［ａ］‐［ｂ］、［ａ］×［ｂ］と記載する。

また、暗号文［ａ］を平文ｂで割る処理は、［ａ］／ｂのような記法とする。入力がベクトルや行列で、要素毎にこれらの処理を行う場合も同様に［→ａ］＋［→ｂ］、［Ａ］＋［Ｂ］のような記法とする。

加減乗算において入力が行列Ａと列ベクトル→ｂの場合は、行列の各列ベクトルに対して→ｂとの要素毎の演算を実施し、入力が行列Ａと行ベクトル^ｔ→ｂの場合は、行列の各行ベクトルに対して^ｔ→ｂとの要素毎の演算を実施するものとする。

■総和
ベクトル［→ａ］の要素の総和を求める処理をｓｕｍ（［→ａ］）と記述する。またｓｕｍ（［Ａ］）のようにｍ×ｎの行列がｓｕｍの入力の場合は、列方向の総和を計算し、長さｎの行ベクトル［^ｔ→ｃ］を出力するものとする。

■ｐｒｅｆｉｘ ｓｕｂ
ベクトル［→ａ］：＝（［ａ_１］，［ａ_２］，・・・，［ａ_ｎ］）とスカラー［ｂ］から（［ｂ］，［ｂ］－［ａ_１］，［ｂ］－（［ａ_１］＋［ａ_２］），・・・，［ｂ］－Σ［→ａ］）となるベクトルを計算することを、ｐｒｅｆｉｘＳｕｂ（［→ａ］，［ｂ］）と記載する。

■逆数
暗号文［ａ］の逆数１／［ａ］を計算することを、［ｃ］←ｒｅｃｉｐｒｏｃａｌ（［ａ］）のように記載する。入力がベクトルの場合も同じ記法とする。

■指数
暗号文［ａ］を入力とし、ネイピア数ｅの［ａ］乗を計算することを［ｃ］←ｅｘｐ（［ａ］）のように記載する。入力がベクトルの場合も同じ記法とする。

■Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎ
Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎは、Ｇｒｏｕｐ－ｂｙ ｓｕｍやＧｒｏｕｐ－ｂｙ ｃｏｕｎｔといった様々なＧｒｏｕｐ－ｂｙ演算で共通的に用いることができる中間データを生成する処理である。中間データは置換表［→π］と、キーの値の境目かどうかを表すフラグ［→ｅ］を含み、これらを使いまわすことで、同じキーを用いた様々なＧｒｏｕｐ－ｂｙ演算を効率良く行うことができる。

図３、図４を参照して、Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎの演算で得られる置換表→πと、フラグ→ｅの例を説明する。ここでは、説明の便宜上、平文の処理として説明する。

置換表→πは、Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎに入力されたベクトルにおける各要素が、何番目の要素に移動すれば、当該ベクトルの要素がソートされるかを表すベクトルである。

図３において、ベクトル→ｋをＧｒｏｕｐ－ｂｙ ｃｏｍｍｏｎに入力すると、当該ベクトルがソートされ、図３に示すとおりの置換表→πが得られる。例えば、置換表→πの最初の要素である３は、入力ベクトル→ｋの最初の要素が、３番目の要素になればソートできることを示している。

フラグ→ｅは、Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎに入力されたベクトルのソート済みのベクトルの各要素に対して、その要素の下の値と比較して同じならば、その要素の位置に０を入れ、異なる場合は１を入れたベクトルである。図４に、入力とフラグの例を示す。図４に示すように、フラグの一番最後については、下の値と比較できないので１を入れる。

キーのベクトル［→ｋ］を入力してＧｒｏｕｐ－ｂｙ ｃｏｍｍｏｎを行うことを式（８）のように記述する。置換表［→π］を用いてベクトル［→ａ］や行列［Ａ］（行数が［→π］の長さと等しい）をソートする処理を式（９）、式（１０）のように記述する。

ソート済みのベクトル［→ａ´］や［Ａ´］と［→ｅ］を用いてＧｒｏｕｐ－ｂｙ ｓｕｍを行うことを式（１１）（１２）のように記述する。［→ａ´］、［Ａ´］などのプライムは、ソート済みであることを表し、以降も同様の記法とする。

ｓｏｒｔ、ｇｒｏｕｐＢｙＳｕｍの入力が行列の場合、処理は列ごとに行われる。また一般的にＧｒｏｕｐ－ｂｙ ｓｕｍを行うと出力のサイズは入力のサイズ以下になるが、本実施の形態では入力と出力のサイズは同じであり、不要な分は末尾を０でパディングする。これによって時点キー属性の数を秘匿することができる。なお、後述する処理の説明において、便宜的にＧｒｏｕｐ－ｂｙ ｓｕｍの結果を「長さが時点数のベクトル」や「時点数×ｎの行列」と記載しているが、実際には不要な部分を０でパディングした「長さがレコード数のベクトル」や「レコード数×ｎの行列」を扱っている。

例えば、図４に示すフラグ→ｅを（０，１，０，０，１，１）と記載するとして、仮に→ａ´が（２，１，３，５，１，２）であるとすると、→ａ´と→ｅを入力とするＧｒｏｕｐ－ｂｙ ｓｕｍは、（２＋１，３＋５＋１，２，０，０，０）＝（３，９，２，０，０，０）になる。

（パラメータ推定装置１００の動作について）
以下、パラメータ推定装置１００の動作例を説明する。パラメータ推定装置１００の演算部１２０は、データ格納部１４０（データベース）に格納された暗号文のデータを読み出し、秘密計算により、前述した式（５）、式（６）、式（７）を計算して、ｃｏｘ比例ハザード回帰のパラメータ推定を行う。以下では、まず、特徴的な動作について説明する。

＜全レコードをまとめて処理する＞
式（５）、式（６）をそのまま実装した場合、時点毎に反復計算していき、計算結果を順番に足していくような処理になる。秘密計算ｃｏｘ比例ハザード回帰を計算する場合には時点数を秘匿するため、前述のｇｒｏｕｐＢｙＣｏｍｍｏｎ、ｇｒｏｕｐＢｙＳｕｍを用いる。

ｇｒｏｕｐＢｙＣｏｍｍｏｎ、ｇｒｏｕｐＢｙＳｕｍでは時点の値を秘匿したまま時点毎の集計を行い、不要な部分は０でパディングされるため、時点数に関する情報が漏れない。

また、時点数毎の反復処理ではなく全レコードをまとめて処理することにより、秘密計算での処理コストが大きい演算の回数を減らすことができて処理効率も良い。つまり、例えば、時点１でΣ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）を計算し、時点２でΣ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）を計算し、...といった計算ではなく、全時点分のΣ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）を一度に計算する。より具体的には、各時点の値がスカラーの場合は長さが時点数のベクトルとして扱ってまとめて計算し、各時点の値が長さｎのベクトルの場合は時点数×ｎの行列として扱って、まとめて計算する。また、各時点の値がｎ×ｎ行列の場合は時点数×ｎ×ｎの３階のテンソルとして扱ってまとめて計算する。

＜処理コストの大きい演算の削減＞
式（５）、式（６）にはｅｘｐや除算が多く含まれ、また、Σ_ｊ∈Ｒｉの処理はＧｒｏｕｐ－ｂｙ ｓｕｍであるため、秘密計算での計算コストが非常に大きい。

本実施の形態に係るパラメータ推定装置１００は、ｃｏｘ比例ハザード回帰におけるｅｘｐ、除算、Ｇｒｏｕｐ－ｂｙ ｓｕｍといったコストの大きい処理を最小限に抑え、効率良く計算する。単純に式（５）、式（６）の通りに計算した場合、Ｎｅｗｔｏｎ法の１反復あたりｅｘｐは７回×時点数、除算は３回×時点数分必要になるが、本実施の形態では下記の通り最小限に抑えた。

・ｅｘｐの計算が１反復あたり１回
・逆数の計算が１反復あたり１回
また、Ｇｒｏｕｐ－ｂｙ ｓｕｍについても、処理をｇｒｏｕｐＢｙＣｏｍｍｏｎ、ｇｒｏｕｐＢｙＳｕｍに分けることで、より効率的に処理を行うこととしている。これらの演算の削減についてより詳細に説明する。

■ｅｘｐの削減
式（５）、式（６）に示すとおり、ｅｘｐの引数が全て→β^Ｔ→ｚ_ｊであるため、一度計算したら後は使い回せば良い。これに加えて、前述の全レコードをまとめて処理することで処理が並列化され、１反復あたり１度で済むようになる。

■除算の削減
除算を逆数の計算＋乗算で行う場合、Σ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）の逆数を、式（５）の第２項と式（６）の第１項で使い回せるため、除算２回ではなく逆数の計算１回＋乗算２回で済む。式（６）の第２項は前述の２つとは除数が異なるが、この項は除算をしなくても求めることができる。式（５）の第２項と式（６）の第２項を見比べてみると、式（５）の第２項のｄ_ｉを除いた部分をＡとしたとき、式（６）の第２項はＡＡ^Ｔで表せるため、積のみで式（６）の第２項は計算できる。これに加えて、前述の全レコードをまとめて処理することで処理が並列化され、１反復あたり１度で済むようになる。

■Ｇｒｏｕｐ－ｂｙ ｓｕｍの削減
式（５）、式（６）の通りに実装するとＧｒｏｕｐ－ｂｙｓｕｍを処理の中で何度も行うことになり、処理効率が低下する。そこで、本実施の形態では、全てキーが同一であることに着目し、最初にキーに対して一度だけ行ったＧｒｏｕｐ－ｂｙｃｏｍｍｏｎで得た［→ｅ］を使い回すこととしている。

本実施の形態のパラメータ推定装置１００では、Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎを活用し、秘密計算ｃｏｘ比例ハザード回帰を効率良く計算する。ｇｒｏｕｐＢｙＳｕｍの計算では、境目を表すフラグ［→ｅ］を用いた集計を行うのみである。

＜詳細な処理内容について＞
次に、パラメータ推定装置１００が実行する詳細な処理内容を説明する。ここでは、秘匿化された観測データがデータ格納部１４０においてデータベースとして格納されているものとし、演算部１２０は、秘密計算によりそのデータに対して処理を行うことで、ｃｏｘ比例ハザード回帰のパラメータを推定する。処理動作において、上述した特徴的な動作が行われる。

演算部１２０による処理の対象となるデータベースのイメージを図５に示す。図５は、説明の便宜上、データが平文で示されており、また、時点の昇順にソート済の状態を示している。

図５に示すように、データベースには、患者数ｍの患者（観測対象）毎に、ｎ個の特徴量と、観測の時点と、その時点での状態（死亡＝１、打ち切り＝０）が格納されている。ｍはレコード数でもあり、Ｄは時点数である。例えば、時刻ベクトル→ｔ＝（１，１，１，２）であれば、時点数は２である。

図５の例では、１患者に対し、死亡があった場合に、その患者、その患者の特徴量、その時点、その状態がデータベースに記録される。また、その時点で打ち切りされている患者についても、その患者、その患者の特徴量、その時点、その状態がデータベースに記録される。

同時点で、複数患者に対して死亡又は打ち切りが観測される場合がある。そのため、時点数Ｄ≦患者数ｍである。

演算部１２０は、データベースからデータを読み出すことで、患者全員分の特徴量をｍ×ｎの行列Ｚとし、時点を時刻ベクトル→ｔとし、患者全員分の状態を状態ベクトル→ｃとして保持する。

Ｚ、→ｔ、→ｃはいずれも初期状態においてはソートされていない。前述したように、Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎにより、最初に時刻ベクトル→ｔをソートして置換表→πを作り、それを使い回すことで、特徴量Ｚや状態ベクトル→ｃの→ｔをキーとするソートは、置換表→πに基づいた並べ替えのみで済む。すなわち、通常のソートを行うより低コストで済む。

＜処理手順＞
パラメータ推定装置１００の演算部１２０は、上記のデータベースのデータに対して、図６、図７に示すアルゴリズムの手順に従って、パラメータ推定を行う。図６、図７には、説明のために行番号を付している。以下では、処理部分の行番号をステップ番号とみなして説明を行う。

図６のアルゴリズム１のステップ３において、演算部１２０は、ｎ次の重みベクトルである［→β］を［０］で初期化する。ステップ４において、Ｇｒｏｕｐ－ｂｙ ｃｏｍｍｏｎにより、→ｔをソートして置換表→πを作成するとともに、フラグ→ｅを作成する。

ステップ５、６において、演算部１２０は、［Ｚ］、［→ｃ］をそれぞれソートして、［Ｚ´］、［→ｃ´］を作成する。

演算部１２０は、ステップ８において、［Ｚ´］から死亡例の特徴量以外を０とした［Ｚ´_ｄｅａｄ］を作成し、ステップ９において、ｇｒｏｕｐＢｙＳｕｍにより時点毎の死亡例の特徴量の総和からなる［Ｓ］を作成する。ステップ１１において、時点毎の死亡数［→ｄ］を算出する。ステップ１３において、ｍ×ｎ×ｎテンソルである［^{ｔｅｎｓｏｒ}ＺＺ´］を作成し、ステップ１５～１７において［→β］の更新を行う。

ステップ１６のｃａｌｃＧＨの処理については、図７を参照して説明する。ステップ３～６において、演算部１２０は、→ｚｅｘｐ（→β^Ｔ→ｚ）に対応する［Ｗ´］と、→ｚ→ｚ^Ｔｅｘｐ（→β^Ｔ→ｚ）に対応する［^{ｔｅｎｓｏｒ}Ｘ´］等を算出する。ステップ４で計算したｅｘｐ（→β^Ｔ→ｚ_ｊ）に対応する［^ｔ→ｖ´］が以降の計算で使い回されている。

ステップ８～１０において、演算部１２０は、各時点のΣ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）に対応する［→ｖ_ｐｓｕｂ］を算出する。時点数長のベクトルである［→ｖ_ｐｓｕｂ］の各要素が、スカラー値であるΣ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）になっている。つまり、ここでの計算は、時点数毎の反復処理でなく、全レコードをまとめて処理する計算である。下記で説明する［Ｗ_ｐｓｕｂ］、［^{ｔｅｎｓｏｒ}Ｘ_ｐｓｕｂ］の算出においても同様に全レコードをまとめて処理している。

ステップ１２～１４において、演算部１２０は、各時点のΣ_ｊ∈Ｒｉ→ｚ_ｊｅｘｐ（→β^Ｔ→ｚ_ｊ）に対応する［Ｗ_ｐｓｕｂ］を算出する。ステップ１６～１８において、演算部１２０は、各時点のΣ_ｊ∈Ｒｉ→ｚ_ｊ→ｚ_ｊ ^Ｔｅｘｐ（→β^Ｔ→ｚ_ｊ）に対応する［^{ｔｅｎｓｏｒ}Ｘ_ｐｓｕｂ］を算出する。

ステップ２０において、演算部１２０は、Σ_ｊ∈Ｒｉｅｘｐ（→β^Ｔ→ｚ_ｊ）の逆数に対応する［→ｙ］を算出する。逆数の計算はこの部分だけである。

ステップ２２～２５において、演算部１２０は、勾配である式（５）を計算する。［→ｙ］、［→ｄ］は長さが時点数のベクトルであり、［Ｇ］、［Ｗ］、［Ｓ］は時点数×特徴量数の行列であり、計算結果の［^ｔ→ｇ］は長さが特徴量数の横ベクトルである。ステップ２５のｓｕｍにより、全時点の総和を計算する。

ステップ２７～３０において、演算部１２０は、へシアンである式（６）を計算する。ステップ２９の［Ｇ_ｔｍｐ］［Ｇ_ｔｍｐ］^Ｔが、前述したＡＡ^Ｔに対応している。

（実施の形態の効果）
以上説明した本実施の形態に係る技術により、時点の数を復号することなく、効率良くｃｏｘ比例ハザードモデルのパラメータ推定を行うことが可能となる。

すなわち、本実施の形態に係る技術により、従来の平文での処理のような反復処理を行わず全データを一度に処理することで、時点数を復号せずに計算できるようになる。また、反復処理を減らしたことにより、除算、指数、ｇｒｏｕｐ－ｂｙ ｓｕｍと言った秘密計算での処理コストが大きい処理が並列化され、データや時点数を秘匿したまま効率良くｃｏｘ比例ハザード回帰のパラメータ推定が行えるようになる。

（実施の形態のまとめ）
本明細書には、少なくとも下記の各項に記載したパラメータ推定装置、パラメータ推定システム、パラメータ推定方法、及びプログラムが記載されている。
（第１項）
ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置であって、
イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
前記演算部により推定されたパラメータを出力する出力部と
を備えるパラメータ推定装置。
（第２項）
前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数のｅｘｐの計算を、１反復あたり１回のｅｘｐの計算と、その計算結果を利用した演算により実行する
第１項に記載のパラメータ推定装置。
（第３項）
前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数の逆数の計算を、１反復あたり１回の逆数の計算と、その計算結果を利用した演算により実行する
第１項又は第２項に記載のパラメータ推定装置。
（第４項）
前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における時点毎の計算を、ベクトル又は行列又はテンソルを利用して全時点分についてまとめて実行する
第１項ないし第３項のうちいずれか１項に記載のパラメータ推定装置。
（第５項）
ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定システムであって、
イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
前記演算部により推定されたパラメータを出力する出力部と
を備えるパラメータ推定システム。
（第６項）
ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置により実行されるパラメータ推定方法であって、
イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算ステップと、
前記演算ステップにより推定されたパラメータを出力する出力ステップと
を備えるパラメータ推定方法。
（第７項）
コンピュータを、第１項ないし第４項のうちいずれか１項に記載のパラメータ推定装置における各部として機能させるためのプログラム。

以上、本実施の形態について説明したが、本発明はかかる特定の実施形態に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形・変更が可能である。

１００ パラメータ推定装置
１１０ 入力部
１２０ 演算部
１３０ 出力部
１４０ 格納部
１０００ ドライブ装置
１００１ 記録媒体
１００２ 補助記憶装置
１００３ メモリ装置
１００４ ＣＰＵ
１００５ インタフェース装置
１００６ 表示装置
１００７ 入力装置
１００８ 出力装置

Claims (7)

1. ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置であって、
   イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
   前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
   前記演算部により推定されたパラメータを出力する出力部と
   を備えるパラメータ推定装置。
2. 前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数のｅｘｐの計算を、１反復あたり１回のｅｘｐの計算と、その計算結果を利用した演算により実行する
   請求項１に記載のパラメータ推定装置。
3. 前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における複数の逆数の計算を、１反復あたり１回の逆数の計算と、その計算結果を利用した演算により実行する
   請求項１又は２に記載のパラメータ推定装置。
4. 前記演算部は、前記パラメータ推定のための反復計算において使用される計算式における時点毎の計算を、ベクトル又は行列又はテンソルを利用して全時点分についてまとめて実行する
   請求項１ないし３のうちいずれか１項に記載のパラメータ推定装置。
5. ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定システムであって、
   イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースを格納するデータ格納部と、
   前記データベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算部と、
   前記演算部により推定されたパラメータを出力する出力部と
   を備えるパラメータ推定システム。
6. ｃｏｘ比例ハザードモデルのパラメータ推定を、秘密計算により実行するパラメータ推定装置により実行されるパラメータ推定方法であって、
   イベントが観測された時点と、当該時点の観測対象の特徴量と、当該時点の観測対象の状態とを含むレコードを、観測対象毎に有するデータベースから、時点からなるベクトルを読み出し、当該ベクトルをソートすることにより、置換表と、時点の境目を示すフラグとを生成し、前記置換表と、前記フラグとを用いることにより、時点の値を秘匿したまま、前記特徴量の時点毎の集計を行い、集計結果に基づいて前記パラメータ推定を行う演算ステップと、
   前記演算ステップにより推定されたパラメータを出力する出力ステップと
   を備えるパラメータ推定方法。
7. コンピュータを、請求項１ないし４のうちいずれか１項に記載のパラメータ推定装置における各部として機能させるためのプログラム。

Images