統計数理（2009）
第 57 巻 第 2 号 425–442
c 2009 統計数理研究所
［原著論文］

多重母集団寸法指標のノンパラメトリック最尤推定
2 時点の個票データへの適用

佐井 至道

†

（受付 2008 年 9 月 26 日；改訂 2009 年 6 月 1 日；採択 6 月 16 日）

要旨

標本調査で得られた個票データの公開におけるリスク評価では，標本寸法指標を基にして推定された母集団寸法指標を用いることが多い．その推定のために，初期の頃にはポアソンガン
マモデルが，最近ではピットマンモデルなどの超母集団モデルが用いられている．加えて，制約付きノンパラメトリック最尤推定法も提案され，計算時間上の問題も克服されつつある．
官庁統計の標本調査は継続調査として行われる場合が多いが，そのような場合にも，個々の時点の個票データに対して独立にリスク評価を行うのがこれまでは一般的であった．本論文では，継続調査で得られた複数の時点の個票データに対する同時リスク評価を行うために，寸法指標の概念を多重寸法指標に拡張し，多重標本寸法指標からの多重母集団寸法指標のノンパラ
メトリック推定法を提案する．
アメリカにおいて 1990 年と 2000 年に実施されたセンサスの 1％抽出個票データを題材に，多重母集団寸法指標の推定値に対する制約条件と，その条件をペナルティー関数として対数尤
度関数へ取り込む方法を検討する．また，推定を安定させるために，個々の時点の個票データ
からそれぞれ母集団寸法指標を推定した後に，その推定値を制約条件に組み込み，多重母集団
寸法指標を推定する方法も提案する．
キーワード： 多重指標，寸法指標，個票データ，リスク評価，ノンパラメトリック推定，官庁統計．
1.

はじめに個票データが標本調査で得られた場合，その公開のリスク評価のために標本寸法指標を基にして推定された母集団寸法指標を用いることが多い．寸法指標とは，いくつかの変数について
値の組み合わせが同じである個体数を求めた場合の，個体数別の頻度分布である． Bethlehem
et al.（1990）
がポアソンガンマモデルを提案して以来，様々な超母集団モデルと呼ばれる母集団寸法指標推定のためのモデルが提案されてきており，その中で，Pitman（1995）が提案したピットマンモデル（Ewens-Pitman sampling formula（EPSF）とも呼ばれる）の実データへの当て
はまりの良いことが報告されている．一方，佐井（2002, 2003）は，超母集団モデルを用いない
制約付きノンパラメトリック最尤推定法を提案し，佐井（2006）では，母集団寸法指標の推定値†

岡山商科大学 経済学部：〒700–8601 岡山市北区津島京町 2–10–1

426

統計数理

第 57 巻

第 2 号 2009

に対する制約条件をペナルティー関数で置き換えるなど，計算時間上の問題を克服するいくつかの方法について提案を行った．
官庁統計の標本調査では，例えば労働力調査が毎月実施されるように，継続調査として行われる場合がほとんどである．継続調査では，パネル調査のように，いったんサンプリングされた個体が数回に渡って調査対象となる場合と，原則的に毎回異なる個体がサンプリングされる場合がある．
いったんサンプリングされた個体が数回に渡って調査対象となる場合，同一個体のレコードをリンクして 1 つのレコードに加工した個票データを作成すると，その公開のリスク評価としては，その個票データに対応する仮想的な母集団の寸法指標の推定値を利用することが可能である．
一方で，調査時点ごとに異なる個体がサンプリングされる場合，これまでは個々の時点の個
票データからそれぞれ母集団寸法指標を推定することにより，リスク評価を行うのが一般的であった．
リスク評価によって個票データが危険と判断された場合，例えば年齢について，90 歳以上を実年齢でなく 90 歳以上とするトップコーディングや，5 歳階級ごとにまとめる丸めのように，
カテゴリーを併合する秘匿措置を用いて個体を識別しにくくする．そのような秘匿措置を段階的に施しながら母集団寸法指標を推定し，推定値そのもの，あるいはそれから計算される指標が基準を満足するところで，個票データを安全と判断する．指標としては，変数の値の組み合
わせが，母集団でも標本でも他のすべての個体と異なる（一意と呼ばれる）
個体数の推定値などが用いられる．
ところが，複数の時点の個票データに対して同じ秘匿措置を施した場合，ある時点の公開リ
スクが大きく減少しても，他の時点の公開リスクが同じように減少するとは限らない．その構造を把握して効率的な秘匿を行うためには，複数の時点における変数の同じ値の組み合わせに対する個体数の組の情報が必要になる．
本論文では，継続調査で得られた複数の時点の個票データに対する同時リスク評価を行うた
めに，寸法指標の概念を多重寸法指標に拡張し，多重標本寸法指標からの多重母集団寸法指標のノンパラメトリック最尤推定法を提案する．その導入により，どのような秘匿措置を施せば，複数の時点の個票データのリスクを効率的に減少できるのか，あるいは個々の時点の個票デー
タに対して異なる秘匿措置を施した方がよいのか，という秘匿措置の指針が得られるものと考えられる．
多重寸法指標とは，いくつかの変数の値の組み合わせが同じである個体数を複数の母集団について求めた場合の，個体数の組別の頻度分布である．定着しているとは言えないが，本論文では，母集団においては多重母集団寸法指標，標本においては多重標本寸法指標という用語を用いることにする．
母集団寸法指標を推定するための超母集団モデルを多重母集団寸法指標への推定に拡張するのは容易ではなく，これまで妥当なモデルは提案されていない．これに対して，ノンパラメト
リック最尤推定の対象を母集団寸法指標から多重母集団寸法指標に拡張するのは容易である．
まず 2.1 節では寸法指標の概念と，標本寸法指標に基づく母集団寸法指標のノンパラメトリッ
ク推定法について簡単に紹介する． 2.2 節では，寸法指標の概念を多重寸法指標に拡張するとともに，多重標本寸法指標を基にした多重母集団寸法指標のノンパラメトリック推定法を提案
する．
アメリカにおいて 1990 年と 2000 年に実施されたセンサスの 1％抽出個票データを 3.1 節で紹介し，3.2 節ではその個票データを題材として，多重母集団寸法指標の推定値に対する制約
条件と，その条件をペナルティー関数として対数尤度関数へ取り込む方法について提案する．

多重母集団寸法指標のノンパラメトリック最尤推定

427

3.3 節では，種々の条件の下での推定結果を示すが，その改善のために 3.4 節では 2 段階で推定を行う方法を提案する．
2.

多重母集団寸法指標のノンパラメトリック推定

2.1 標本寸法指標と母集団寸法指標
この節では，寸法指標の定義と，標本調査で得られた個票データを公開する際のリスク評価について，その概略を述べる．
母集団の N 個の個体が，個体を特定することが可能な数種類のキー変数の値の組み合わせに基づいて K 個のセルのいずれかに分けられているものとして，第 j 番目のセルのサイズ（含まれる個体数）を Fj (j = 1, 2, . . . , K) とする．ここで，サイズ l のセル数，すなわち Fj = l となるセル数を Sl (l = 0, 1, 2, . . . , L) とし，母集団寸法指標と呼ぶ． L はセルのサイズの最大値である．なお，本論文のノンパラメトリック推定では総セル数 K を用いない．
次に標本（個票データ）の大きさを n として，抽出率を λ = n/N と表す．標本では Fj ，Sl の代わりに fj ，sl (l = 0, 1, . . . , L) という表記を用いるが，後者が標本寸法指標である．
標本調査で個票データが得られている場合，計算された標本寸法指標を基にして母集団寸法指標を推定するのが，個票データのリスク評価としてしばしば用いられる方法である．母集団寸法指標の推定値 Ŝl が求まれば，例えば母集団でも標本でもサイズ 1 のセル数の推定値としてλŜ1 を用いることができる．ある個体のキー変数の値の組み合わせが母集団でも標本でも他のすべての個体と異なれば，その個体は識別可能，すなわち母集団のどの個体であるかを特定することが可能である．例えば個票データに他人に知られたくない情報が含まれていれば，特定された個体についてその情報も分かってしまうため最も危険と考えられる．また母集団でサイズ 2 以上であっても，標本としてそのうちの 1 個以上の個体が取られると，識別されないまでも予測される可能性がある．このように小さいサイズの母集団寸法指標の推定が重要である．
推定にはピットマンモデルなどの超母集団モデルが用いられることが多いが，佐井（2002,
では，標本寸法指標に基づく母集団寸法指標の制約付きノンパラメトリック推定法を提2003）案した．ここではその要点のみを簡潔に述べる．
非復元単純無作為抽出によって得られた標本から標本寸法指標 (s1 , s2 , . . . , sL ) が計算された
ときに，尤度を最大にするような非負の母集団寸法指標 (S1 , S2 , . . . , SL ) を求めるが，その尤度関数は抽出率が十分小さいときに(2.1)

Lap (S1 , S2 , . . . , SL | s1 , s2 , . . . , sL ) =L1
e−µl µl sl
·
n
N−n
sl !
N Cn λ (1 − λ)l =1と近似できる．ただし，
(2.2)

µl =

L






Sl · l Cl λl (1 − λ)l−ll=lであり，尤度関数の添字の ap は近似を意味する．
（2.1）
式を最大にする母集団寸法指標を推定値とするのが自然ではあるが，推定値は非常に不安定で，このままでは実用性に乏しい．そこで，佐井（2003）では，推定の際に母集団寸法指標の推定値にいくつかの制約を置く提案を行い，実データを基に比較を行った．そのうち，下記の 4 つの制約を課す場合の結果が最も良かった．
（a） Sl ≥ 0 （l = 1, 2, . . . , L）
（b） L
l=1 l · Sl = N（c） Sl ≥ Sl+1 （l = 1, 2, . . . , L − 1）

428

統計数理

第 57 巻

第 2 号 2009（d） 2 · log Sl ≤ log Sl−1 + log Sl+1 （l = 2, 3, . . . , L − 1）は母集団寸法指標の推定値の対数が下に凸であることを意味する．これまで扱っ
制約条件（d）
た多くの寸法指標は，概ねこれらの制約条件を満たしていた．
また，標本寸法指標を生成する可能性のあるすべての母集団寸法指標 (S1 , S2 , . . . , SL ) について網羅的に尤度を計算することは計算時間上困難であるため，最大尤度の数種類の探索法について比較を行った．
更に，佐井（2006）
では母集団寸法指標の推定値に対する制約条件をペナルティー関数で置き換えるなど，計算時間上の問題を克服する方法について提案を行った．ここではその詳細については省略し，ペナルティー関数については 3.2 節で紹介する．
2.2 多重標本寸法指標と多重母集団寸法指標
この節では，前節で述べた寸法指標の概念を多重寸法指標に拡張する．多重寸法指標については 渋谷・佐井（2007）, Sibuya and Sai（2008）が紹介して，多変量に拡張したピットマンモデ
ルの適用を検討しているが，実データへの当てはまりは良くない．また，佐井（2007, 2008）では 馬場・坂口（2006）の紹介した学会所属データに対してノンパラメトリック推定を行っている．このデータは多重寸法指標においてサイズの最大値が 1 となる特殊なケースである．
M 個の時点において同じ調査項目についての調査が行われたとして，第 i 時点における母集団を Di ，個票データが得られる標本を di と表す．また，第 i 時点における母集団の個体数を
Ni ，標本の個体数を ni ，抽出率を λi = ni /Ni とする．
ここで，各時点の母集団において前節と同様にセルを構成し，キー変数の同じ値の組み合わせのセルについて，各時点において含まれる個体数を考える．複数時点の母集団 (D1 , D2 , . . . , DM )においてサイズの組 (l1 , l2 , . . . , lM ) となるセル数を S(l1 ,l2 ,...,lM ) と表し，多重母集団寸法指標と呼ぶ．標本でも同様に，(d1 , d2 , . . . , dM ) においてサイズの組 (l1 , l2 , . . . , lM ) となるセル数をs(l1 ,l2 ,...,lM ) と表し，多重標本寸法指標と呼ぶ． Li を第 i 時点でのサイズの最大値とするとき
li = 0, 1, . . . , Li として，ここではサイズとして 0 を含めて議論する．なお，キー変数のとる値の数の積として総セル数 K は得られるが，その中にはあり得ない組み合わせなどが含まれることが多く，実質的な総セル数が確定できないことが多い．本論文で提案する方法でも K を確定しないことにする．そのため，もし K が確定できれば，少なくとも 1 時点でサイズが正となるセル数を K から引くことによって求められる s(0,0,...,0) も，ここでは求めない．また
S(0,0,...,0) は推定対象にしないことにする．
(0, 0, . . . , 0) を除く各時点におけるすべてのサイズの組み合わせに対する多重母集団寸法指標の組と多重標本寸法指標の組を，それぞれ(2.3)

S = (S(0,0,...,0,1) , . . . , S(0,0,...,0,LM ) , . . . ,S(L1 ,L2 ,...,LM −1 ,0) , . . . , S(L1 ,L2 ,...,LM −1 ,LM ) ),(2.4)

s = (s(0,0,...,0,1) , . . . , s(0,0,...,0,LM ) , . . . ,s(L1 ,L2 ,...,LM −1 ,0) , . . . , s(L1 ,L2 ,...,LM −1 ,LM ) )と表す．
各時点において独立に非復元単純無作為抽出された標本から s が得られたときの S の尤度
関数は



S(l1 ,l2 ,...,lM ) !
1

(2.5)
L(S | s) = M ,l ,...,l ) k(l1 ,l2 ,...,lM ),(l ,l ,...,l ) !
C
n
N
(l
i C
i=1 i
1 2
M
1

(l1 ,l2 ,...,lM )

1

2

M

多重母集団寸法指標のノンパラメトリック最尤推定

M




·

(l1 ,l2 ,...,lM )

i=1

k

429

(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,l )Mli Cl
iと表される．ただし，k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) は母集団においてサイズの組 (l1 , l2 , . . . , lM ) のセ
ルのうち，標本においてサイズの組 (l1 , l2 , . . . , lM
) と観測される数を表し，C1 は S から s が生成されるような k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) の組についてのすべての組み合わせを表す．なお，（2.5）式の積において，母集団におけるサイズの組み合わせ (l1 , l2 , . . . , lM ) は (0, 0, . . . , 0) を含ま
ず，標本におけるサイズの組み合わせ (l1 , l2 , . . . , lM
) は (0, 0, . . . , 0) を含む．
時点の数 M ，母集団 Di の大きさ Ni と最大サイズ Li が増加するにつれ C1 の組み合わせ（2002, 2003）と同様の手順の数は膨大になり，時間的に（2.5）式の計算は困難となるため，佐井によりポアソン分布の確率関数の積で近似を行う．
まず，各時点において独立に，しかも母集団 Di (i = 1, 2, . . . , M ) から各個体を他の個体とは独立に確率 λi でサンプリングするベルヌーイ抽出を考える．ベルヌーイ抽出では大きさ ni の標本が得られていても，母集団の大きさは ni 以上のすべての可能性を考えなければならない
が，母集団の大きさを Ni に限定すると，観測された標本の大きさが ni のときに
（2.5）式に対応する尤度関数は(2.6)

LB (S | s)

=





S(l1 ,l2 ,...,lM ) !
 (l1 ,l2 ,...,lM ) k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) !
Mk(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,l )


M
li
li −lili Cl λi (1 − λi )iC1 (l1 ,l2 ,...,lM )·(l1 ,l2 ,...,lM )=M

Ni −niiλn·i (1 − λi )

i=1

i=1









C1 (l1 ,l2 ,...,lM )·(l1 ,l2 ,...,lM )=M
S(l1 ,l2 ,...,lM ) !
k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) !
k(l ,l ,...,l ),(l ,l ,...,l ) (l1 ,l2 ,...,lM )

M


i=1

1

2

M

1

2

M

li Cl

ini
Ni −ni· L(S | s)
NiCni λi (1 − λi )i=1と書ける．尤度関数の添え字の B はベルヌーイ抽出を意味する．多重母集団寸法指標によらず，ベルヌーイ抽出の尤度関数は非復元単純無作為抽出の尤度関数の定数倍となる．従って，各時点において非復元単純無作為抽出された標本に対して最大尤度をとる多重母集団寸法指標
は，ベルヌーイ抽出を想定した場合に，第 i 時点での母集団の大きさが Ni (i = 1, 2, . . . , M ) となる多重母集団寸法指標の中で最大尤度をとるものに一致する．
ベルヌーイ抽出では母集団におけるサイズの組 (l1 , l2 , . . . , lM ) についてのすべての k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) の組 k(l1 ,l2 ,...,lM ) が，他のサイズの組とは独立に多項分布に従い，その確率関数は(2.7)

f (k(l1 ,l2 ,...,lM ) )=S(l1 ,l2 ,...,lM ) !
k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) !
Mk(l ,l ,...,l ),(l ,l ,...,l )
1 2
M


1 2
M
li −lilili Cl λi (1 − λi )i(l1 ,l2 ,...,lM )·(l1 ,l2 ,...,lM )

i=1

統計数理

430

第 57 巻

第 2 号 2009

である．積において，標本におけるサイズの組み合わせ (l1 , l2 , . . . , lM
) は (0, 0, . . . , 0) を含む．
ここで，抽出率 λi (i = 1, 2, . . . , M ) が十分小さいと仮定する．このとき，k(l1 ,l2 ,...,lM ) は，k(l1 ,l2 ,...,lM ),(0,0,...,0) を除くと各要素の期待値は小さく，互いの共分散も小さい．そこでk(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) ((l1 , l2 , . . . , lM
) = (0, 0, . . . , 0)) を，それぞれ独立なポアソン分布の確率関数で(2.8) f (k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) )−{S(l ,l ,...,l

=e


1

2

M)

·

M

· S(l1 ,l2 ,...,lM ) ·li(1−λi )li −li }lili −lii=1 li Cl λi
i

M

i=1

li Cl λii(1 − λi )k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,l )Mk(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) !

と近似すると，ポアソン分布の再生性から k(l1 ,l2 ,...,lM ),(l1 ,l2 ,...,lM ) の (l1 , l2 , . . . , lM ) に関する和である多重標本寸法指標 s(l1 ,l2 ,...,lM ) も，他の多重標本寸法指標とは独立に次のようにポアソ
ン分布に従う．

(2.9)

f (s(l1 ,l2 ,...,lM ) ) =e−µ(l ,l ,...,l )
1 2
M

1

µ(l1 ,l2 ,...,lM ) =2

Ms(l1 ,l2 ,...,lM ) !

ただし(2.10)s(l ,l ,...,l )µ(l 1,l 2,...,l M)(l1 ,l2 ,...,lM )(≥(l1 ,l2 ,...,lM ))S(l1 ,l2 ,...,lM ) ·

M

i=1

.

li
li −lili Cl λi (1 − λi )iである．したがって尤度関数（2.5）式は(2.11)Lap (S | s)
= M

1

ni
Ni −nii=1 Ni Cni λi (1 − λi )·(l1 ,l2 ,...,lM )e−µ(l ,l ,...,l )1 2Ms(l ,l ,...,l )µ(l 1,l 2,...,l M)1

2Ms(l1 ,l2 ,...,lM ) !


と近似される．積において，標本におけるサイズの組み合わせ (l1 , l2 , . . . , lM
) は (0, 0, . . . , 0) を
含まない．

3.

多重母集団寸法指標の推定方法

3.1 アメリカのセンサスデータ
3 節では，アメリカにおいて 1990 年と 2000 年に実施されたセンサスの 1％抽出個票データ（U.S. Census Bureau, 1993, 2003）
を用いて数値実験を行う．数値実験では推定を安定させるために，多重母集団寸法指標の推定値に制約条件を課す検討を行う．
ワシントン州に在住している 20 歳以上の就業者の個人レコードを用いるが，個票データに含まれる項目のうち 2 回の調査において共通の 10 項目（年齢，実子の有無，血縁の子の有無，
性別，結婚，通勤手段，通勤時間，職業，労働週数，週労働時間）
をキー変数として選ぶ．この.
場合，総セル数はキー変数の取り得るすべての値の組み合わせとして K = 2.739 · 1011 である．
しかし，この中には電車通勤で年間の労働週数が 0 週のようなあり得ない組み合わせが含まれている．また，例えば週労働時間としては，記録可能な 0 時間から 99 時間（正確には 99 時間
以上を意味する）までを用いて計算している．したがって，上記の値は名目上の総セル数であ多重母集団寸法指標のノンパラメトリック最尤推定

431

る．なお，このデータには既に秘匿措置が施されており，ここで行う実験は提案する推定方法の評価を目的とするもので，この個票データのリスク評価を目的とするものではない．
これらの 1％抽出個票データは標本調査で得られていると見なすことが可能で，計算された
多重寸法指標を基にして母集団の多重寸法指標を推定することができる．しかし，実際の多重母集団寸法指標は得られないため推定結果を評価することができない．そこで数値実験では，これらの個票データを母集団と見なし，1990 年のデータを D1 ，2000 年のデータを D2 とする．
個体数はそれぞれ N1 = 24846，N2 = 30234 である．
表 1 に多重母集団寸法指標 S を示す．例えば S(0,1) = 21971，S(1,0) = 18154，S(1,1) = 1008のように，サイズ l1 を縦に，l2 を横にとっている．サイズの最大値は 1990 年が 19，2000 年
が 17 であるが，推定では最大サイズを L1 = L2 = 19 と大きい方に合わせる．なお実際の推定
では母集団における最大サイズ L は未知であるため，サイズの最大値として十分大きい値を設定するべきであるが，計算時間上の問題もあり，推定結果に影響を与えない程度の比較的小さい値を選ぶことになろう．
母集団と見なした個票データ Di (i = 1, 2) から，それぞれ抽出率 λi = 1/2 で標本 di を非復元単純無作為抽出した．n1 = 12423，n2 = 15117 である．抽出率 1/2 は十分小さいとは言えず，（2.8）によるポアソン近似の精度に問題がある．しかし，本数値実験では母集団が小さいため，抽出率を小さくとった場合には多重標本寸法指標がほとんど情報を持たなくなってしまうので，
ここでは抽出率としては 1/2 を用いた．
この標本から求められた多重標本寸法指標 s を表 2 に示す．
表 1，2 において，行和と列和はそれぞれ 1990 年と 2000 年の寸法指標を表すが， 1 時点のデータからはサイズ 0 の頻度は観測できない．例えば表 1 において，行和の値 23204 は 1990
年にサイズ 0 のセル数であるが，2000 年のデータがなければ観測することはできない．なお，
これらの値は他の時点において個体が入っていたセルのみの数である．
3.2 多重母集団寸法指標の推定値に対する制約条件とペナルティー関数の導入
多重標本寸法指標を基に，（2.11）式を最大とする多重母集団寸法指標の近似的な最尤推定値を求めるが，推定を安定させるために，寸法指標の推定と同様に次のような制約条件の一部，
表 1.

多重母集団寸法指標（ワシントン州の 20 歳以上の就業者）．出典：U.S. Census Bureau（1993, 2003）．

432

統計数理
表 2.

第 57 巻

第 2 号 2009

多重標本寸法指標（抽出率：1/2）．

またはすべてを多重母集団寸法指標の推定値に課すことを検討する．
（A） S(l1 ,l2 ) ≥ 0 （l1 = 0, 1, . . . , L1 , l2 = 0, 1, . . . , L2 ）（B） (l1 ,l2 ) l1 · S(l1 ,l2 ) + (l1 ,l2 ) l2 · S(l1 ,l2 ) = N1 + N2
（C） S(l1 ,l2 ) ≤ S(l1 −1,l2 ) （l1 = 1, 2, . . . , L1 , l2 = 0, 1, . . . , L2 ），S(l1 ,l2 ) ≤ S(l1 ,l2 −1) （l1 = 0, 1, . . . , L1 , l2 = 1, 2, . . . , L2 ）（D） 2 · log S(l1 ,l2 ) ≤ log S(l1 −1,l2 ) + log S(l1 +1,l2 ) （l1 = 1, 2, . . . , L1 − 1, l2 = 0, 1, . . . , L2 ）,2 · log S(l1 ,l2 ) ≤ log S(l1 ,l2 −1) + log S(l1 ,l2 +1) （l1 = 0, 1, . . . , L1 , l2 = 1, 2, . . . , L2 − 1）
制約条件（A）は多重母集団寸法指標の推定値が非負を表し，（B）は 2 時点での母集団の大きさの合計に関する制約条件で，課すことは妥当であろう．また（C），（D）については， 2.1 節で述べた母集団寸法指標の推定値に対する制約条件（c），（d）を，各行と各列に適用したもので自然な拡張と考えられる．
なお，すべての制約条件で S(0,0) を除いて考える．また，制約条件（D）については S(l1 ,l2 ) ≥ 2についてのみ適用する．
ここで 佐井（2006）に倣って，渋谷（2005）が提案した次のような滑らかなペナルティー関数を導入する．
(3.1)

x

Pnlt(x; ε) = ε log(e ε + 1).

ただし，ε はパラメータである．ペナルティー関数の導入により制約条件を尤度関数に取り込めるため，計算時間の短縮を図ることができる．また，このペナルティー関数はすべての x
について微分可能であるため，後で述べる勾配法による探索も容易である．（3.1）式の微分は
x

(3.2)

dPnlt(x; ε)
eε
= x
dx
eε + 1となり，ロジスティック分布の分布関数を横軸方向に ε 倍だけ拡大したものとなる．
図 1 に ε = 0.01, 0.001 としたペナルティー関数 y = Pnlt(x; ε) のグラフを示す．ε → 0 のとき，グラフは y = 0 (x < 0), y = x (x ≥ 0) に近づくため，ε は原点付近での滑らかさを表すパラメータと考えることができる．

多重母集団寸法指標のノンパラメトリック最尤推定

図 1.

433

x

滑らかなペナルティー関数 Pnlt(x; ε) = ε log(e ε + 1)．

（B）を除く制約条件をペナルティー関数で表し，次のように対数尤度関数へ取り込む．
(3.3)

F(S | s) = log Lap (S | s) − c1
− c2





Pnlt(−S(l1 ,l2 ) ; ε1 )(l1 ,l2 )Pnlt(−(S(l1 −1,l2 ) − S(l1 ,l2 ) ); ε2 )(l1 ,l2 )− c2

Pnlt(−(S(l1 ,l2 −1) − S(l1 ,l2 ) ); ε2 )(l1 ,l2 )− c3

Pnlt(−(log S(l1 −1,l2 ) + log S(l1 +1,l2 ) − 2 · log S(l1 ,l2 ) ); ε3 )(l1 ,l2 )− c3

Pnlt(−(log S(l1 ,l2 −1) + log S(l1 ,l2 +1) − 2 · log S(l1 ,l2 ) ); ε3 ).

(l1 ,l2 )

（A）
，
（C）
，
（D）
の制約条件に関するペナルティー関数において，ε1 ，ε2 ，ε3 はそれぞれの滑らかさを表すパラメータであり，c1 ，c2 ，c3 はそれぞれのペナルティー関数の関数 F への影響力を表す係数である．
関数 F を S(0,0) を除く各 S(l1 ,l2 ) (l1 = 0, 1, . . . , L1 ; l2 = 0, 1, . . . , L2 ) で偏微分し，0 と置いた連立方程式を解くことによって関数 F の最大値を得ることが最善であるが，解析的に解くことは困難である．
そこで適当な初期値から，ステップごとに，F が増加する方向に各 S(l1 ,l2 ) を独立に移動させる勾配法によって最大値を探索する．なお S(l1 ,l2 ) として実数の範囲を考える．
Ŝ(l1 ,l2 )(p) を p 番目のステップにおける多重母集団寸法指標の推定値とするとき，∂FŜ(l1 ,l2 )(p+1) = Ŝ(l1 ,l2 )(p) + t(l1 ,l2 )(p) · g(3.4)∂ Ŝ(l1 ,l2 )(p)によって，p + 1 番目のステップの推定値 Ŝ(l1 ,l2 )(p+1) を求める．ただし，g(·) は微係数を調整する関数である．また t(l1 ,l2 )(p) は移動する距離を調整する係数で，(l1 , l2 ) と Ŝ(l1 ,l2 )(p) などを基に各ステップで変化させる．
なお，このようにして多重母集団寸法指標の推定値を変化させた場合，一般に p + 1 番目の推

定値は制約条件（B）が表す超平面 (l1 ,l2 ) l1 · S(l1 ,l2 ) + (l1 ,l2 ) l2 · S(l1 ,l2 ) = N1 + N2 上では得られないが，推定値と原点とを通る直線とこの超平面との交点へ移動させる補正をステップごとに行
う．母集団の大きさに関する制約条件としては， (l1 ,l2 ) l1 · S(l1 ,l2 ) = N1 ， (l1 ,l2 ) l2 · S(l1 ,l2 ) = N2

434

統計数理

第 57 巻

第 2 号 2009という 2 つの条件を用いるべきであるが，上記のような補正が容易ではないため，ここではやや緩い制約条件としている．
3.3 アメリカのセンサスデータについての推定結果
この節では，多重母集団寸法指標の推定を種々の設定で行い，それらの結果の比較を行う．
数値実験では，原則として（3.3）式で， c1 = 1.0，c3 = 0.0，ε1 = 0.001，ε2 = ε3 = 0.01 とし，
c2 については 0, 0.1, 1.0, 10.0, 100.0, 1000.0 の 6 通りを用いる．これに加えて，c2 = 1000.0，
c3 = 10.0 の場合も考える．また（3.4）式において
(3.5)
(3.6)

g(x) = sign(x)· | x |0.5 ,t(l1 ,l2 )(p) = b(p) · Ŝ(l1 ,l2 )(p) / (l1 + l2 )0.5

として，p 番目のステップにおける多重母集団寸法指標の推定値が大きいほど大きく変化させ，サイズの和 l1 + l2 が大きいほど小さく変化させる．また，b(p) は 10−4 を初期値として，3 ス
テップ連続して F が増加しなかった場合には 0.9 倍し，F が増加した場合に 1.1 倍し，10−7を下回った時点で探索を終了する．
結果を表 3 ∼ 9 に示す．なお，推定値と周辺頻度である行和，列和は，それぞれ整数に四捨五入しているため，合計が一致していない場合がある．また表中の 0 は実際には 0.5 未満の正の値である．
表 3 は制約条件として
（A）
（
，B）
のみを課した場合の結果である．母集団寸法指標の推定において，2.1 節の制約条件
（a）
（
，b）
のみを課した場合の結果と同様に，多くの推定値が 0 となる一方で，いくつかのサイズの組み合わせについては大きな値が得られている．
表 4 ∼ 8 では，制約条件（C）の影響を徐々に強くした推定結果を示しており，表 9 では制約条件（D）も追加している．表 4（c2 = 0.1），表 5（c2 = 1.0）の推定結果が実際の多重母集団寸法指標に近いが，c2 の値によって結果は大きく異なっている．
3.4 推定方法の改善 ∼ 周辺母集団寸法指標を先行して推定する 2 段階法
残念ながら多重寸法指標の推定に関する経験はそれほど多いとは言えず，適切な制約条件，
係数，パラメータを選択できる保証はない．これに対して，寸法指標の推定に関しては経験的知識の蓄積があり，ある程度安定した推定結果を得ることが可能である．多重母集団寸法指標の表 3. 多重母集団寸法指標の推定値（c2 = 0.0，c3 = 0.0）．

多重母集団寸法指標のノンパラメトリック最尤推定
表 4. 多重母集団寸法指標の推定値（c2 = 0.1，c3 = 0.0）．

表 5. 多重母集団寸法指標の推定値（c2 = 1.0，c3 = 0.0）．

表 6.

多重母集団寸法指標の推定値（c2 = 10.0，c3 = 0.0）．

435

統計数理

436

第 57 巻

第 2 号 2009表 7.

多重母集団寸法指標の推定値（c2 = 100.0，c3 = 0.0）．

表 8.

多重母集団寸法指標の推定値（c2 = 1000.0，c3 = 0.0）．

表 9.

多重母集団寸法指標の推定値（c2 = 1000.0，c3 = 10.0）．

多重母集団寸法指標のノンパラメトリック最尤推定437

サイズ 0 を除く行和と列和は，それぞれ 1990 年と 2000 年の母集団寸法指標であるが，表 3∼9で推定されたこれらの値も c2 の値によって大きく異なっている．
そこでこの節では，推定の改善のために，まず各年の母集団寸法指標をそれぞれ推定した後に，その推定値によって行和と列和を固定しながら多重母集団寸法指標を推定する 2 段階の推定方法を提案する．
母集団寸法指標の推定では， 3.2 節と同様に，2.1 節で紹介した制約条件
（a），
（c），
（d）をペナルティー関数で表して，（3.7）式のように対数尤度関数に取り込んだうえで，勾配法による最大値の探索を行う．なお，探索の各ステップの推定値は一般に制約条件（b）を満たさないため，
3.2 節と同様に，推定値と原点とを通る直線と（b）が表す超平面との交点へ移動させる補正をス
テップごとに行う．
(3.7)

F(S1 , S2 , . . . , SL | s1 , s2 , . . . , sL )= log Lap (S1 , S2 , . . . , SL | s1 , s2 , . . . , sL )
− c1

L


Pnlt(−Sl ; ε1 ) − c2

l=1

− c3

L−1


L


Pnlt(−(Sl−1 − Sl ); ε2 )l=2Pnlt(−(log Sl−1 + log Sl+1 − 2 · log Sl ); ε3 ).

l=2c1 ，c2 ，c3 は係数，ε1 ，ε2 ，ε3 はパラメータであるが，これまで扱った母集団寸法指標の推定の多くでは， c1 = c2 = 10.0，c3 = 1.0，ε1 = ε2 = 0.0001， ε3 = 0.001 を標準の組み合わせとしてきた．ここでは，この組み合わせを含む c2 ，c3 のいくつかの組み合わせに対して母集団寸法指標の推定を行った．
1990 年の推定結果を表 10 に，2000 年の推定結果を表 11 にそれぞれ示す．第 1 行がサイズを示し，その下に推定された母集団寸法指標を表示している．例えば S1 が表 10 では 19642 ∼ 19782，
表 11 では 23367 ∼ 23405 と推定されており，係数 c2 ，c3 が変化しても安定した推定結果が得られている．
上記で標準と呼んだパラメータの組み合わせを用いて推定された 1990 年と 2000 年の母集団寸法指標を，それぞれ Ŝl1 ，Ŝl2 とするとき，多重母集団寸法指標の最尤推定値の勾配法による
表 10.

c2 ，c3 の組み合わせに対する母集団寸法指標の推定値の比較（1990 年）．

表 11.

c2 ，c3 の組み合わせに対する母集団寸法指標の推定値の比較（2000 年）．

統計数理

438
表 12.

第 57 巻

第 2 号 2009

周辺頻度を固定した場合の多重母集団寸法指標の推定値（c2 = 1000.0，c3 = 10.0）．

探索における各ステップで，制約条件（B）に変えて
(B )

L2

l2 =0

S(l1 ,l2 ) = Ŝl1 (l1 = 1, 2, . . . , L1 ),L1

S(l1 ,l2 ) = Ŝl2 (l2 = 1, 2, . . . , L2 )l1 =0を用いた補正を行う．ここでは上式の左辺と右辺との差が 0.01 以下になるまで，各ステッ
プにおいて補正を繰り返す．前節の検討で良い推定結果が得られなかった係数の組み合わせ
c2 = 1000.0， c3 = 10.0 を用いた場合の，多重母集団寸法指標の推定結果を表 12 に示す．制約条件（B’）による補正において誤差を許しているため，表 12 の周辺頻度は表 10, 11 の値とは一致していない場合がある．
推定結果は改善されている．この方法を用いることにより，（3.3）式におけるパラメータと係数の決め方によらずに，ある程度，安定した推定が可能と考えられる．
4.

おわりに

本論文では，継続調査で得られた複数の時点の個票データに対する同時リスク評価を行うた
めに，寸法指標の概念を多重寸法指標に拡張し，多重標本寸法指標に基づく多重母集団寸法指標の制約付きノンパラメトリック推定法を提案した．
多重母集団寸法指標の推定結果で，例えば S(1,1) の値が相対的に大きいことがわかれば，第1 時点でサイズ 1 のセルを減らす秘匿措置を行うと，第 2 時点においてもサイズ 1 のセルを減らす効果の高いことがわかる．逆に各時点では S1 の値が大きいにも関わらず S(1,1) の値が小さければ，各時点の個票データに対して独立に秘匿措置を行った方が効果的であろう．
これまで扱った官庁統計などの個票データでは，個体数と比較してセル数が極めて多い場合がほとんどである．寸法指標では観測サイズ 0 のセルは観測できないものの，サイズ 1 以上のセル数と比較して数多く存在し，多重寸法指標ではその一部が S(l1 ,0) ，S(0,l2 ) として観測される．S(1,1) の評価では各時点の S1 との比較が基本と考えられるが，多重寸法指標内で比較する場合には S(l1 ,0) ，S(0,l2 ) ではなく S(l1 ,l2 )（l1 , l2 = 0）と比較するべきと考える．
本論文では，アメリカにおいて 1990 年と 2000 年に実施されたセンサスの 1％抽出個票デー
タを題材として，提案した推定方法を適用し，その有効性について検討も行った．多重母集団
寸法指標の推定についてはこれまでの経験が少ないことを勘案し，初めに多重母集団寸法指標

表 13.

多重母集団寸法指標（ワシントン州の 20 歳以上のすべての個人，就業と年間収入の項目を含む）．出典：U.S. Census Bureau（1993, 2003）．

多重母集団寸法指標のノンパラメトリック最尤推定

439

統計数理

440

第 57 巻

第 2 号 2009の行和と列和である母集団寸法指標を推定し，次にその値を制約条件として多重母集団寸法指標を推定するという 2 段階の推定法も提案した．
これらの数値実験は抽出率が 1/2 と高かったこともあり，結果は概ね満足できるものであったが，一般の官庁統計で用いられている 1/1000 程度の抽出率で得られたデータへ適用するためには，多重母集団寸法指標の推定値に対する制約を強くするなど，更なる検討が必要と考えられる．ただし，官庁統計では一般に標本の大きさもここでの実験よりは大きい．また，ここではキー変数の値についてトップコーディングや丸めなどの秘匿措置を行わなかったが，通常は公開のリスクを小さくするため秘匿措置を行うことになる．その場合，多重母集団寸法指標の各時点のサイズの最大値は大きくなり，多重標本寸法指標の各時点のサイズの最大値もまた大きくなる．母集団寸法指標の推定では，標本寸法指標のサイズの最大値が大きいほど安定した推定を行うことができることが経験的に知られている．そのため抽出率の減少に比べて精度の低下は小さい可能性もある．
多重寸法指標の形状についても多様なパターンが存在する．同じセンサスの 1％抽出個票デー
タのうち，ワシントン州に在住している 20 歳以上の全個人レコードについて，3.1 節で用いた
10 項目に，就業，年間収入の 2 項目を追加して多重寸法指標を求めると表 13 のようになる．
上部と左部に延びている非負の頻度のセルの多くには，非就業者が入っている．このような形状の母集団寸法指標を推定するためには，（E） l1 ≥ l2 のとき S(l1 ,l2 ) ≤ S(l1 +1,l2 −1) （l1 + l2 = 2, 3, . . . , min(L1 , L2 )），l1 ≤ l2 のとき S(l1 ,l2 ) ≤ S(l1 −1,l2 +1)（l1 + l2 = 2, 3, . . . , min(L1 , L2 )）（F） 2 · log S(l1 ,l2 ) ≤ log S(l1 −1,l2 +1) + log S(l1 +1,l2 −1)（l1 + l2 = 2, 3, . . . , min(L1 , L2 )）のような制約条件を追加しなければ，良い推定結果は得られない．多重母集団寸法指標の推定値に課す制約条件についても，今後検討していく必要があろう．
なお，本論文では M 時点の多重母集団寸法指標に対する最尤推定法を提案したが，数値実
験では 2 時点の個票データしか扱っていない．時点数が多い場合には，多重母集団寸法指標の推定値に対して，より強い制約条件を用いなければ，推定が不安定になると想像される．この点についても今後の課題になる．特に，相関構造を表現できるような制約条件が必要となるかもしれない．
また時点数が多い場合には，推定目標を多重母集団寸法指標全体ではなく，S(1,1,...,1) やS(1,0,...,0) などに絞ることや，各時点におけるサイズが 2 以下の個体数のように，特定の範囲に含まれるサイズの組み合わせに対する個体数を推定目標とすることも必要になろう．
謝辞

本論文で扱った多重寸法指標に関して，慶應義塾大学名誉教授の渋谷政昭先生からは数々の有益な助言をいただきました．また，金沢大学の星野伸明先生，東京大学の竹村彰通先生からは本研究に対するアドバイスをいただきました．ここに感謝いたします．2 名の査読者と編集委員の方々からも貴重なご意見をいただきました．ここに感謝いたします．本論文は科学研究費補助金
（課題番号 19300098, 18200019），統計数理研究所共同利用研究プログラム（20-共研-2029）の研究成果に基づくものである．

多重母集団寸法指標のノンパラメトリック最尤推定441

参 考 文 献馬場康維，坂口尚文（2006）. 複数名簿のマッチングによる共通集合の推定，2006 年度統計関連学会連合大会予稿集，p. 117.
. Disclosure control of microdata, Journal
Bethlehem, J. G., Keller, W. J. and Pannekoek, J.（1990）
of the American Statistical Association, 85, 38–45.
Pitman, J.（1995）. Exchangeable and partially exchangeable random partitions, Probability Theory
and Related Fields, 102, 145–148.
佐井至道（2002）. サイズインデックスの制約付き最尤推定，岡山商大論叢, 37（3）, 61–79.
佐井至道（2003）. 母集団寸法指標のノンパラメトリック推定，統計数理, 51（2）, 183–198.
佐井至道（2006）. ペナルティ関数を利用した母集団寸法指標の制約付きノンパラメトリック推定，岡山
商大論叢, 42（1）
, 1–21.
佐井至道（2007）. 多重寸法指標のノンパラメトリック推定，2007 年度統計関連学会連合大会予稿集，186–186.
佐井至道（2008）. 母集団多重寸法指標のノンパラメトリック推定，岡山商大論叢, 43（3）, 1–18.
渋谷政昭（2005）. 滑らかな罰金関数，日本学術振興会科学研究費補助金研究会資料，1–3.
渋谷政昭，佐井至道（2007）. 多重指標の確率分割，応用統計学会 2007 年度年会講演予稿集，13–18.
Sibuya, M. and Sai, S.（2008）. Analysis of a dataset for statistical disclosure control by random
partition of a multi-index, Cherry Bud Workshop 2008, 1–13.
. 1990 Census of Population and Housing, Public Use Microdata Samples
U. S. Census Bureau（1993）（microdata）, Washington, D. C.
U. S. Census Bureau（2003）. Census 2000 Public Use Microdata Sample（microdata）, Washington,
D. C.

442

Proceedings of the Institute of Statistical Mathematics Vol. 57, No. 2, 425–442(2009)

Nonparametric Maximum Likelihood Estimationof Multi-population Size Indices—Application to Microdata Sets of Two Occasions—Shido Sai
Faculty of Economics, Okayama Shoka University

Population size indices estimated from observed sample size indices are often used toassess the disclosure risk of a microdata set sampled from a population. Parametric methods with superpopulation models such as the Poisson gamma model or the Pitman model
are usually used to estimate population size indices. A restricted nonparametric maximum
likelihood estimation method has also been proposed, and its problem of computing time
is being improved.
Though most sample surveys of oﬃcial statistics are repeated periodically, it is usual
to carry out risk assessment for each microdata set independently. In order to assess thedisclosure risk for several microdata sets simultaneously, the size indices are extended tothose of the multi-indices, and a nonparametric method for estimating multi-population
size indices from multi-sample size indices is proposed.
Restrictions for the estimates of multi-population size indices are expressed as penalty
functions and are taken into a log likelihood function. To improve the estimation, a twostage method is also proposed, in which the multi-population size indices are estimatedunder the condition that their marginal frequencies are ﬁxed on the population size indicesestimated in advance. These proposed methods are applied to 1990 and 2000 1-Percent
Public Use Microdata Sample Files in the U. S.

Key words: Multi-indices, size indices, microdata, risk assessment, nonparametric estimation, oﬃcialstatistics.