統計数理（2003）
第 51 巻 第 2 号 297–319
c 2003 統計数理研究所


特集「個票開示問題の統計理論」［総合報告］

超母集団モデルによる個票開示リスク評価
星野 伸明

†

（受付 2003 年 1 月 31 日；改訂 2003 年 4 月 18 日）

要旨

個票データは，各個体が分割表のどのセルに所属するかを示すものである．そして秘匿処理
は，分割表の解像度を粗くする操作に他ならない．分割表を細かくすればセル内の個体数（頻度）が減る事を考えると，頻度の頻度（寸法指標）推測は開示リスク評価において，重要な役割を果たす．しかし母集団に何の仮定もおかないと，実用的な推測は不可能である．このため有限母集団解析の定石どおり，超母集団モデルを仮定する．経験的に母集団は Zipf の法則に従う
ので，右裾が長い分布で混合したポアソン分布を利用する．つまり基本モデルでは，各セルでの頻度を独立同一な混合ポアソン分布とみなす．個体総数（の期待値）
を一定のままセル数を無限大とする極限操作を，小数法則と呼ぼう．無限分解可能な混合ポアソン分布の基本モデルに小数法則を適用すると，モデルとして意味のある極限分布が得られる．無限分解可能な混合ポ
アソン分布は多いので，この法則によって便利なモデルを新規に導出できるかもしれない．な
お幅広い母集団の記述の為には，多くのモデルを使い分ける必要が有る．従って新規モデルを提供していく事は重要である．
キーワード： 母集団一意，無限分解可能，複合ポアソン，混合ポアソン，自然数の確
率分割．
1.

導入

個票が特定の調査客体のものだと明らかになる事を，「個体識別」または「個体開示」という．「（個票）開示リスク」とは，そのような危険性の事である．この概念自体は自明かもしれないが，リスク評価の統計学的定式化は必ずしも自明ではない．本稿では超母集団モデルという方法論が，個票開示リスク評価で果たす役割を考察する．本論文の構成は以下のとおりである．1.1 節では計測可能なリスク概念が導入される．1.2 節ではリスク評価における超母集団モデルの必要性を説明する．1.3 節では既存のモデル同士の関係を明らかにする．その中で無限分解可能性と条件付けが，重要な意味を持つ．2 章の各節では，既存のモデルの使用法が解説される．3 章では今後の課題について述べたい．
1.1 個票開示リスクの定式化
ある統計調査で，個体の同定・識別に役立つ情報を含む項目が L 個有るとしよう．そのような項目を「キー変数」という．調査項目の中には，個体の識別に使えないものも存在する．
例えば事業所の特定品目在庫などは，外からはうかがい知れないだろう．リスク評価について
は，そのような項目は無視して良い．第 l, l = 1, 2, . . . , L, 変数について，レンジは cl 個のカテ†

金沢大学 経済学部：〒920–1192 石川県金沢市角間町

298

統計数理

第 51 巻

第 2 号 2003

ゴリーに分類されるだろう．例えば性別という変数なら，男と女の二分類 (c = 2) になる．個票データセットは第 l キー変数が第 l フィールドを構成し，調査客体（統計単位）
毎にレコードをつくる．調査客体は通常個人や世帯，事業所である．
ではどのようなデータセットが危険と考えられるか．極端なケースとして，名前や住所を含
むようなデータセットが挙げられる．この場合調査客体の特定は容易であり，このようなリス
クは明らかに許容する事が出来ない．何故「名前」や「住所」をデータとして含んではならないのかというと，特定の名前や住所という条件を満たす母集団での個体が非常に少ない（たい
てい一意に定まる）からである．小さいグループの中から個体を特定するのは，比較的容易かもしれない．一般に考えて，キー変数の特定の組み合わせ条件を満たす個体が少なければ，それらは比較的危険な個体と考えられる．キー変数について組み合わせの総数は各カテゴリー数の積，すなわち c1 × · · · × cL であり，これを J と書く．この組み合わせのそれぞれを「セル」と呼び，各個体はいずれかのセルに所属する．母集団において第 j セルに所属する個体数を Fjと書く（j = 1, 2, . . . , J ）．ここで Fj が小さければ小さい程，そのセルに所属する個体データの公開はより危険であろう．
直感的に考えて，調査客体の情報を粗くすればそのデータセットはより安全になる．そのような秘匿処理（修正）については，いくつか選択肢が有る．任意の変数の組について，全てのレ
コードに同じ修正をする場合「大域的」という．そうでなければ「局所的」と呼ぶ．以下では実務的に重要な，大域的修正のみ考える．また修正の技術としては，「隠蔽」と「再符号化」を考える．前者は変数の組を，公開対象から外してしまう．「名前」などは大域的隠蔽の自明な対象である．再符号化を説明しよう．これは変数のカテゴリーを再編成する事を言う．「学歴」
を例にとる．「大卒」と「大学院卒」が別物として調査されていても，再符号化して両者を同一の「大卒以上」で扱える．また再符号化には，数値の丸め等も含まれる．隠蔽はその変数のカテゴリー数を 1 にする事と解釈出来るので，再符号化の極端なケースと同等である．なお再符号化を進める事で，危険性は単調に減少するはずである．セルの大きさをリスクの指標とす
れば，これをうまく説明する事が出来る．すなわち再符号化によるカテゴリー数 cl の減少は，
J の減少を意味する．この時，セルあたりの個体数は増える方向で変化する．
以上の議論から，セルの大きさはリスクの指標としてふさわしい性質を持っている事が理解される．従って考え方としては，秘匿処理を用いて分割表を確定，この表の危険性をセルの大きさで計る，という事になる．特定のセルの大きさではなく危険なセルの総数に興味が有るので，「頻度の頻度」（Good（1953））を数える．すなわち

J(定義 1.1)

Si =I(Fj = i),i = 0, 1, . . . ,j=1

ただし I(·) は指示関数であり，I(Fj = i) =

1, Fj = i,0, Fj = i.

なお Fj , j = 1, 2, . . . , J, は非負整数である．以下では Sibuya（1993）にならい，Si を「寸法指標」と呼ぶ．これは個票開示リスク評価の分野では，Greenberg and Zayatz（1992）が言う equivalenceclass の概念に対応する．ここで記法をまとめておこう．空でないセルの総数を∞(定義 1.2)U=

i=1

Si = J − S 0299

超母集団モデルによる個票開示リスク評価

であらわす．また母集団サイズは

 J(定義 1.3)

N=

∞

Fj =

j=1

iSi < ∞i=1

とする．標本については，セル内個体数を fj と書く．標本寸法指標も同様に定義され，si で表す．u は標本において空でないセルの総数である．また標本サイズは n = Jj=1 fj としよう．
母集団で一意な事が知られている個体のレコードが公開されると，理屈としては個体識別が可能である．故に S1 は「母集団一意」と呼ばれ，Marsh et al.（1991）によればリスクの指標として，米国・英国の統計当局が用いている．また，二意な事が知られている個体のレコードが公開されたとしよう．このレコードで表される属性を持つ二者は，自分が公開されているのでな
ければ他方が公開されている事が分かる．論理的帰結として識別が可能なケースは，この二つである．しかし母集団で一意や二意な事が知られているような個体は，非常に稀である．従って一意や二意の個体の公開が，即危険という事にはならない．最低限言えるのは，S1 や S2 が大きいようなデータはより危険だという事である．また S3 以下なら即安全とも言い切れない．
このように考えると，ある母集団のリスクを寸法指標の重み付き和で表すのは自然である．
すなわちリスクを N
i=1 wi iSi , ただし wi , i = 1, . . . , N, は非負，のように考えれば計測可能な概念となる．例えば Bethlehem et al.（1990）は



 
N

2

i
N

i=1

Si

の逆数を “resolution” と呼び，リスクの指標として用いている．また Greenberg and Zayatz（1992）では，エントロピーから類推して

  
N

−

log

i=1

i
N

i
Si
N

が用いられた．母集団一意のみ考慮する場合は，w1 = 1, w2 = w3 = · · · = 0 である．なおSkinner et al.（1994）によれば，カナダの統計当局が用いているリスクの指標は，標本一意かつ母集団一意な個体数の，標本一意に占める割合である．これはnS1
N · s1で推定できる．私見を述べれば，重み wi の選択は利便性の観点から行えば良い．そして許容できるリスクの範囲は，重みの付け方に依存して決まる．
1.2 リスクの推測
全数調査については以上の議論で済む．しかし標本調査では母集団の寸法指標が未知であり，推定の対象となる．ここで利用可能な情報は，標本の寸法指標 (s0 , s1 , . . . , sn ) である．単純化のために，標本が単純無作為非復元抽出でとられたとしよう．この場合

(1.1)






E(s1 )
E(s2 )
E(s3 )..
.
E(sn )

=

1
N
n







1
1

N−1
n−1

2
1

N−2
n−1

3
1

N−3
n−1

···

q
1

N−q
n−1

0

2
2

N−2
n−2

3
2

N−3
n−2

···

q
2

N−q
n−2

0

3
3

N−3
n−3

···

q
3

N−q
n−3

0
..
.

..

0

···

..
.

.
q
q

N−q
n−q






S1
S2
S3
..
.
Sq

300

統計数理

第 51 巻

第 2 号 2003という関係が成立し，ここから Si の推定量を構成できる．Engen（1978）の定理 2.1 によれば，もし i ≤ n ならば（標本寸法指標に基づく）
唯一の不偏推定量が存在して，(1.2)

ŜiEngen =

 
n

k

k=i

j=i

N
j

k
j
n
j

j
i

(−1)j−iskと書ける．しかし Si , i > n について，不偏推定量は存在しない．
残念ながら推定量（1.2）は，きわめて不安定な事が知られている（例えば渋谷（2002）の数値実
験を参照せよ）．理由として Shlosser（1981）は，（1.1）の右辺の行列がほとんど特異に近いと指摘する．また渋谷（（2000）, p. 159）は推定が不安定である事について，直感的な説明を与えている．このように母集団に何の仮定もおかない推定がうまくいかない時，超母集団という構造を仮定するのが一つの方法である．すなわち母集団の性質を分布で記述する．以下では我々も，
この方法論を採用しよう．なお超母集団アプローチについては，Smith（1976）の整理が分かりやすい．
なお筆者は超母集団を，「議論を明確にする為の数学的道具に過ぎない」と解釈する．超母集団を，「実世界の確率的機構・過程の記述」とみなすべきではない．Feller（1966）の 2.4 節における，ロジスティック分布関数に関する警告を見よ．多くの相容れない理論モデルが，同じデータによって支持されてきた実例が紹介されている．データが表面的にあてはまるからと言って，背後の構造は一意に特定されない．だとすればモデルの価値は，実機構・過程と切り離されて存在する．なお仮定される超母集団の母数は，データより推定する．すなわち経験ベ
イズ法が採用される．Lehmann and Casella（（1998）, Chapter 4）によれば，経験ベイズ法はモデルの誤特定について比較的頑健である．従って構造の仮定による一般性の喪失を，ある程度緩和する意味がある．とはいえ真の母集団と仮定される超母集団が遠ければ，寸法指標の推定が良いはずがない．観測されている様々な母集団が，頻度的な裏づけとなるような超母集団が望ましい．
実は母集団の寸法指標は，経験的に一定の傾向，右への歪みを示す．例えば単語の使用頻度，
個人所得や都市の人口を，大きさ順に並べるとする．一般に，順序統計量 x(1) ≥ x(2) ≥ · · · ≥ x(n).
で考えよう．ここで r = 1, 2, . . . , n について，ある定数 c が存在して r · x(r) = c と Zipf（1949),
個票開示リスク評価も例外ではないでMandelbrot（1983）等は主張した．様々な分野
このようなデータが観測される為，これは「Zipf の法則」と呼ばれている．この経験則は寸法指標に関して，Si ∝ 1/i2 と主張している（Urzúa（2000））．同様に Pareto（1897）は，所得分布の研究から Si ∝ 1/ia , ただし a は適当な定数，としている．個票の寸法指標実データは，例えば竹村（1998）が与えている．これを見ると，Zipf 等の主張は不当とは言い難い．つまり，このような「逆 J 字」を記述出来る分布がモデルとして望ましい．
「逆 J 字」型分布は，右裾の長い連続分布とポアソン分布を混合して得られる．すなわち，
第 j セルの度数 Fj が平均 λ のポアソン分布に従い，λ がガンマ分布や対数正規分布などに従うと考える．そしてもっとも単純なのは，Fj が独立に同一の分布に従う場合だろう．この時寸法指標の同時分布は


∞

(1.3)

P(S0 = t0 , S1 = t1 , . . . ) = J!

P(Fj = i)ti
ti !
i=0

という多項分布で書ける．従って Fj の分布を規定すれば寸法指標の挙動は定まる．また既知の母集団サイズ N0 を，モデルに取り込む必要がある．ここで二通りの考え方が有りうる．一
つ目は母集団サイズを固定し，N = N0 と制約する．二つ目は母集団サイズを確率変数とみなし，E(N ) = N0 と制約する．一番目の方法はよりリアルだが，モデルが複雑になり条件付分布

301

超母集団モデルによる個票開示リスク評価

P(F1 , F2 , . . . , FJ | N = N0 ) を評価しなければならない．二番目の場合は独立同一分布を扱うので，評価が単純となる．すなわち P(F1 , F2 , . . . , FJ ) = Jj=1 P(Fj )，ただし E(Fj ) = N0 /J, j =1, 2, . . . , J, と制約する．
次に標本分布を考察する．単純化の為に単純無作為非復元抽出，すなわちP(f1 = y1 , f2 = y2 , . . . , fJ = yJ | F1 , F2 , . . . , FJ ) =およびベルヌーイ抽出，すなわち

F2
y2

N
n

· · · FyJJ    JP(f1 = y1 , f2 = y2 , . . . , fJ = yJ | F1 , F2 , . . . , FJ ) =

F1
y1

j=1

Fj
yj

n0
N0

yj

，

N0 − n0N0Fj −yjのみ使用する．なお本稿では，周辺分布 P(f1 , f2 , . . . , fJ ) を「標本分布」と呼ぶ．モデル
P(F1 , F2 , . . . , FJ ) = pN が所与の時，抽出方法に応じて標本分布を得たい．もし標本分布がN0 を既知の標本サイズ n0 に置き換えるだけで得られる，すなわち P(f1 , f2 , . . . , fJ ) = pn と表されるなら議論は簡潔である．言い換えれば，超母集団から直接標本を抽出する場合の分布と，超母集団から抽出された母集団より標本を抽出した場合の分布が同じ状況が望ましい．まず母
集団サイズが固定されたモデルの場合，Takemura（1999）がそのようになる条件を示した．すなわち個体のラベルにモデルが依存しない時，pN から単純無作為非復元抽出された標本分布は
pn となる．また母集団サイズが確率変数の場合，本稿で考察する独立同一混合ポアソン分布モ
デル pN は，ベルヌーイ抽出をすれば標本分布が pn となる．何故なら Feller（（1957）, p. 268）によれば，平均 N0 λ のポアソン分布に従う確率変数から抽出率 n0 /N0 でベルヌーイ抽出をす
れば，標本は平均 n0 λ のポアソン分布に従う．
以上の議論を踏まえると，適当な標本抽出法を仮定すれば，寸法指標の推測は次のように行うことが出来る．まず仮定されたモデルの母数を標本から推測する（θ̂）．そして推測された母数の下で，寸法指標の期待値 Eθ̂ (Si ) を評価する．これを寸法指標の推定値 Ŝi と考える．
1.3 混合ポアソン分布と小数法則使用するモデルが決まれば，モデルに応じてリスクの推定値が得られる．ではいかなるモデルを用いるか．出来るだけ多様な母集団を記述出来る，柔軟なモデルが望ましい．しかし私見
では，単独のモデルで全ての母集団をカバーする必要は無い．多くのモデルを用意し，データに依存して AIC 等の基準でモデルを選択すれば良い．そして選択されたモデルによるリスク
評価を，採用すべきである．従って多くのモデルを提供する事は，記述出来る母集団の範囲を広げる為に重要である．まず既存の超母集団モデルの挙動をふまえた上で，モデル開発の指針を導きたい．以下では個票開示リスク評価における超母集団モデルの文脈を整理する．そしてモデル同士の関係を議論しよう．なお本節ではモデルの名前に言及するにとどめ，その詳細は次節で記述される．
本分野では Bethlehem et al.（1990）がガンマ＝ポアソンモデル（またはポアソン・ガンマモデル）を超母集団モデルとして利用したのが嚆矢となる．これは独立同一な Fj , j = 1, . . . , J, の分布をガンマ分布で混合したポアソン分布，すなわち負の二項分布とみなす．佐井（1998）は Jと E(N ) の大小関係について場合分けをし，ガンマ＝ポアソンモデルの挙動を解析した．また佐井（2000）では，ガンマ＝ポアソンモデルを利用してリスク評価を議論している．しかし，ガ
ンマ＝ポアソンモデルの評判は芳しくない．Zayatz（1991）はガンマ＝ポアソンモデルをコルモ
ゴロフ＝スミルノフ検定にかけたが，有意に当てはまりを欠くと報告している．Skinner et al.
（1994）はガンマ=ポアソンモデルが S1 を過小推定するとし，モデルの限界ではないかと指摘した．小さなセルが支配的なデータセットは，経験的にガンマ＝ポアソンモデルで記述しきれ

302

統計数理

第 51 巻

第 2 号 2003

ないようである．Hoshino（2001）は Pitman モデルの含意として，ガンマ＝ポアソンモデルは.
安全（つまり S1 /U = 0）なデータセットしか記述できないだろう，と述べている．一意な個体割合が無視できないデータセットが興味の対象だとすれば，ガンマ＝ポアソンモデルが不十分という評価を受けているのは納得できる．それからガンマ＝ポアソンモデルを修正する試みを紹介する．確率変数 X が負の二項分布に従う時，Chen and Keller-McNulty（1998）は X + 1 の分布を Slide Negative Binomial と呼ぶ．彼らは Fj の分布を SNB と仮定，S1 の推定が改善されたと報告する．佐井・竹村（2000）では，セル間に相関が有る場合を考察した．
他に当然，ガンマ＝ポアソン以外のモデルも使用されている．Skinner and Holmes（1993）
は，
Fisher et al.（1943）の対数級数分布および対数正規＝ポアソン分布を，米国及びイタリアのセ
ンサスデータへ当てはめている．ここでは対数正規＝ポアソン分布の当てはまりが良好であった．なお Skinner and Holmes（1998）は，対数正規＝ポアソンモデルの母数を対数線形モデルで記述する事を考察している．Hoshino and Takemura（1998）は Anscombe（1950）の対数級数モデルを再発見し，母数推測を議論した．また Takemura（1999）は，ディリクレ＝多項（多変量負の超幾何分布）モデルを提案した．Omori（1999）はディリクレ＝多項モデルを仮定し，母集団一意の事後確率を議論した．Ewens モデルは Takemura（1999）で言及されているが，同分布の使用は Samuels（1998）にも見られる．Hoshino（2001）は対数正規＝ポアソンモデル，Ewens モ
デル，ディリクレ＝多項モデル，Pitman モデルを 1995 年の労働力調査のデータにあてはめた．
AIC による比較では，Pitman モデルが優越する結果となった．Hoshino（2003）
では Conditional（2000）が用いた 1997 年の労働力調査
Inverse Gaussian-Poisson（CIGP）モデルを，佐井・竹村
データにあてはめている．
前節の議論では，モデルとして混合ポアソン分布を用いるのが自然であった．実はこれら既存のモデルは全て，混合ポアソン分布またはその極限と関係付けられる．従って混合ポアソン
分布の性質を理解する事で，モデルの意味が一層明確になる．以下では混合ポアソン分布の重要な性質を述べる．なおここで触れられない性質については，Johnson et al.（（1993）, Chapter
等を見よ．
3）
確率変数 F が平均 λ のポアソン分布に従い，λ の密度関数は f (λ) とする．つまり我々は混合ポアソン分布，(1.4)

P(F = y) =

1
y!



∞

exp(−λ)λy f (λ)dλ ,y = 0, 1, 2, . . . ,0を考察する．Maceda（1948）によれば，分布（1.4）のモメントが存在するとして，E(F (r) ) = E(λr ) ,(1.5)

ただし F (r) = F (F − 1) · · · (F − r + 1) である．また逆に全ての次数 r について
（1.5）が成立するならば，（1.4）を満たす．言いかえれば，階乗積率母関数が混合される分布の積率母関数と等しい事が，混合ポアソン分布の必要十分条件という事である．
次に，J → ∞ という極限操作によるモデル導出を考察する．本分野では J が非常に大きくなるので，このような極限に意味がある．ただし期待値制約が有るので(1.6)

J → ∞,

ただし E(N ) = N0 固定 ,という操作を適用する．仮に Fj , j = 1, . . . , J, が独立で同一な分布に従うとしよう．Engen（1977）によれば，もしci= P(Fj = i) , i = 1, 2, . . . ,
(1.7)
lim
J →∞ J

超母集団モデルによる個票開示リスク評価

303

ならば J → ∞ の時，寸法指標 Si , i = 1, 2, . . . , の極限分布は各 i 独立に，平均 ci のポアソン分布である．つまり多項分布（1.3）の周辺の二項分布は，独立なポアソン分布となる．これはいわゆる「小数法則」に他ならない．E(N ) = N0 という制約が（1.7）を意味するのなら，極限分布として意味のあるモデルが得られる．
では極限操作（1.6）が意味を持つような分布族は存在するのだろうか．少なくとも複合ポア
ソン分布がそれにあたる．正の整数上の確率変数 X について，P(X = i) = qi , i = 1, 2, . . . , とする．X の確率母関数は

∞g(z) =

qi z i ,i=1で表される．ここで確率母関数(1.8)

G(z) = exp(a(g(z) − 1)) ,a > 0,

で定義される非負整数上の分布を，複合ポアソン分布と言う．さて F1 , . . . , FJ が互いに独立に同一の分布（1.8）に従うとしよう．この時 N の分布の確率母関数は(1.9)

GN (z) = G(z)J= exp(Ja(g(z) − 1)) ,a > 0,

である．この事からも分かるように，複合ポアソン分布において母数 a は平均に比例する．従って制約条件 E(N ) = N0 は，適当な定数 A について Ja = A のように書ける．Hoshino（（2002b）,
Theorem 1）によれば，ここで J → ∞ (a → 0) の時，Si は独立に平均 qi A のポアソン分布となる．なお Kemp（1978）によれば，g(z) は（1.8）のゼロ切り落とし分布の a → 0 での極限分布を定義する．この g(z) の二重性は，後述されるように混乱を生んだ．
複合ポアソンの独立同一分布モデルの場合，N の分布は Fj の分布の母数 a を Ja に置き換えれば得られる．従って N = N0 の条件付分布を導出しやすい．g(z) がべき級数分布の場合，
Hoshino（（2002b）, Theorem 2）は混合ポアソンの条件付分布 P(F1 , F2 , . . . , FJ | N ) とその（1.6）
による極限 P(S1 , S2 , . . . , SN | N ) を与えている．また J = 1, 2, . . . について Ja = A と固定される場合，N の分布は極限操作（1.6）を適用しても変化しない．この時（1.9）は，複合ポアソン分
布が無限分解可能という事を意味する．逆に Lévy の定理によれば，非負整数上の無限分解可能分布は複合ポアソン分布となる（Feller（1957）, Section 12.3 を見よ）
．すなわち N の分布が極限操作によって変化しない独立同一分布のモデルは，複合ポアソンに限る．
次に複合ポアソン分布と混合ポアソン分布の関係を考察する．Maceda（1948）は，混合ポア
ソンかつ複合ポアソンとなる例がいくらでも作れる事を指摘した．また混合ポアソン分布が無限分解可能である必要十分条件は，λ の分布が無限分解可能である事を示した．この時 Lévy の定理より，混合ポアソンは複合ポアソンである．Steutel（1983）によれば，よく使われる分布の多くが無限分解可能である．従って極限操作（1.6）と条件付け (N = N0 ) の組み合わせは，モデ
ル構築の方法としてかなりの一般性を持つ．
但し混合ポアソンかつ複合ポアソンとなる分布であっても，解析的に扱いやすいとは限らない．例えば対数正規分布は Thorin（1977）によれば，無限分解可能である．従って対数正規=
ポアソン分布は複合ポアソンとなるが，その解析的操作は困難である．本稿では混合ポアソン
かつ複合ポアソンであり解析的操作が容易な具体例として，λ がガンマ分布および逆ガウシア
ン分布に従う場合を紹介する．Anscombe（1950）は極限操作（1.6）を，ガンマ＝ポアソンモデルに適用し，対数級数モデルを得た．なお負の二項分布のゼロ切り落とし分布の極限が，対数級数分布になる．これは g(z) の二重性の一例であり，Fisher の記述が曖昧だったので混乱を招いた．Sibuya et al.（1964）によれば，ガンマ＝ポアソンモデルを母集団サイズで条件付けると

304

統計数理
表 1.

第 51 巻

第 2 号 2003

混合ポアソン分布モデル．

ディリクレ＝多項モデルを得る．そして Watterson（1976）, Takemura（1999）によれば，ディリ
クレ＝多項モデルに（1.6）を適用すると Ewens モデルが得られる．また対数級数モデルを母集
団サイズで条件付けると Ewens モデルが得られる（Watterson（1974）, Hoshino and Takemura（1998））．次に逆ガウシアン分布の場合だが，Hoshino（2003）は逆ガウシアン＝ポアソンモデルの条件付き分布として CIGP モデルを提案した．そして Hoshino（2002a）は逆ガウシアン＝ポ
アソンモデルに小数法則を適用し，Engen（1974）の拡張負の二項分布モデルの特殊ケースを得た．なお逆ガウシアン＝ポアソン分布のゼロ切り落とし分布の極限は，（切り落とし）拡張負の二項分布の特殊な場合である．また Hoshino（2002a）によれば，CIGP モデルに小数法則を適用すると，Limiting CIGP（LCIGP）
モデルが得られる．これは拡張負の二項分布モデルの特殊ケースの条件付分布になっている．以上の議論は表 1 のように整理される．
2.

超母集団モデル各論

本章ではこれまでに言及したモデルを紹介する．なお全てのモデルについて，定義式の引数は（定義 1.1）から（定義 1.3）を満たす．また本章では対数尤度関数を，標本寸法指標 (s1 , s2 , . . . , sn )の関数 L で表す．
2.1 ディリクレ＝多項モデル
本節の議論は Johnson et al.（1997）の 35.13.1 節，及び Takemura（1999）による．

（セル同士が交換可
定義. ディリクレ＝多項モデルでは母集団サイズ N が固定されている．
能な場合の）
モデルは γ > 0 についてP(F1 = y1 , . . . , FJ = yJ ) =N !Γ(Jγ) Γ(γ + y1 )Γ(γ + yJ )···Γ(Jγ + N ) Γ(γ)y1 !
Γ(γ)yJ !

のように表される．この時寸法指標に関しては次式で表される．

N(定義 2.1)

P(S0 = t0 , . . . , SN = tN ) =N !J!Γ(Jγ)
Γ(Jγ + N ) i=0Γ(γ + i)
Γ(γ)i!



ti

1
.
ti !

リスクの推定確率の総和が 1 になる事を利用して，階乗モメントを導出出来る．すなわち

 
J

(2.1)

E

j=1

(r )
Fj j

=



J
[rj ]
j=1 γ
,
(Jγ)[R]

N (R)



305

超母集団モデルによる個票開示リスク評価

但し R = Jj=1 rj , x(r) = x(x − 1) · · · (x − r + 1) , x[r] = x(x + 1) · · · (x + r − 1) である．寸法指標の期待値については

Ji+γ −1E(Si ) =P(Fj = i) = J
i
j=1





N − i + (J − 1)γ − 1
N −i

N + Jγ − 1.
N

母数 γ の推定については，Mosimann（1962）, Takemura（1999）等で議論されている．単純無作為非復元抽出で，大きさ n の標本が取られたとする．標本分布は（定義 2.1）において，N をn で置き換えて得られる．
最尤法から考察しよう．Levin and Reeds（1977）は，尤度関数が γ に関して単峰と証明した．
最尤推定量は dL/dγ = 0 の解である．

      
n

i−1

dL
=
si
dγ
i=1
n

d2 L
=
si
dγ 2
i=1

1
γ
+
j
j=0

n−1

−

J
,
Jγ
+j
j=0

i−1

−1
(γ
+ j)2
j=0

n−1

+

J2
(Jγ + j)2
j=0

より，ニュートン＝ラフソン法など高速に収束する数値解法が利用できる．
次にモメント法を試す．
1
T =i(i − 1)sin(n − 1) i=2の時，Takemura（1999）はγ̂T akemura =を推定量としている．また標本分散を
(2.2)

2

s =



1−T
JT − 1i=0 si (i − n/J)

2

J −1

と書く．この時近似的に Bethlehem et al.（1990）の推定量nγ̂Bethlehem =(2.3)
J(Js2 /n − 1)を得る．ただし元々これは，ガンマ＝ポアソン分布の推定量として提案された．
2.2 ガンマ＝ポアソンモデル
Greenwood and Yule（1920）によれば，（1.4）で λ がガンマ分布に従う場合，(2.4)

P(F = y) =Γ(y + γ) γ y
p ·q ,
Γ(γ)y!

y = 0, 1, 2, . . . ,但し 0 < γ, 0 < p < 1, q = 1 − p であり，これは負の二項分布である．負の二項分布の性質については Johnson et al.（1993）の Chapter 5 を見よ．
定義. ガンマ＝ポアソンモデルでは，母集団サイズが確率変数である．すなわち正の γ, β について，JΓ(yj + γ) γ yj
p ·q ,
P(F1 = y1 , . . . , FJ = yJ ) =Γ(γ)yj !
j=1







ただし p = 1/(N0 β + 1), q = 1 − p のように表される．なお寸法指標で表せば∞(定義 2.2)

P(S0 = t0 , S1 = t1 , . . . ) = J!

i=0

Γ(i + γ) γ i
p ·q
Γ(γ)i!

ti

1
ti !

306

第 51 巻統計数理

第 2 号 2003となる．なお期待値制約 E(N ) = N0 を満たすため，(2.5)

γβ = 1/J .

リスクの推定周辺の負の二項分布（2.4）については，階乗モメント
E(F



(r)

)=



が評価できる．寸法指標の期待値は
J

(2.6)

E(Si ) =P(Fj = i) = J

j=1

r

q
p

Γ(i + γ)Γ(γ)i!

γ [r]1N0 β + 1

 
γ

·

N0 β
N0 β + 1

i.

負の二項分布に関する母数の推定については Johnson et al.（1993）の 5.8 節，Engen（1978）の 3.4 節等が詳しい．標本分布については母集団サイズがランダムなので，抽出率 n0 /N0 のベ
ルヌーイ抽出を前提とする．この時（定義 2.2）について N を n で置き換えれば，標本分布が得られる．
ここで最尤法を考えよう．制約式（2.5）の下では


∞

∂L
=
si
∂γ
i=0

1
1
+ ··· +− log
γ
γ +i−1(γ + i)n0n0+1 +− 1 − log γJγγ(n0 + Jγ)=0を数値的に解けば良い．なお
∞

∂2L
=
si
∂γ 2
i=0

−1
−1
n0
+··· ++
γ2(γ + i − 1)2γ(n0 + Jγ)
n0 γ(n0 + Jγ) − (n0 + 2Jγ)(n0 γ + in0 )1+−
γ 2 (n0 + Jγ)2γである．なおモメント推定量（2.3）は反復数値解法において，初期値として用いる事も出来るだ
ろう．
2.3 対数正規＝ポアソンモデル
Preston（1948）のアイデアに基づき，（1.4）において λ が対数正規分布に従うとする．この時(2.7)

P(F = y) =1√
y! 2πV∞0

λy−1 exp(−λ − (log λ − M )2 /2V )dλ,y = 0, 1, 2, . . . ,のようになる．これを対数正規＝ポアソン分布と呼ぼう．なお本節の議論は，Shaban（1988）に大部分含まれている．
定義.







J

P(F1 = y1 , . . . , FJ = yJ ) =

1
√
y
!
2πV
j
j=1

0

∞

λyj −1 exp −λ −ただし V > 0, ∞ > M > −∞ である．寸法指標について表せば(定義 2.3)

対数正規＝ポアソンモデルにおいて，母集団サイズは確率変数である．すなわちP(S0 = t0 , S1 = t1 , . . . )(log λ − M )2
2V

dλ ,







超母集団モデルによる個票開示リスク評価
∞

= J!

i=0

1
√
i! 2πV∞0λi−1 exp −λ −(log λ − M )2
2V



307

ti

1
ti !

となる．期待値制約 E(N ) = N0 を満たすため，本稿ではM = log N0 − log J − V /2 .

リスクの推定階乗モメントは Bulmer（1974）が示している．すなわちE(F (r) ) = exp rM +

(2.8)



1 2
r V
2



である．寸法指標の期待値は（2.7）を利用して
1
E(Si ) = J √i! 2πV0∞λi−1 exp(−λ − (log λ − M )2 /2V )dλ .

母集団サイズが確率変数のモデルについては，抽出率 n0 /N0 のベルヌーイ抽出を適用するという事であった．母集団分布が（定義 2.3）で表される時，標本分布は N を n で置き換えて得られる．なおこの場合の期待値制約を
m = log n0 − log J − V /2

(2.9)

と書く．
最尤法から考察しよう．


∞

∂L
si
=
∂m
i=0



P(F = i)P(F = i)

∞


 

1
P(F = i) + 2V
P(F = i)
1
∂L
=
+
si −
∂V
2V
P(F = i)
i=0

= 0,

∞

=

si

i=0

P(F = i)P(F = i)=0を数値的に解けば最尤推定値が得られる．ここで P(F = i) の評価が必要になるが，Bulmer（1974）によれば期待値制約（2.9）が無いものとして，∂P(F = i)= iP(F = i) − (i + 1)P(F = i + 1) ,∂m
1
∂P(F = i)= {i2 P(F = i) − (i + 1)(2i + 1)P(F = i + 1) + (i + 1)(i + 2)P(F = i + 2)}∂V2となる．いずれにしても数値積分は避けられない．なお期待値制約（2.9）の下ではラグランジュ
乗数法を用いれば，数値解が得られるだろう．
モメント推定量としては，V̂M oment = log(s2 − n/J) − 2 log n + 2 log Jを得る事が出来る．ただし s2 は標本分散（2.2）である．
2.4 逆ガウシアン＝ポアソンモデル
逆ガウシアン（Inverse Gaussian）
分布については，Seshadri（1999）のモノグラフが応用も含めて網羅的である．他に Seshadri（1993），Johnson et al.（1994）も参考になる．
（1.4）において λ が IG 分布に従う時，Holla（1966）によれば(2.10)

P(F = y) =√(αθ/2)y2αexp(α 1 − θ)
Ky−1/2 (α) ,πy!

y = 0, 1, 2, . . . ,308

第 51 巻統計数理



第 2 号 2003



ただし 0 < θ ≤ 1, α > 0 となる．ここで Ky−1/2 (·) は y − 1/2 次第三種変形ベッセル関数であり，Ky−1/2 (α) =πexp(−α)2αy−1(y − 1 + i)!
(2α)−i
(y
−
1
−
i)!i!
i=0

,

y = 1, 2, . . . ,また K−1/2 (α) = K1/2 (α) となる．本稿では（2.10）を IG ＝ポアソン分布と呼ぶ．本分布は
Seshadri（（1999）, Section 7.1）において，研究史も含めて詳しく解説されている．また Johnson（1993）, Section 11.15）も参考になる．
et al.（
定義. 逆ガウシアン＝ポアソンモデルでは母集団サイズが確率変数であり，以下のように
書ける．


J

P(F1 = y1 , . . . , FJ = yJ ) =j=1√(αθ/2)yj2αexp(α 1 − θ)Kyj −1/2 (α) ,πyj !

ただし 0 < θ ≤ 1, α > 0 である．寸法指標については(定義 2.4)

 

P(S0 = t0 , S1 = t1 , . . . )
∞

= J!

i=0

√
(αθ/2)i2αexp(α 1 − θ)Ki−1/2 (α)
π
i!



ti

1
ti !

となる．期待値制約 E(N ) = N0 は，次式と同値である．
√
2N0 1 − θ.
α=Jθ
リスクの推定IG ＝ポアソン分布の r 次階乗モメントは，IG 分布の r 次モメントである．すなわち√
rKr−1/2 (α 1 − θ)αθ√√E(F (r) ) =.
2 1−θ
K−1/2 (α 1 − θ)







寸法指標の期待値は次のように書ける．
E(Si ) = J√(αθ/2)i2αexp(α 1 − θ)Ki−1/2 (α) .
πi!

次に標本分布を考察する．母集団サイズが確率変数なので，抽出率 n0 /N0 のベルヌーイ抽出を考えよう．この時母集団分布が（定義 2.4）ならば，標本分布は N を n で置き換えて得られる．
最尤法から考察する．期待値制約が無いものとして，

∞√Ki−1/2 (α)∂L
2n + Jsi=+J 1−θ+,∂α2α
Ki−1/2 (α)
i=0
Jα
n
∂L
=− √
+ .
∂θ
θ
2 1−θ

実は良く知られているように，γKγ (α)
= −Rγ (α) + ,Kγ (α)
α
ただし Rγ (α) = Kγ+1 (α)/Kγ (α) と書ける（例えば Seshadri（1999）, p. 125 を見よ）．これを利用して，尤度方程式を数値的に解く事が出来る．期待値制約√
2n0 1 − θ
α=
Jθ

309

超母集団モデルによる個票開示リスク評価が有る場合はラグランジュ乗数法を用いれば良い．もしくは制約式を代入して α を消去すれ
ば，微分を必要としない最大化アルゴリズムも使える．
期待制約無しの場合，モメント推定は以下の通りになる．
1,θ̂ = 1 −2V(x)/x̄ − 1

α̂ =

2x̄



1 − θ̂

θ̂

.

Anscombe（1950）によれば，観測された頻度 0 の標本割合が寸法指標の推定の効率性に大きく影響するという．Sichel（1973）はこの議論を受けて，以下の推定量を提案した．
(2.11)

θ̂ = 1 −− log φˆ0
2x̄ + log φˆ0



2

,



1
x̄
α̂ = − log φˆ0 1 +2x̄ + log φˆ0,

ただし x̄ は標本頻度の平均，φˆ0 は観測された頻度 0 の標本割合である．
2.5 対数級数モデル
負の二項分布（2.4）から度数 0 を切り落としたP(F = y) =Γ(y + γ) γ yp q
Γ(γ)y!

(1 − pγ ),y = 1, 2, . . . ,について，γ → 0 の極限を適用するとP(F = y) =

c qy
,
y

y = 1, 2, . . . ,ただし c = −{log(1 − q)}−1 ，を得る．これを「対数級数分布」と言う．我々はこれとは区別して，Anscombe（1950）の「対数級数モデル」を考察する．
定義.

N0 > 0, 0 < p < 1, q = 1 − p について，p = 1/(N0 β + 1)，
λi = N0p · q i−1,i

i = 1, 2, . . . ,

とする．この時母集団サイズが確率変数である対数級数モデルは∞(定義 2.5)

P(S1 = t1 , S2 = t2 , . . . ) =λtii exp(−λi )
ti !
i=1

のように定義される．この時 E(N ) = N0 となっている．
リスクの推定対数級数モデルでは，Si は平均 λi のポアソン分布に従う．従って(r)

E(Si ) = λi rと書ける．Johnson et al.（1993）の 4 章を見よ．なお各 i 独立であり，特に E(Si ) = λi .
母数の推定は Hoshino and Takemura（1998）で議論された．抽出率 n0 /N0 のベルヌーイ抽出を考えよう．1.2 節で確認したように，ポアソン分布からのベルヌーイ標本はポアソン分布に従う．従って対数級数モデルの標本分布も対数級数モデルであり，母集団分布の N0 を n0 に置
き換えればよい．なお応用上は n = n0 とみなす事に注意．
対数尤度 L を β で微分して，最尤方程式を得る．すなわち最尤推定量 β̂ は，−uβ + log(n0 β + 1) = 0

310

統計数理

第 51 巻

第 2 号 2003の解を数値的に評価すれば良い．上式の左辺をもう一度微分してn0−u +n0 β + 1が得られるので，ニュートン＝ラフソン法等が使える．
モメント推定量は，例えば
2s2β̂ =n(s1 − 2s2 )等が考えられる．
2.6 Ewens モデル
Ewens（1972）が導入した分布を Ewens 分布，または Ewens Sampling Formula と呼ぶ．本分布については Johnson et al.（1997）の 41 章で，詳しく解説されている．

定義.

母数 θ > 0 について，母集団サイズ N が固定されたモデル(定義 2.6)

P(S1 = t1 , S2 = t2 , . . . , SN = tN ) =

θU
θ[N]



N!
N
ti
i=1 i ti !

,ただし θ[N] = θ(θ + 1)(θ + 2) · · · (θ + N − 1), を Ewens 分布モデルとする．
リスクの推定寸法指標の同時階乗モメントは Sibuya（1993）が与えている．すなわち

 

N

E

Si (ri )i=1

ただし r = r1 + · · · + rN , s =

 
N

=

θr θ[N−s] N (s)
θ[N]
i=1

1
i!

ri

,

N
i=1 iri ≤ N ．特にi

E(Si ) =θN −j+1.
i j=1 θ + N − j

Ewens 分布は個体のラベルに依存しない．故に大きさ N の母集団から n 個の標本を単純無作為非復元抽出する場合，標本分布は（定義 2.6）について N, U を n, u で置き換えて得られる．
母数 θ の最尤推定量は渋谷（1991）が示したように，


n

u
1
−
=0
θ j=1 θ − 1 + j



の解である．二次の微分係数

n

u
∂2L
1
=− 2 +
∂θ2
θ
(θ
−
1
+ j)2
j=1

を利用してニュートン＝ラフソン法等を用いれば良い．なお渋谷（1991）が指摘するように，
Ewens 分布は指数分布族に所属する．従って最尤推定は一意な解を持つ．
簡単なモメント推定量を挙げよう．E(s1 ) = θn/(θ + n − 1) より
θ̂ =を導く事が出来る．

s1 (n − 1)n − s1

311

超母集団モデルによる個票開示リスク評価

2.7 Pitman モデル
Pitman（1995）
は Ewens 分布を二母数へ拡張し，Pitman 分布を得た．これは Pitman Sampling
Formula とも呼ばれる．Ewens 分布と同様に寸法指標の同時確率を与える為，モデルとして使用できる．

定義. 母数 0 ≤ α < 1, θ > −α またはある自然数 m について α < 0, θ = −mα を満たすような組み合わせについて，Pitman モデルは以下の様に定義される．すなわち固定された母集団
サイズ N について，

N(定義 2.7)

P(S1 = t1 , S2 = t2 , . . . , SN = tN ) = N !

θ[U :α]θ[N] j=1(1 − α)[j−1]
j!



tj

1
,
tj !

ただし θ[U :α] = θ(θ + α) · · · (θ + (U − 1)α), θ[N] = θ(θ + 1) · · · (θ + N − 1).
もし α が 0 ならば，
（定義 2.7）は Ewens モデル（定義 2.6）と一致する．また α < 0 の場合 θ =−Jα > 0, γ = −α > 0 とおけば，（定義 2.7）はディリクレ＝多項モデル（定義 2.1）となる．
リスクの推定Yamato and Sibuya（2000）が，寸法指標の同時階乗モメントを評価している．すなわち

 

N

E

Si (ri )Nθ[r:α] (θ + rα)[N−s] N (s)
θ[N]
i=1

=

i=1

ただし r = r1 + · · · + rN , s =N
i=1 iri ≤ N ．特にE(Si ) =(1 − α)[i−1] N (i)θi!

(1 − α)[i−1]i!

(θ + α)[N−i]
θ[N]





ri

,

.

Pitman 分布は個体のラベルに依存しない．故に大きさ N の母集団から n 個の標本を単純無作為非復元抽出した場合，標本分布は母集団分布（定義 2.7）について N, U を n, u で置き換えて得られる．
最尤推定量は，以下の同時方程式の解である．



u−1



n−1

∂L
1
1
=
−
= 0,
∂θ
θ
+
iα
θ
+
i
i=1
i=1



u−1

n

i−1

∂L
i
1
=
−
= 0.
si
∂α
θ
+
iα
J
−
α
i=1
i=2
j=1






これは数値的に評価する必要が有る．例えば Hoshino（2001）は二次の微分係数，
u−1

n−1

∂2L
1
1
=−
+
,
2
(∂θ)2
(θ
+
iα)
(θ
+
i)2
i=1
i=1



u−1

n

i−1

∂2L
i2
1
=−
−
si
< 0,
2
2
(∂α)
(θ
+
iα)
(j
−
α)2
i=1
i=2
j=1



u−1

i
∂2L
<0
=−
∂θ∂α
(iα
+
θ)2
i=1

を利用し，ニュートン＝ラフソン法を用いた．

312

統計数理

第 51 巻

第 2 号 2003

次に近似的モメント推定量を紹介する．
θ̂ =nuc − s1 (n − 1)(2u + c),
2s1 u + s1 c − ncα̂ =θ̂(s1 − n) + (n − 1)s1,
nu

ただし c = s1 (s1 − 1)/s2 ．導出は Hoshino（2001）の Appendix を見よ．
2.8 条件付逆ガウシアン＝ポアソン（CIGP）モデル
Hoshino（2003）は逆ガウシアン＝ポアソンモデル（定義 2.4）について，母集団サイズ N を固定した条件付分布を CIGP 分布と呼び，その性質を議論した．これをモデルとして用いる．

定義. 母数 α > 0 について，母集団サイズ N 所与の CIGP モデルは以下の様に定義される．
(定義 2.8)

 

P(S0 = t0 , . . . , SN = tN )
2α
π

=

J−1
2

N

J!N !
J N+1/2 KN−1/2 (Jα) i=0

Ki−1/2 (α)
i!

リスクの推定寸法指標の同時階乗モメントは Hoshino（2003）によれば，   



N

E

(r )
Sj j

j=1

ただし r ==2απN
j=1 rj , R =E(Si ) =

r
2


N

N !J!KN−R−1/2 ((J − r)α)(J − r)N−R+1/2(N − R)!(J − r)!J N+1/2 KN−1/2 (Jα) j=1



ti

1
.
ti !

Kj−1/2 (α)
j!



rj

,

(rj )
N
= Sj (Sj − 1) · · · (Sj − rj + 1)，となる．特に
j=1 jrj , Sj

2α Ki−1/2 (α) N !KN−i−1/2 ((J − 1)α)(J − 1)N−i+1/2.
πi!
(N − i)!J N−1/2 KN−1/2 (Jα)

CIGP 分布は個体のラベルに依存しない．故に（定義 2.8）に従う大きさ N の母集団から n 個の標本を単純無作為非復元抽出した場合，標本分布は（定義 2.8）の N を n に置き換えたものとなる．
最尤推定量は尤度方程式 dL/dα = 0 の解である．ここで


n



dL
= JRn−1/2 (Jα) −si Ri−1/2 (α) ,
dα
i=0

 

ただし Rγ (α) = Kγ+1 (α)/Kγ (α) である．また二次の微分係数は
d2 L
= J2
dα2

2
Rn−1/2(Jα) +

2n
Rn−1/2 (Jα)
Jα2si Ri−1/2
(α) +

−

i=0

2i
Ri−1/2 (α)α

− J2 + Jで与えられる．
モメント推定については，変形ベッセル関数の評価が問題になる．すなわち厳密な推定量は計算上不便である．従って Hoshino（2003）は，IGP 分布のモメント推定量を利用して近似的モ
メント推定量を提案した．標本分散（2.2）を利用すれば，
α̃ =



n

n(2Js2 − n)J(Js2 − n)

が一つの近似的推定量である．また s0 を用いた IGP 分布の推定量（2.11）に対応して
1n/Jᾱ = − (log s0 − log J) 1 +2
n/J + log s0 − log J313

超母集団モデルによる個票開示リスク評価が提案された．Hoshino（2003）による実データへの当てはめでは，後者の近似的推定値が前者
より最尤推定値に近い傾向が見られた．
2.9 拡張負の二項分布モデル
本節の議論は Hoshino（2002b）に詳しい．ガンマ＝ポアソンモデルでの寸法指標 Si の期待値
（2.6）は，期待値制約（2.5）を代入すると(2.12)

E(Si ) =1 Γ(i + γ)β Γ(1 + γ)i!



1
N0 β + 1


γ

N0 β
N0 β + 1iと書ける．対数級数モデルでは γ → 0 の挙動を考察したが，Engen（1974）
は（2.12）
式の Γ(1 + γ)に着目し，i ≥ 1 の時 −1 < γ < 0 と出来ると指摘した．このようなモデルを，「拡張負の二項分
布（ENB）モデル」と呼ぶ．しかし寸法指標の同時分布は与えられず，Engen のアイデアは「拡張（切り落とし）負の二項分布」として利用されてきた．確率関数で表すとP(F = i) =−γΓ(i + γ) iθ ∝ E(Si ) ,1 − (1 − θ)−γ Γ(1 + γ)i!

i = 1, 2, . . . ,ただし −1 < γ < 0, 0 < θ ≤ 1 である．これは ENB モデルとは区別しなければならない．より
詳しくは Engen（1978），または Johnson et al.（1993）の 5.12.2 節を見よ．
定義. −1 < γ (γ = 0), 0 < θ < 1 について，拡張負の二項分布モデルは次のように定義さ
れる．


∞

(定義 2.9)

P(S1 = t1 , S2 = t2 , . . . ) =exp(−τ (i; γ, θ))τ (i; γ, θ)ti
,
ti !
i=1

ただし

N0 (1 − θ)γ+1 θi Γ(i + γ)θΓ(1 + γ)i!
である．このモデルを ENB(γ) と表す．期待値制約 E(N ) = N0 を満たす．
τ (i; γ, θ) =θ = 1 の時，逆ガウシアン＝ポアソンモデルの極限として得られるのが ENB(−1/2) である．
また逆ガウシアン＝ポアソン分布のゼロ切り落とし分布で α → 0 とすると，拡張（切り落とし）負の二項分布 (γ = −1/2) を得る．なお Hoshino（2002b）によれば，一般化された逆ガウシア
ン＝ポアソンモデルの極限として一般の ENB(γ) が得られる．

リスクの推定Si がポアソン分布に従う事より，(r)E(Si ) = τ (i; γ, θ)rのように書ける．なお各 i 独立であり，特に E(Si ) = τ (i; γ, θ).
母数推定の議論について，詳しくは Hoshino（2002b）を見よ．拡張負の二項分布モデルでは，
ベルヌーイ抽出が用いられる．抽出率が n0 /N0 の場合，標本分布は（定義 2.9）について N を nで置き換えたものになる．
対数尤度について，一次の微分係数は以下の式で与えられる．
n0
∂L
=−
∂θ
γ



(1 + γ)(1 − θ)γ(1 − θ)γ+1 − 1
+
2
θ
θ

∂L
n0
1−θ
= n0
+
∂γ
θγ 2
θ



−



+ (n − u)(γ + 1)
1
−u
,
θ
1−θ

(1 − θ)γ+1(1 − θ)γ+1log(1 − θ)
+
γ2
γ



314

統計数理

第 51 巻i−1∞+ u log(1 − θ) + S1 +

si

i=2

第 2 号 2003

1
.
γ
+
j
j=1

また二次の微分係数は以下の式で与えられる．
n0 ((1 − θ)γ+1 − 1)n0∂2L=
− 2 (1 − θ)γ+1 log(1 − θ)∂γ∂θ
θ2 γ 2
θ γ
n0
n0
u
n0
(1 − θ)γ log(1 − θ) −(1 − θ)γ log(1 − θ) −.
+ 2 (1 − θ)γ −
θγ
θγ
θ
1−θ
∂2L
n0 2 − 2(1 − θ)γ+12(1 + γ)(1 − θ)γ(1 + γ)γ(1 − θ)γ−1
=−
−
−
2
3
2
∂θ
γ
θ
θ
θ
n−u
u(γ + 1)−−.
θ2
(1 − θ)2






∂2L
n0
=
∂γ 2
θ

2(1 − θ)γ+1(1 − θ)γ+1−2(1 − θ) + 2(1 − θ)γ+1
2
−
log(1
−
θ)
+
log
(1
−
θ)
γ3
γ2
γ

∞

−





i−1

si

i=2

1
.
(j
+
γ)2
j=1

最尤推定は，同時方程式 ∂L/∂θ = ∂L/∂γ = 0 を解いて得られる．ただし応用の際は n0 を nの実現値で置きかえるので，n = n0 となる．この時尤度方程式は簡略化されて，以下の様になる．
n0 (1 − (1 − θ)γ )(1 − θ)
∂L
=0⇔u=
,
∂θ
γθ
∂L
n0 (1 − θ)
=0⇔
∂γ
γθ



 
∞

i−1

1 − (1 − θ)γ1si+ log(1 − θ) +
= 0.
γ
γ
+
j
i=2
j=1

γ が所与でない限り，モメント推定であっても数値的に方程式を解く事になる．簡便な近似
推定の余地は有ると思われるが，議論はなされていない．なお γ が −1 に近い程，小さなセルが支配的となる．
2.10 極限 CIGP モデル
Hoshino（2002a）
は CIGP 分布（定義 2.8）
に極限操作
（1.6）
を適用した．ここで得られた分布をLimiting CIGP（LCIGP）
分布と呼び，モデルとして用いる事が出来る．

定義.

母数 A > 0 について，LCIGP モデルは次のように表される．

(定義 2.10)

P(S1 = t1 , . . . , SN = tN )
=


N

π
N ! exp(−A)
2 AN−U +1/2 KN−1/2 (A) i=1(2i − 3)!!
i!



ti

1
.
ti !

なお本モデルでは母集団サイズが固定されている．
リスクの推定Hoshino（2002a) によれば，寸法指標の同時階乗モメントは次式で与えられる．

 


N

i=1

ただし r =

(r )

Si i

E

N
i=1 ri , R =


N

=

KN−R−1/2 (A)Ar−R N !
KN−1/2 (A)(N − R)! i=1(2i − 3)!!
i!



ri

N
(R)
= n(n − 1) · · · (n − R + 1) である．
i=1 iri , n,

超母集団モデルによる個票開示リスク評価

315

特に i = 1, 2, . . . , N についてE(Si ) =KN−i−1/2 (A)A1−i (2i − 3)!!N !
.
KN−1/2 (A)i!(N − i)!

単純無作為非復元抽出の場合，大きさ n の標本分布は母集団分布（定義 2.10）において N, U
を n, u で置き換えて得られる．
対数尤度について一次の微分係数は次のとおり．
∂L1= −1 − (2n − u) + Rn−1/2 (A) ,∂A
A

ただし Rγ (ω) = Kγ+1 (ω)/Kγ (ω) である．最尤推定量は，∂L/∂A = 0 の解である．これを解くには二次の微分係数
2n
1
∂2L
2
Rn−1/2 (A) + (2n − u) 2 − 1= Rn−1/2
(A) −
∂A2
A
A

を用いて，ニュートン＝ラフソン法が利用できる．LCIGP 分布は指数族に所属するので，最尤推定の挙動は良好である．
LCIGP 分布について厳密なモメント推定量は，変形ベッセル関数の評価を伴い実用的ではない．従って Hoshino（2002a）は，以下の近似的モメント推定量を提案した．
Ã =

3.



u/(1 −u,

1 − 4s2 /s1 ) ,もし 4s2 /s1 ≤ 1 ,その他の場合 .

展望

経験的に s1 − s2 は，si − si+1 , i = 2, 3, . . . , に比べて変化率が大きい．つまり寸法指標は滑らかに変化しないので，s1 の「跳び」を記述出来る分布が望ましい．そして「跳び」という異質性を記述するには，追加の母数が必要と思われる．従って，二母数のモデルが有利だろう．2
章で紹介されたモデルのうち，二母数モデルは Pitman モデルと ENB モデルのみである．故に Pitman モデルを ENB モデルが改善するかどうかは興味深い．ところが Hoshino（2002b）によれば，二つのモデルは N, U 所与の条件付分布が同じになる．この場合，あてはまりは同程度だろう．結局，これらの分布とは別の強みを持つモデルが望まれる．
モデルの開発に試行錯誤は付き物だが，1.3 節で述べた Hoshino（2002b）の一般論を利用すれば，錯誤は減るだろう．基本的には正の整数上の分布 P(i), i = 1, 2, . . . , が，si /u を記述出来れば良い．この時，E(Si ) ∝ P(i), i = 1, 2, . . . , となる独立ポアソン分布のモデルが，小数法則の極限として正当化される．そのようなモデルで N の分布は複合ポアソンなので，N = N0 という条件付けの際に組み合わせ論的困難は存在しない．
このような正の整数上の分布によるモデル構築を考えると，我々は Zipf 以来の古典的議論に回帰した事に気付く．すなわち，計量言語学や統計的生態学においてゼロ切り落とし分布をあてはめる議論である．これら単純なあてはめは発展性が無いため，その議論は忘れられかけている．しかし本稿で示した観点からは，この文脈で用いられた正の整数上の分布が，新たな意
味を持つ．例えばべき級数分布族の中でのあてはまりの比較は，Wani and Lo（1986）によって議論された．これは望ましい新規モデルについての指針を与える．
形式的に本分野の超母集団モデルは離散多変数分布である．Johnson et al.（1997）の序文で“less has been done” と述べられているように，実用的な離散多変数分布はそれ程多く知られていない．従って超母集団モデルの開発は，確率論的にも価値がある．例えば P(F1 , F2 , . . . , FJ | N )の下での U の挙動は，壺モデルでの占有問題に対応する．なお U は計量言語学では語彙数，統

316

統計数理

第 51 巻

第 2 号 2003

計的生態学では種の総数を表す重要な統計量である（Bunge and Fitzpatric（1993）の総合報告を見よ）．また P(S1 , . . . , SN | N ) は，自然数の確率分割と同等である．つまり小数法則と条件付けを組み合わせる事で，体系的にこれらの離散多変数分布が得られる．
本分野における超母集団モデルの議論が，著者にとっては思いもよらなかった一般性を見せている事を，ここまでで報告したつもりである．本稿がこの分野に興味を持ってもらうきっか
けになる事を期待しつつ，とりあえず筆をおくことにする．
謝

辞

本研究は統計数理研究所の共同利用研究プロジェクト「個票データの開示におけるリスクの評価と官庁統計データの公開への応用（代表：佐井至道）」の成果に基づくものであり，科学研
究費の援助も受けている．

参 考 文 献
Anscombe, F. J.（1950）. Sampling theory of the negative binomial and logarithmic series distributions, Biometrika, 37, 358–382.
Bethlehem, J. G., Keller, W. J. and Pannekoek, J.（1990）. Disclosure control of microdata, J. Amer.
Statist. Assoc., 85, 38–45.
Bulmer, M. G.（1974）. On ﬁtting the Poisson lognormal distribution to species-abundance data,
Biometrics, 30, 101–110.
. Estimating the number of species: A review, J. Amer. Statist.
Bunge, J. and Fitzpatrick, M.（1993）
Assoc., 88, 364–373.
Chen, G. and Keller-McNulty, S.（1998）. Estimation of identiﬁcation disclosure risk in microdata,
Journal of Oﬃcial Statistics, 14, 79–95.
. On species frequency models, Biometrika, 61, 263–270.
Engen, S.（1974）
Engen, S.（1977）. Comments on two diﬀerent approaches to the analysis of species frequency data,
Biometrics, 33, 205–213.
. Stochastic Abundance Models, Chapman and Hall, London.
Engen, S.（1978）
Ewens, W. J.（1972）. The sampling theory of selectively neutral alleles, Theoretical Population Biology, 3, 87–112.
. An Introduction to Probability Theory and Its Applications, Vol. 1, 2nd ed., Wiley,
Feller, W.（1957）
New York.
Feller, W.（1966）. An Introduction to Probability Theory and Its Applications, Vol. 2, Wiley, NewYork.
Fisher, R. A., Corbet, A. S. and Williams, C. B.（1943）. The relation between the number of species
and the number of individuals in a random sample of an animal population, Journal of Animal
Ecology, 12, 42–58.
. The population frequencies of species and the estimation of population parameters,
Good, I. J.（1953）
Biometrika, 40, 237–264.
Greenberg, B. V. and Zayatz, L. V.（1992）. Strategies for measuring risk in public use microdata ﬁle,Statist. Neerlandica, 46, 33–48.
. An inquiry into the nature of frequency distribution repreGreenwood, M. and Yule, G. U.（1920）
sentative of multiple happenings with particular reference to the occurrence of multiple attacks
of disease of or repeated accidents, J. Roy. Statist. Soc., 83, 255–279.
Holla, M. S.（1966）. On a Poisson-inverse Gaussian distribution, Metrika, 11, 115–121.

超母集団モデルによる個票開示リスク評価

317

Hoshino, N.（2001）. Applying Pitman’s sampling formula to microdata disclosure risk assessment,
Journal of Oﬃcial Statistics, 17, 499–520.
Hoshino, N.（2002a）. On limiting random partition structure derived from the conditional inverse
Gaussian-Poisson distribution, Tech. Report, CMU-CALD-02-100, School of Computer Science, Carnegie Mellon University, Pittsburgh.
Hoshino, N.（2002b）. Engen’s extended negative binomial model revisited, Discussion Paper, No.
2002-1, Faculty of Economics, Kanazawa University, Kanazawa.
. Random clustering based on the conditional inverse Gaussian-Poisson distribuHoshino, N.（2003）
tion, J. Japan Statist. Soc., 33, 105–117.
. Relationship between logarithmic series model and otherHoshino, N. and Takemura, A.（1998）
superpopulation models useful for microdata disclosure risk assessment, J. Japan Statist. Soc.,
28（2）, 125–134.
. Univariate Discrete Distributions, 2nd ed., Wiley,Johnson, N. L., Kotz, S. and Kemp, A. W.（1993）
New York.
. Continuous Multivariate Distributions, Vol.
Johnson, N. L., Kotz, S. and Balakrishnan, N.（1994）
1, 2nd ed., Wiley, New York.
Johnson, N. L., Kotz, S. and Balakrishnan, N.（1997）. Discrete Multivariate Distributions, Wiley,
New York.
. Cluster size probabilities for generalized Poisson distributions, Comm. Statist.
Kemp, A. W.（1978）
Theory Methods, 7, 1433–1438.
Lehmann, E. L. and Casella, G.（1998）. Theory of Point Estimation, 2nd ed., Springer, New York.
Levin, B. and Reeds, J.（1977）. Compound multinomial likelihood functions are unimodal: Proof ofa conjecture of I. J. Good, Ann. Statist., 5, 79–87.
. On the compound and generalized Poisson distributions, Ann. Math. Statist.,
Maceda, E. C.（1948）
19, 414–416.
Mandelbrot, B. B.（1983）. The Fractal Geometry of Nature, W. H. Freeman and Company, NewYork.
Marsh, C., Skinner, C., Arber, S., Penhale, B., Openshaw, S., Hobcraft, J., Lievesley, D. and Walford,
N.（1991）. The case for a sample of anonymised records from the 1991 census, J. Roy. Statist.
Soc. Ser. A, 154, 305–340.
. On the compound multinomial distribution, the multivariate β-distribution
Mosimann, J. E.（1962）
and correlations among proportions, Biometrika, 49, 65–82.
Omori, Y.（1999）. Measuring identiﬁcation disclosure risk for categorical microdata by posterior
population uniqueness, Statistical Data Protection—Proceedings of the Conference, Lisbon,
25 to 27 March 1998-1999 Edition, 59–76, Oﬃce for Oﬃcial Publications of the European
Communities, Luxembourg.
. Cours d’Économie Politique, F. Rouge, Lausanne.
Pareto, V.（1897）
Pitman, J.（1995）. Exchangeable and partially exchangeable random partitions, Probab. Theory Related Fields, 102, 145–158.
. The commonness, and rarity, of species, Ecology, 29, 254–283.
Preston, F. W.（1948）
佐井至道（1998）. 個票データにおける個体数とセル数との関係，応用統計学，27, 127–145.
佐井至道（2000）. 予測個体数の期待値に基づく個票データのリスク評価，統計数理，48, 229-251
佐井至道，竹村彰通（2000）. 個票データにおける分類の併合モデル，応用統計学，29, 63–82.
Samuels, S. M.（1998）. A Bayesian, species-sampling-inspired approach to the uniques problem in
microdata disclosure risk assessment, Journal of Oﬃcial Statistics, 14, 373–383.
Seshadri, V.（1993）. The Inverse Gaussian Distribution, Clarendon Press, Oxford.
Seshadri, V.（1999）. The Inverse Gaussian Distribution, Lecture Notes in Statist., 137, Springer,

318

統計数理

第 51 巻

第 2 号 2003

New York.
. Poisson-lognormal distributions, Lognormal Distributions: Theory and AppliShaban, S. A.（1988）
cations (eds. E. L. Crow and K. Shimizu), 195–210, Marcel Dekker, New York.
Shlosser, A.（1981）. On estimation of the size of the dictionary of a long text on the basis of a sample,
Engineering Cybernetics, 19, 97–102.

渋谷政昭（1991）. あるクラスタ数分布と，その同音語の解析への応用，応用統計学, 20, 139–153.
. A random clustering process, Ann. Inst. Statist. Math., 45, 459–465.
Sibuya, M.（1993）渋谷政昭（2000）. 調査データ公有化における理論的技術的課題, 統計調査制度とミクロ統計の開示（松田芳郎，濱砂敬郎，森 博美 編著）, 145–167，日本評論社，東京．
渋谷政昭（2002）. 母集団と標本で孤立している個体数，文部科学省科学研究費補助金研究成果報告書「調査データの公有化における理論的問題」（課題番号 11480055）, 25–34.
Sibuya, M., Yoshimura, M. and Shimizu, R.（1964）. Negative multinomial distribution, Ann. Inst.
Statist. Math., 16, 409–426.
. The density and size distribution of diamonds, Bull. Int. Statist. Inst., 45,
Sichel, H. S.（1973）
420–427.
Skinner, C. J. and Holmes, D. J.（1993）. Modelling population uniqueness, Proceedings of the International Seminar on Statistical Conﬁdentiality, Dublin, 175–199.
Skinner, C. J. and Holmes, D. J.（1998）. Estimating the re-identiﬁcation risk per record in microdata,
Journal of Oﬃcial Statistics, 14, 361–372.
Skinner, C. J., Marsh, C., Openshaw, S. and Wymer, C.（1994）. Disclosure control for census microdata, Journal of Oﬃcial Statistics, 10, 31–51.
Smith, T. M. F.（1976）. The foundations of survey sampling: A review, J. Roy. Statist. Soc. Ser. A,
139, 183–195.
. Inﬁnite divisibility, Encyclopedia of Statistical Sciences, Vol. 4, 114–116, Wiley,
Steutel, F. W.（1983）
New York.
竹村彰通（1998）. 労働力調査に見られるサイズインデックス，文部科学省科学研究費補助金研究成果報
告書「統計データの個票開示における開示制限の決定理論的評価」（課題番号 09206102）, 95–104.
Takemura, A.（1999）. Some superpopulation models for estimating the number of population uniques,
Statistical data protection—Proceedings of the Conference, Lisbon, 25 to 27 March 1998–1999
Edition, 45–58, Oﬃce for Oﬃcial Publications of the European Communities, Luxembourg.
. On the inﬁnite divisibility of the lognormal distribution, Scand. Actuarial J.,
Thorin, O.（1977）
121–148.
Urzúa, C. M.（2000）. A simple and eﬃcient test for Zipf’s law, Economics Letters, 66, 257–260.
Wani, J. K. and Lo, H. P.（1986）. Selecting a power-series distribution for goodness of ﬁt, Canad. J.
Statist., 14, 347–353.
Watterson, G. A.（1974）. Models for the logarithmic species abundance distributions, Theoretical
Population Biology, 6, 217–250.
Watterson, G. A.（1976）. The stationary distribution of the inﬁnitely-many neutral alleles diﬀusion
model, J. Appl. Probab., 13, 639–651.
. Moments of some statistics of Pitman sampling formula, Bull.
Yamato, H. and Sibuya, M.（2000）
Inform. Cybernet., 32, 1–10.
Zayatz, L. V.（1991）. Estimation of the percent of unique population elements in a microdata ﬁleusing the sample, Statistical Research Division Report, RR-91/08, U.S. Bureau of the Census,
Washington, D.C.
Zipf, G. K.（1949）. Human Behavior and the Principle of Least Eﬀort, Addison-Wesley, Cambridge,
Massachusetts.

Proceedings of the Institute of Statistical Mathematics Vol. 51, No. 2, 297–319 (2003)

319

Microdata Disclosure Risk Evaluation with Superpopulation Models:
A Review
Nobuaki Hoshino(Faculty of Economics, Kanazawa University)

Microdata identify each record’s position in the corresponding contingency table.
Hence, the anonymization of microdata is nothing but coarsening the resolution of acontingency table. Because a finer table results in the decrement of individuals in a cell,
the frequency of cells of the same frequency of individuals plays an important role in theevaluation of disclosure risk. However, the estimation of the frequencies of frequencies ispractically impossible without an assumption on a population. Based on the standard
theory of finite population analysis, we employ superpopulation models as an assumption
for the estimation. Then Zipf’s law empirically validates the use of a Poisson distribution
mixed by a distribution with a heavy tail. A basic population model assumes that thefrequency of individuals in a cell is subject to an independent identical mixed Poisson
distribution. Let the law of small numbers indicate a limiting argument that lets thenumber of cells be infinity, where (the expectation of) the total number of individuals isfixed. A proper model arises by applying the law of small numbers to a basic model thatconsists of infinitely divisible mixed Poisson distributions. Because there are many kindsof infinitely divisible mixed Poisson distributions, it may be possible to derive useful newmodels with the law. In order to describe various populations, developing new models isof great importance.

Key words: Population unique, inﬁnite divisibility, compound Poisson, mixed Poisson, random
partitioning of the natural numbers.