(57)【特許請求の範囲】
【請求項１】
Ｎを撹乱対象となるデータのデータ数、Ｍを撹乱対象となるデータの属性の数、Ａｊをｊ番目の属性（ただし、ｊは１≦ｊ≦Ｍなる整数）とし、
撹乱対象となるデータの確率密度関数と撹乱データの確率密度関数の比である密度比モデルｗ（ｚ）を用いて、撹乱対象となるデータの統計値を推定する撹乱データ再構築装置であって、
属性Ａｊがカテゴリ属性である場合は維持確率ρｊと濃度｜Ａｊ｜を用いて算出される属性Ａｊの条件付き確率
【数２１】
0006505574.tif 000022
に従うランダムな処理を施す維持置換撹乱、属性Ａｊが数値属性である場合は有界ラプラス分布のパラメータφｊと属性Ａｊの値域[ａｊ，ｂｊ]（ａｊ，ｂｊはａｊ≦ｂｊを満たす実数）を用いて算出される属性Ａｊの条件付き確率
【数２２】
0006505574.tif 000023
に従うランダムな処理を施す有界ラプラスノイズ加算が施された撹乱データの標本｛ｙｉ｝ｉ＝１…Ｎと、前記維持確率ρｊ、前記濃度｜Ａｊ｜、前記パラメータφｊ、前記値域[ａｊ，ｂｊ]とから、Ｍ個の属性Ａｊの条件付き確率ＰＹ｜ＸＡｊ（ｖ’｜ｖ）（１≦ｊ≦Ｍ）の積である条件付き確率ＰＹ｜Ｘ（ｙ｜ｘ）を用いて、１≦ｉ≦Ｎ、１≦ｊ≦Ｎを満たす整数ｉ、ｊの組に対して条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）を算出する条件付き確率算出部と、
前記条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）から、
【数２３】
0006505574.tif 000024
を条件として
【数２４】
0006505574.tif 000025
を満たすｗを前記密度比モデルｗ（ｚ）として算出する密度比モデル算出部と
を含む撹乱データ再構築装置。
【請求項２】
請求項１に記載の撹乱データ再構築装置であって、
α（ｉ）（ただし、ｉは１≦ｉ≦Ｎなる整数）を実数とし、
前記密度比モデルｗ（ｚ）はカーネル関数の族｛ｋ（ｚｊ，ｚｉ）｝（ただし、ｉ、ｊは１≦ｉ≦Ｎ、１≦ｊ≦Ｎなる整数）を用いて
【数２５】
0006505574.tif 000026
と表現されることを特徴とする撹乱データ再構築装置。
【請求項３】
請求項２に記載の撹乱データ再構築装置であって、
前記密度比モデル算出部は、
前記族｛ｋ（ｚｊ，ｚｉ）｝（ただし、ｉ、ｊは１≦ｉ≦Ｎ、１≦ｊ≦Ｎなる整数）を用いて、前記標本｛ｙｉ｝ｉ＝１…Ｎから、ｋ（ｙｉ，ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列であるカーネル行列Ｋを算出する密度比行列算出部と、
前記ＰＹ｜Ｘ（ｙｉ｜ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列である条件付き確率行列Ｐと前記カーネル行列Ｋから、前記実数α（ｉ）を第ｉ成分とするＮ次元ベクトルである線形結合ベクトルαを算出する線形結合ベクトル算出部と
を含む撹乱データ再構築装置。
【請求項４】
請求項１に記載の撹乱データ再構築装置であって、
α（ｉ）（ただし、ｉは１≦ｉ≦Ｒなる整数）を実数とし、
前記密度比モデルｗ（ｚ）は基底関数の族｛ζｉ（ｚｊ）｝（ただし、ｉ、ｊは１≦ｉ≦Ｒ、１≦ｊ≦Ｎなる整数）を用いて
【数２６】
0006505574.tif 000027
と表現されることを特徴とする撹乱データ再構築装置。
【請求項５】
請求項４に記載の撹乱データ再構築装置であって、
前記密度比モデル算出部は、
前記族｛ζｉ（ｚｊ）｝（ただし、ｉ、ｊは１≦ｉ≦Ｒ、１≦ｊ≦Ｎなる整数）を用いて、前記標本｛ｙｉ｝ｉ＝１…Ｎから、ζｉ（ｙｊ）をｊ行ｉ列の要素とするＮ×Ｒ行列である基底関数行列Ｘを算出する密度比行列算出部と、
前記ＰＹ｜Ｘ（ｙｉ｜ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列である条件付き確率行列Ｐと前記基底関数行列Ｘから、前記実数α（ｉ）を第ｉ成分とするＲ次元ベクトルである線形結合ベクトルαを算出する線形結合ベクトル算出部と
を含む撹乱データ再構築装置。
【請求項６】
Ｎを撹乱対象となるデータのデータ数、Ｍを撹乱対象となるデータの属性の数、Ａｊをｊ番目の属性（ただし、ｊは１≦ｊ≦Ｍなる整数）とし、
撹乱データ再構築装置が、撹乱対象となるデータの確率密度関数と撹乱データの確率密度関数の比である密度比モデルｗ（ｚ）を用いて、撹乱対象となるデータの統計値を推定する撹乱データ再構築方法であって、
前記撹乱データ再構築装置が、属性Ａｊがカテゴリ属性である場合は維持確率ρｊと濃度｜Ａｊ｜を用いて算出される属性Ａｊの条件付き確率
【数２７】
0006505574.tif 000028
に従うランダムな処理を施す維持置換撹乱、属性Ａｊが数値属性である場合は有界ラプラス分布のパラメータφｊと属性Ａｊの値域[ａｊ，ｂｊ]（ａｊ，ｂｊはａｊ≦ｂｊを満たす実数）を用いて算出される属性Ａｊの条件付き確率
【数２８】
0006505574.tif 000029
に従うランダムな処理を施す有界ラプラスノイズ加算が施された撹乱データの標本｛ｙｉ｝ｉ＝１…Ｎと、前記維持確率ρｊ、前記濃度｜Ａｊ｜、前記パラメータφｊ、前記値域[ａｊ，ｂｊ]とから、Ｍ個の属性Ａｊの条件付き確率ＰＹ｜ＸＡｊ（ｖ’｜ｖ）（１≦ｊ≦Ｍ）の積である条件付き確率ＰＹ｜Ｘ（ｙ｜ｘ）を用いて、１≦ｉ≦Ｎ、１≦ｊ≦Ｎを満たす整数ｉ、ｊの組に対して条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）を算出する条件付き確率算出ステップと、
前記撹乱データ再構築装置が、前記条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）から、
【数２９】
0006505574.tif 000030
を条件として
【数３０】
0006505574.tif 000031
を満たすｗを前記密度比モデルｗ（ｚ）として算出する密度比モデル算出ステップと
を含む撹乱データ再構築方法。
【請求項７】
請求項１ないし５のいずれか１項に記載の撹乱データ再構築装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、データベースにおける個別データを確率的手法により秘匿化したデータから個別データの統計値を推定する技術に関し、特に個別データの統計値を推定する際に用いる確率密度関数を算出するものに関する。
【背景技術】
【０００２】
データベースにおける個別データを確率的手法により秘匿化することを撹乱という。また、秘匿化したデータ（以下、撹乱データという）から元データの統計値（以下、再構築データという）を推定することを再構築という。つまり、撹乱とは撹乱データからなる秘匿データベースを作成することであり、再構築とはこの秘匿データベースに対し統計分析を行い、再構築データを得ることである。
【０００３】
このような撹乱・再構築の技術として、維持置換撹乱（非特許文献１、非特許文献２、非特許文献４）、有界ラプラスノイズ加算（非特許文献３）がある。これらの技術における再構築処理では、撹乱データから再構築データを推定する際、様々な統計分析を可能にするために、元データの生成分布を表す確率密度関数の推定を行っている。確率密度関数を知ることはデータの生成規則を知ることと等価であり、あらゆる統計分析が可能となるからである。
【０００４】
これらの技術では、確率密度関数を多次元ヒストグラム（多次元度数表）として表現し、撹乱前の元データの生成分布の推定を行っている。撹乱、再構築の各処理のイメージを示すために、元データ、撹乱データ、再構築データを図１に示す。この例では再構築データの縦軸を身長、横軸を性別とする２次元ヒストグラムで表現しており、ヒストグラムの各値は百分率の頻度を表している。
【０００５】
非特許文献４では、以下に示す尤度関数と呼ばれる関数を最大化するＰＸ（ｘ）を求める問題を解くことにより、再構築処理を行っている。
【０００６】
【数１】
0006505574.tif 000002
【０００７】
ＰＸ（ｘ）は元データの確率密度関数、ＰＹ（ｙ）は撹乱データの確率密度関数を表す。また、ＰＹ｜Ｘ（ｙ｜ｘ）は撹乱前のデータベース中のデータｘが撹乱処理により撹乱データｙとなる確率を表す。なお、変数ｘ、ｙを省略して、ＰＸ（ｘ）、ＰＹ（ｙ）、ＰＹ｜Ｘ（ｙ｜ｘ）を単にＰＸ、ＰＹ、ＰＹ｜Ｘと表すこともある。
【０００８】
撹乱データとして与えられた標本から算出される撹乱データの度数分布をＰＹ（ｙ）とし、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ  Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムにより、ＰＸ（ｘ）を多次元ヒストグラムとして求めていた。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】五十嵐大、千田浩司、高橋克巳、“多値属性に適用可能な効率的プライバシー保護クロス集計”、コンピュータセキュリティシンポジウム２００８  論文集、2008年10月、第2008巻、pp.497-502
【非特許文献２】五十嵐大、千田浩司、高橋克巳、“k-匿名性の確率的指標への拡張とその適用例”、コンピュータセキュリティシンポジウム２００９ 論文集、2009年10月、第2009巻、pp.1-6
【非特許文献３】五十嵐大、長谷川聡、納竜也、菊池亮、千田浩司、“数値属性に適用可能な, ランダム化によりk-匿名性を保証するプライバシー保護クロス集計”、コンピュータセキュリティシンポジウム２０１２  論文集、2012年10月、pp.639-646
【非特許文献４】Rakesh Agrawal, Ramakrishnan Srikant, and Dilys Thomas, “Privacy Preserving OLAP”, In Proceedings of the 2005 ACM SIGMOD International Conference on Management of Data, 2005, pp.251-262
【発明の概要】
【発明が解決しようとする課題】
【００１０】
上述したように、非特許文献１〜４の技術ではいずれも再構築処理を行う際、元データの確率密度関数を表現する多次元ヒストグラムを算出する。このため、データベース中の属性の数が増えるにつれて指数関数的に空間計算量が増えてしまい、計算が困難になっていた。例えば、１０個の属性があり、それぞれの属性が取りうる値が１０であるとすると、１０１０個の要素のためのメモリ空間を確保しなければならず、たちまち計算が困難になってしまう。その結果、重回帰分析やロジスティック回帰分析といった属性を多く用いる多変量解析を行うことが困難であった。
【００１１】
そこで本発明では、元データの確率密度関数を多次元ヒストグラムとして算出する代わりに、撹乱データの標本から元データの確率密度関数と撹乱データの確率密度関数の比を算出する撹乱データ再構築装置を提供することを目的とする。
【課題を解決するための手段】
【００１２】
本発明の一態様は、Ｎを撹乱対象となるデータのデータ数、Ｍを撹乱対象となるデータの属性の数、Ａｊをｊ番目の属性（ただし、ｊは１≦ｊ≦Ｍなる整数）、｛ｙｉ｝ｉ＝１…Ｎを属性Ａｊがカテゴリ属性である場合は維持置換撹乱、属性Ａｊが数値属性である場合は有界ラプラスノイズ加算が施された撹乱データの標本、ρｊを属性Ａｊがカテゴリ属性である場合の維持確率、｜Ａｊ｜を属性Ａｊがカテゴリ属性である場合の濃度、φｊを属性Ａｊが数値属性である場合の有界ラプラス分布のパラメータ、[ａｊ，ｂｊ]を属性Ａｊが数値属性である場合の値域（ａｊ，ｂｊはａｊ≦ｂｊを満たす実数）、ｗ（ｚ）を撹乱対象となるデータの確率密度関数と撹乱データの確率密度関数の比である密度比モデルとし、前記標本｛ｙｉ｝ｉ＝１…Ｎから前記密度比モデルｗ（ｚ）を算出する撹乱データ再構築装置であって、属性Ａｊがカテゴリ属性である場合は前記維持確率ρｊと前記濃度｜Ａｊ｜を用いて算出される条件付き確率
【００１３】
【数２】
0006505574.tif 000003
【００１４】
、属性Ａｊが数値属性である場合は前記パラメータφｊと前記値域[ａｊ，ｂｊ]を用いて算出される条件付き確率
【００１５】
【数３】
0006505574.tif 000004
【００１６】
を属性Ａｊの条件付き確率ＰＹ｜ＸＡｊ（ｖ’｜ｖ）としてＭ個の条件付き確率ＰＹ｜ＸＡｊ（ｖ’｜ｖ）（１≦ｊ≦Ｍ）の積である条件付き確率ＰＹ｜Ｘ（ｙ｜ｘ）と前記標本｛ｙｉ｝ｉ＝１…Ｎを用いて、１≦ｉ≦Ｎ、１≦ｊ≦Ｎを満たす整数ｉ、ｊの組に対して条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）を算出する条件付き確率算出部と、
【００１７】
【数４】
0006505574.tif 000005
【００１８】
を条件として
【００１９】
【数５】
0006505574.tif 000006
【００２０】
を満たすｗ（ｚ）を密度比モデルとして算出する密度比モデル算出部とを含む。
【発明の効果】
【００２１】
本発明によれば、元データの確率密度関数の代わりに元データの確率密度関数と撹乱データの確率密度関数の比を算出することにより、再構築計算に必要となる空間計算量を削減することが可能となる。
【図面の簡単な説明】
【００２２】
【図１】撹乱・再構築の処理イメージを示すための、元データ・撹乱データ・再構築データの一例を示す図。
【図２】密度比モデル決定アルゴリズム（カーネル密度比モデル）を示す図。
【図３】密度比モデル決定アルゴリズム（線形密度比モデル）を示す図。
【図４】実施例１の撹乱データ再構築装置１００の構成を示すブロック図。
【図５】実施例１の撹乱データ再構築装置１００の動作を示すフローチャート。
【図６】実施例１の密度比モデル算出部１２０の構成を示すブロック図。
【図７】実施例１の密度比モデル算出部１２０の動作を示すフローチャート。
【発明を実施するための形態】
【００２３】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００２４】
撹乱対象となるデータ、つまり元データのデータ数をＮ、元データの属性の数をＭとする。また、属性はカテゴリ属性と数値属性が混在して含まれているものとする。カテゴリ属性の例としては、図１の性別があり、その属性値の集合は（男、女）となる。数値属性の例としては、図１の身長があり、その属性値の集合は（０ｃｍ〜２００ｃｍ）となる。ｊ番目の属性（つまり、属性値の集合）をＡｊと表し（１≦ｊ≦Ｍ）、すべての属性の集合をＡ＝Ａ１×…×ＡＭとする。
【００２５】
先ほどと同様、ＰＸ（ｘ）は元データの確率密度関数、ＰＹ（ｙ）は撹乱データの確率密度関数を表すものとする。ＰＸ（ｘ）、ＰＹ（ｙ）はいずれも定義域をＡ＝Ａ１×…×ＡＭ、値域を[０，１]とする関数である。なお、変数ｘ、ｙはＭ次元列ベクトルとして扱うこととする。
【００２６】
元データは確率密度ＰＸ（ｘ）を持つ確率分布から得られるＮ個の標本データ｛ｘｉ｝ｉ＝１…Ｎで与えられるとする（ただし、ｘｉはＡの要素）。また、撹乱データは、条件付き確率ＰＹ｜Ｘ（ｙ｜ｘ）に従って元データがランダムに撹乱されているものとし、撹乱後の撹乱データを｛ｙｉ｝ｉ＝１…Ｎとする。また、ｘｉ（ｊ）は、ｉ番目のデータのｊ番目の属性の値を表す。元データｘｉ、撹乱データｙｉはいずれもＭ次元列ベクトルとして表している。
【００２７】
図１を例にすると、ｘ１＝（１７０．５  男）Ｔとなり、ｘ１（１）＝１７０．５、ｘ１（２）＝男である（Ｔは転置を表す）。
【実施例１】
【００２８】
実施例１の発明は、撹乱データから元データの確率密度関数と撹乱データの確率密度関数の比を得るものである。撹乱データの生成には、カテゴリ属性に対しては非特許文献１、２の維持置換撹乱が、数値属性に対しては非特許文献３の有界ラプラスノイズ加算が用いられているものとする。そこで、まずカテゴリ属性の撹乱方法、数値属性の撹乱方法について説明する。
【００２９】
（カテゴリ属性の撹乱方法）
カテゴリ属性に対しては、維持確率ρで属性値を維持し、１−ρの確率で属性値をランダムに変更することで、データを秘匿化する処理が行われる。すなわち、あるカテゴリ属性Ａｊの属性値ｖが属性値ｖ’に変わる条件付き確率ＰＹ｜ＸＡｊ（ｖ’｜ｖ）は、属性Ａｊの維持確率ρｊを用いて
【００３０】
【数６】
0006505574.tif 000007
【００３１】
と表される（ｊは１≦ｊ≦Ｍを満たす整数）。なお、｜Ａｊ｜は属性Ａｊの取り得る値の数（つまり、集合Ａｊの濃度）を示す。
【００３２】
カテゴリ属性に対する撹乱では、この条件付き確率に従うランダムな処理が施される。また、属性Ａｊの維持確率ρｊは公開されるものとする。
【００３３】
（数値属性の撹乱方法）
数値属性に対しては、属性値に有界ラプラス分布に従うノイズを付与する（つまり、有界ラプラスノイズ加算する）ことで撹乱を行う。有界ラプラス分布とは、台の上限と下限が定まっているラプラス分布のことであり、この有界ラプラス分布に従うノイズを付与することで、元データを秘匿化する。ある数値属性Ａｊの属性値ｖが属性値ｖ’に変わる条件付き確率密度ＰＹ｜ＸＡｊ（ｖ’｜ｖ）は、属性Ａｊの有界ラプラス分布のパラメータφｊを用いて
【００３４】
【数７】
0006505574.tif 000008
【００３５】
となる（ｊは１≦ｊ≦Ｍを満たす整数）。ここでγｊ（ｖ）はラプラス分布を有界にしたことによって生じた、有界ラプラス分布を調整するための関数、[ａｊ，ｂｊ]は属性Ａｊの値域（ただし、ａｊ，ｂｊはａｊ≦ｂｊを満たす実数）である。
【００３６】
【数８】
0006505574.tif 000009
【００３７】
数値属性に対する撹乱では、この条件付き確率密度に従うランダムな処理が施される。また、属性Ａｊの有界ラプラス分布のパラメータφｊは公開されるものとする。
【００３８】
以上の撹乱処理が施された撹乱データ｛ｙｉ｝ｉ＝１…Ｎ、カテゴリ属性Ａｊの維持確率ρｊ、カテゴリ属性Ａｊの濃度｜Ａｊ｜、数値属性Ａｊの有界ラプラス分布のパラメータφｊ、数値属性Ａｊの値域[ａｊ，ｂｊ]を用いて、元データの確率密度関数と撹乱データの確率密度関数の比ＰＸ（ｚ）／ＰＹ（ｚ）を推定する。以降、この確率密度関数の比をｗ（ｚ）＝ＰＸ（ｚ）／ＰＹ（ｚ）と表し、ｗ（ｚ）を密度比モデルということとする。したがって、ｗ（ｚ）は定義域をＡ＝Ａ１×…×ＡＭとする実数値関数となる。なお、変数ｚをＭ次元列ベクトルとして扱うこととする。また、ｗ（ｚ）を単にｗと表すこともある。
【００３９】
非特許文献４の技術では、撹乱データから元データの統計値を得る際に、式（１）を解くことにより元データの確率密度関数（多次元ヒストグラム）を推定していた。実施例１の発明では、元データの確率密度関数の代わりに密度比モデルｗ（ｚ）を推定する。その際、式（１）と同様の以下の式（４）を解く。
【００４０】
【数９】
0006505574.tif 000010
【００４１】
ここで、密度比モデルｗ（ｚ）の定義および撹乱データの標本｛ｙｉ｝ｉ＝１…Ｎから算出される度数分布でＰＹ（ｙ）を近似することより、
【００４２】
【数１０】
0006505574.tif 000011
【００４３】
となる。
式（４）および式（５）より、
【００４４】
【数１１】
0006505574.tif 000012
【００４５】
なお、定数倍はａｒｇｍａｘに影響しないため、予め式から除いている。
ＰＸ（ｘ）が確率密度関数であることから、この最適化問題を解く際に以下の制約を設ける。
【００４６】
【数１２】
0006505574.tif 000013
【００４７】
これらを用いて、密度比モデルｗ（ｚ）を推定する問題を定式化すると、以下のようになる。つまり、式（７）を条件として式（６）の尤度関数を最大化する密度比モデルｗ（ｚ）を算出する問題として定式化される。
【００４８】
【数１３】
0006505574.tif 000014
【００４９】
ここで、ＰＹ｜Ｘ（ｙ｜ｘ）は属性ごとの条件付き確率の積であり、属性がカテゴリ属性である場合は式（２）を、数値属性である場合は式（３）を用いて以下の式（８）を計算する。
【００５０】
【数１４】
0006505574.tif 000015
【００５１】
以下、密度比モデルｗ（ｚ）はカーネル密度比モデルまたは線形密度比モデルのいずれかであるとして、式（６）の最適化問題を解くアルゴリズムについて説明する。
【００５２】
まず、密度比モデルｗがカーネル密度比モデルである場合における式（６）の最適化問題を解くアルゴリズムについて説明する。
【００５３】
（密度比モデルｗがカーネル密度比モデルである場合）
密度比モデルｗがカーネル密度比モデルであるとは、密度比モデルｗが以下の式（９）にて表現されることをいう。
【００５４】
【数１５】
0006505574.tif 000016
【００５５】
ここで、ｋ（ｚｊ，ｚｉ）はカーネル関数と呼ばれるものである。また、α（ｉ）は実数である。つまり、カーネル密度比モデルはカーネル関数の線形結合として表現される密度比モデルである。
【００５６】
なお、参考非特許文献５の５章にあるように、カーネル関数には様々なものがある。
（参考非特許文献５： Trevor Hastie, Robert Tibshirani and Jerome Friedman, “The Elements of Statistical Learning: Data Mining, Inference, and Prediction.”, Second Edition, Springer-Verlag, 2009）
例えば、
【００５７】
【数１６】
0006505574.tif 000017
【００５８】
などがある。これらのカーネル関数の中からユーザ（密度比モデルｗを決定する者）がどの関数を用いるかを自ら選択する。また、必要に応じてカーネル関数で用いるパラメータ（例えば、上記式（１１）中のσ）もユーザが設定する。
【００５９】
｛ｋ（ｚｊ，ｚｉ）｝（ただし、ｉ、ｊは１≦ｉ≦Ｎ、１≦ｊ≦Ｎなる整数）をＮ２個のカーネル関数を持つ関数族（以下、カーネル関数の族ともいう）とし、式（６）、式（７）に式（９）を適用したものが以下の式（１２）、式（１３）となる。
【００６０】
【数１７】
0006505574.tif 000018
【００６１】
ただし、α＝（α（１），…，α（Ｎ））Ｔである。
【００６２】
つまり、式（６）の最適化問題を解いて密度比モデルｗを求めることが、式（１３）を条件として式（１２）の最適化問題を解いてパラメータα＝（α（１），…，α（Ｎ））Ｔを求めることに還元される。式（１２）は、参考非特許文献６に記載されている方法（勾配法によるパラメータ更新および制約条件の充足）と同様の方法により、解くことが可能である。
（参考非特許文献６： Masashi Sugiyama, Shinichi Nakajima, Hisashi Kashima, Paul von Buenau and Motoaki Kawanabe, “Direct Importance Estimation with Model Selection and Its Application to Covariate Shift Adaptation”, In Advances in Neural Information Processing Systems, 2008, pp.1433-1440）
【００６３】
カーネル密度比モデルを特徴づけるパラメータα＝（α（１），…，α（Ｎ））Ｔを求めるアルゴリズムを図２に示す（以下、αを線形結合ベクトルという）。ここで、ＰはＰＹ｜Ｘ（ｙｉ｜ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列（以下、Ｐを条件付き確率行列という）、Ｋはｋ（ｙｉ，ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列（以下、Ｋをカーネル行列という）である。また、βは学習率と呼ばれる値（０＜β＜１）、εは収束の判定のための閾値（ε＞０）を表す。なお、βが１に近づくほどアルゴリズムの収束が早くなり、０に近づくほどアルゴリズムの収束が遅くなる。高くしすぎると収束しなくなることがあるため、適切な値をユーザが決める必要がある。
【００６４】
ｃｏｌＳｕｍ（・）は列単位の総和を求める関数、ｍａｘ（０，・）は列ベクトルの要素ごとに０と比較をし、その要素が０よりも小さい場合はその要素を０で置き換えたベクトルを返す関数、．／は要素ごとの除算（element-wise division）、１Ｎはすべての要素が１であるＮ次元列ベクトルである。Ａ、Ｂをそれぞれａｐ，ｑ、ｂｐ，ｑをｐ行ｑ列の要素とするＰ×Ｑ行列（ｐ、ｑは１≦ｐ≦Ｐ、１≦ｑ≦Ｑなる整数）とすると、Ａ．／Ｂは、ａｐ，ｑ／ｂｐ，ｑをｐ行ｑ列の要素とするＰ×Ｑ行列となる。要素ごとの除算．／は、結合力が弱く最後に計算されるため、図２の６行目の１Ｎ．／ＰＫαｔは１Ｎの各要素をＰＫαｔの各要素で除したＮ次元列ベクトルを示す。
【００６５】
なお、アルゴリズム中に出てくるｂはＮ次元ベクトルである。
【００６６】
次に、密度比モデルｗが線形密度比モデルである場合における式（６）の最適化問題を解くアルゴリズムについて説明する。
【００６７】
（密度比モデルｗが線形密度比モデルである場合）
密度比モデルｗが線形密度比モデルであるとは、密度比モデルｗが以下の式（１４）にて表現されることをいう。
【００６８】
【数１８】
0006505574.tif 000019
【００６９】
ここで、ζｉ（ｚｊ）は基底関数と呼ばれるものである。また、α（ｉ）は実数である。つまり、線形密度比モデルは基底関数の線形結合として表現される密度比モデルである。
【００７０】
なお、参考非特許文献５にあるように、基底関数にも様々なものがある。例えば、
【００７１】
【数１９】
0006505574.tif 000020
【００７２】
などがある。なお、ｚｊ（ｉ）は先述の通りベクトルｚｊの第ｉ成分を表す。
【００７３】
これらの基底関数の中からユーザがどの関数を用いるかを自ら選択する。また、必要に応じて基底関数で用いるパラメータもユーザが設定する。
【００７４】
｛ζｉ（ｚｊ）｝（ただし、ｉ、ｊは１≦ｉ≦Ｒ、１≦ｊ≦Ｎなる整数）をＲ×Ｎ個の基底関数を持つ関数族（以下、基底関数の族ともいう）とし、式（６）、式（７）に式（１４）を適用したものが以下の式（１７）、式（１８）となる。
【００７５】
【数２０】
0006505574.tif 000021
【００７６】
ただし、α＝（α（１），…，α（Ｒ））Ｔである。
【００７７】
つまり、式（６）の最適化問題を解いて密度比モデルｗを求めることが、式（１８）を条件として式（１７）の最適化問題を解いてパラメータα＝（α（１），…，α（Ｒ））Ｔを求めることに還元される。式（１７）も、式（１２）を解くためのアルゴリズムと同様のアルゴリズムで解くことが可能である。
【００７８】
線形密度比モデルを特徴づけるパラメータα＝（α（１），…，α（Ｒ））Ｔを求めるアルゴリズムを図３に示す（以下、αを線形結合ベクトルという）。Ｘはζｉ（ｙｊ）をｊ行ｉ列の要素とするＮ×Ｒ行列である（以下、Ｘを基底関数行列という）。図３のアルゴリズムは図２のアルゴリズム中のＫをＸで書き換えたものとなっている。Ｐ、β、εについては図２のそれと同様である。
【００７９】
なお、図３のアルゴリズム中に出てくるｂは、図２のｂとは異なり、Ｒ次元ベクトルとなる。
【００８０】
密度比モデルｗとしてカーネル密度比モデルまたは線形密度比モデルを仮定することにより、勾配法を用いて元データの確率密度関数と撹乱データの確率密度関数の比（密度比モデルｗ）を特徴付けるパラメータである線形結合ベクトルαを算出することができる。この際、算出に必要な空間計算量は、カーネル密度比モデルであればＯ（Ｎ２）、線形密度比モデルであればＯ（ＮＲ）である。
【００８１】
また、この線形結合ベクトルαを用いて表現される密度比モデルｗを用いることで、多変量解析を行うことが可能となる。例えば、参考非特許文献６に記載されている重点サンプリングを用いることで、多変量解析だけでなく、機械学習といった分析も可能となる。
【００８２】
以下、図４〜図５を参照して実施例１の撹乱データ再構築装置１００について説明する。図４は、実施例１の撹乱データ再構築装置１００の構成を示すブロック図である。図５は、実施例１の撹乱データ再構築装置１００の動作を示すフローチャートである。
【００８３】
図４に示すように撹乱データ再構築装置１００は、条件付き確率算出部１１０と、密度比モデル算出部１２０を含む。
【００８４】
条件付き確率算出部１１０は、撹乱データの標本｛ｙｉ｝ｉ＝１…Ｎとこの撹乱データの生成を特徴付けるパラメータであるカテゴリ属性Ａｊの維持確率ρｊ及び濃度｜Ａｊ｜、数値属性Ａｊの有界ラプラス分布のパラメータφｊ及び値域[ａｊ，ｂｊ]から条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）を算出する（Ｓ１１０）。ｉ、ｊを１≦ｉ≦Ｎ、１≦ｊ≦Ｎを満たす整数として各（ｉ，ｊ）についてＰＹ｜Ｘ（ｙｉ｜ｙｊ）は式（８）、式（２）、式（３）を用いて算出することができる。
【００８５】
密度比モデル算出部１２０は、条件付き確率ＰＹ｜Ｘ（ｙｉ｜ｙｊ）（１≦ｉ≦Ｎ、１≦ｊ≦Ｎ）を用いて、密度比モデルｗを算出する（Ｓ１２０）。具体的には、式（７）を条件として式（６）の尤度関数を最大化する密度比モデルｗを求めればよい。
【００８６】
以下、図６〜図７を参照して実施例１の密度比モデル算出部１２０について説明する。図６は、実施例１の密度比モデル算出部１２０の構成を示すブロック図である。図７は、実施例１の密度比モデル算出部１２０の動作を示すフローチャートである。
【００８７】
図６に示すように密度比モデル算出部１２０は、密度比行列算出部１２１と、線形結合ベクトル初期化部１２２と、線形結合ベクトル算出部１２３と、出力結果生成部１２４とを含む。なお、密度比モデル算出部１２０は図６には図示しない記録部を有しており、記録部には学習率βと閾値εが記録されているものとする。また、以下ではカーネル行列Ｋ、基底関数行列Ｘのことを密度比行列Ｄということとする。
【００８８】
以下、密度比モデルがカーネル密度比モデルまたは線形密度比モデルであるものとして説明する。また、条件付き確率算出部１１０の算出結果であるＰＹ｜Ｘ（ｙｉ｜ｙｊ）（１≦ｉ≦Ｎ、１≦ｊ≦Ｎ）は、ＰＹ｜Ｘ（ｙｉ｜ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列である条件付き確率行列Ｐとして表現されているものとする。
【００８９】
密度比行列算出部１２１は、撹乱データの標本｛ｙｉ｝ｉ＝１…Ｎから密度比行列Ｄを算出する（Ｓ１２１）。密度比モデルがカーネル密度比モデルの場合、カーネル行列ＫをＤとして算出する。Ｋはｋ（ｙｉ，ｙｊ）をｉ行ｊ列の要素とするＮ×Ｎ行列であるので、例えば、式（１０）や式（１１）の定義にしたがって算出すればよい。また、密度比モデルが線形密度比モデルの場合、基底関数行列ＸをＤとして算出する。Ｘはζｉ（ｙｊ）をｊ行ｉ列の要素とするＮ×Ｒ行列であるので、例えば、式（１５）や式（１６）の定義にしたがって算出すればよい。なお、カーネル行列Ｋを特徴付けるカーネル関数の族｛ｋ（ｚｊ，ｚｉ）｝、基底関数行列Ｘを特徴付ける基底関数の族｛ζｉ（ｚｊ）｝、これらの関数の族で用いるパラメータについては、先述の通り、ユーザが事前に自ら選択しているものとする。
【００９０】
線形結合ベクトル初期化部１２２は、線形結合ベクトルを初期化する（Ｓ１２２）。初期化した線形結合ベクトルを初期線形結合ベクトルといい、α０で表す。密度比モデルがカーネル密度比モデルの場合、α０はＮ次元ベクトルである。また、密度比モデルが線形密度比モデルの場合、α０はＲ次元ベクトルである。いずれの場合もこの列ベクトルの各要素に０以上の値を設定することで初期化する。なお、Ｓ１２１とＳ１２２の処理の順番は入れ替わってもよい。
【００９１】
線形結合ベクトル算出部１２３は、条件付き確率行列Ｐ、密度比行列Ｄ、学習率βを用いて、線形結合ベクトルαを算出する（Ｓ１２３）。条件付き確率行列Ｐ、密度比行列Ｄ、学習率βを用いて初期線形結合ベクトルα０から順次線形結合ベクトルα１、α２、…を算出する。具体的な算出方法は、密度比モデルがカーネル密度比モデルの場合は図２、線形密度比モデルの場合は図３にある通りであり、いずれの場合もαｔ＋１とαｔの差分ベクトルのノルムの自乗（||αｔ＋１−αｔ||２）が閾値εより小さくなったとき計算を停止し、このときのαｔ＋１を線形結合ベクトルαとして算出する。
【００９２】
出力結果生成部１２４は、線形結合ベクトルαと、カーネル関数の族｛ｋ（ｚｊ，ｚｉ）｝あるいは基底関数の族｛ζｉ（ｚｊ）｝を用いて密度比モデルｗを生成し、出力する（Ｓ１２４）。
【００９３】
なお、密度比モデルｗがカーネル密度比モデルまたは線形密度比モデルである場合は線形結合ベクトルαで特徴付けられることから、ｗとαを同一視することによりｗの代わりにαを出力するのでもよい。この場合、図６においては出力結果生成部１２４が、図７においてはＳ１２４が不要となる。また、図４の密度比モデル算出部１２０の出力が密度比モデルｗでなく、線形結合ベクトルαとなる。
【００９４】
実施例１の発明では、あえてあらゆる分析に適用可能な元データの確率密度関数の算出を行わず、幅広い多変量解析に適用可能な密度比モデルｗ（密度比モデルを特徴付けるパラメータである線形結合ベクトルα）を算出する。
【００９５】
従来の方法（元データの確率密度関数を直接求める方法）では、空間計算量がＯ（｜Ａ１｜×…×｜ＡＭ｜）であった。このため、属性の数が増える（つまり、Ｍが大きくなる）につれ指数関数的に空間計算量が増え、メモリ容量が足りなくなり計算が困難となっていた。実施例１の発明では、属性の数が増えるにつれ指数関数的に空間計算量が増えることはない。つまり、計算に必要となる空間計算量を抑えることが可能となる。具体的には、密度比モデルがカーネル密度比モデルである場合はＯ（Ｎ２）、線形密度比モデルである場合はＯ（ＮＲ）の空間計算量に抑えることが可能となる。したがって、実施例１の発明を用いれば属性の数を増加した場合においても多変量解析が可能となる。
【００９６】
＜補記＞
本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ−ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【００９７】
ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。
【００９８】
ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。
【００９９】
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【０１００】
既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【０１０１】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１０２】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１０３】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１０４】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
0006505574.tif 000032
【図２】
0006505574.tif 000033
【図３】
0006505574.tif 000034
【図４】
0006505574.tif 000035
【図５】
0006505574.tif 000036
【図６】
0006505574.tif 000037
【図７】
0006505574.tif 000038