Thompson, M.E., Wu, C. (2008) Simulation-based randomized systematic PPS sampling under substitution of units.　Survey Methodology, 34(1), 2-10.

　仕事の都合で読んだ奴。標本抽出デザインがややこしくて包含確率が解析的に出せないとき(ここではProportion-to-size標本抽出で回答拒否があるという場面を想定している)、計算機パワーで無理矢理なんとかしちゃうという話である。
　掲載誌はカナダ統計局が出している雑誌で、著者らはITC China Surveyというタバコについてのコホート調査の中の人らしい。

　いわく。
　complex survey [なんて訳すんですかね。多段抽出調査のことを指すと思う] では、調査ウェイトを求めるために一次包含確率[ある個体が標本に含まれる確率のこと] を求める必要があるよね。
　PPS抽出 [個体がサイズを持っていて、包含確率がサイズに比例している抽出のこと]では、一次包含確率はサイズをリスケールした値である。標本サイズ固定の非復元PPS抽出で実装が一番楽なのは無作為化系統PPSである。

　無作為化系統PPSの手順はいろいろあるが、たとえばこういうのである。
　母集団サイズを\(N\), 個体\(i\)のサイズを\(x_i\)とする。\(z_i \equiv \frac{x_i}{\sum_i^N x_i}\)とする。すべての\(i\)で\(nz_i < 1\)と仮定する。
　\(N\)個の個体をランダムに並べて、\(A_0 = 0, \ A_j = \sum_i^j nz_i\)とする。当然\(A_N = n\)となる。さて、0以上1以下の一様乱数を\(u\)として、\(k = 0, 1, \ldots, n-1\)について、\(A_{j-1} \leq u+k < A_j\)を満たす\(j\)を持つ\(n\)個を選ぶ。こうすれば、標本\(s\)への1次包含確率\( \pi_i \equiv P(i \in s)\)は無事\(\pi_i = n z_i\)となる。
　[要は、長さ\(n\)の線分を描いて、その長さを母集団メンバーのサイズに応じて分割しておき、均等に\(n\)個の点を打って、点が落ちたエリアの母集団メンバーを選ぶわけだ]

　なお、PPS抽出についての研究のほとんどは母集団合計\(T \equiv \sum_i^N y_i\)のHorvitz-Thompson推定量$$\hat{T} \equiv \sum_{i\in s} \frac{y_i}{\pi_i} $$に基づいている。\(w_i\)が\(i\)のみに依存している場合、HT推定量は\(T\)の線形推定量のなかで唯一のデザイン不偏推定量である。

　実際にはPPS抽出を実現するのは難しい。個体の回答拒否とか個体の置き換えとかのせいで、結果として\(\pi_i = n z_i\)にならないことが多い。
　例を挙げよう。国際タバコ規制(ITC)では中国の政策評価調査を行っている(ITC China Survey)。確率不均一の多段抽出デザインで、7都市から喫煙者(都市あたり800人)・非喫煙者(200人)を抽出する。都市(7)→街区(各10)→ブロック(各2)→世帯→個人という抽出が行われる。街区とブロックは無作為化系統PPS抽出、世帯と個人は基本的には単純無作為抽出である。ところが、実施上の難しさがいろいろある。街区やブロックが人手不足を理由に拒否してきたりする。仕方ないので代替するわけだけど、すると\(\pi_i\)はもうなんだかわからなくなる。拒否がランダムに起きると仮定すればなんとか求められるけど、でも…[中略]。
　そこでご提案です。デザインはわかってんだから、\(\pi_i\)はモンテカルロシミュレーションで求めましょう。

　話はかんたんである。確率抽出デザインを\(p\)とし、これに従って\(K\)個の独立な標本を抽出する。個体\(i\)が含まれた標本の数を\(M_i\)とする。\(\pi_i\)を\(\pi_i^{*} \equiv M_i/K\)で推定する。
　\(M_i\)は二項分布するから\(E(\pi_i^{*}) = \pi_i, \ Var(\pi_i^{*}) \leq 1/4K\)である。

　手法の正確さは合計のHT推定量 \(\hat{T}^{*} = \sum_{i \in s} y_i/\pi_i^{*}\)のパフォーマンスで調べることにしよう[原文では\(\tilde{T}\)だが見分けにくいので\(\hat{T}^{*}\)とメモする]。相対バイアス\(\frac{\hat{T}-\hat{T}_{MC}}{\hat{T}}\)の絶対値の下界は下式となる。任意の\(\epsilon > 0\)について$$ P \left( \frac{|\hat{T}-\hat{T}^{*}|}{\hat{T}} \leq \epsilon \right) \geq 1-\frac{2-(1+\epsilon^2)}{K \epsilon^2} \left( \sum_{i \in s} \frac{1}{\pi_i} – n \right) $$　[付録に証明が載っている。チェビシェフの不等式を使うほかは式変形だけで示せるらしいのだが、単純そうな式変形のところでつまずいてしまった… しくしく… まあとにかく、このくだりは「必要な\(K\)はどのくらいか」という話である。中略…] というわけで、\(K\)は\(10^7\)とか\(10^6\)とかあればたいてい大丈夫。
　母集団分散とか線形推定量の分散とかを求めるには二次の包含確率[2つの個体がともに標本に含まれる確率]が必要になる。たとえば\(Q = \sum_i \sum_j q(y_i, y_j)\)のHT推定量は… [さっきと似た議論の繰り返し。中略] というわけで、ここでも\(K\)は\(10^7\)とか\(10^6\)とかでたいてい大丈夫。
　残る問題は、それだけのシミュレーションが現実的な時間でできるかという話である。

　数値例。ITC China Surveyと同じデザインを考えて… [… 面倒くさいので省略]

　このように包含確率をシミュレーションで求めるやり方は、理屈のうえでは、抽出デザインが既知でありさえすればどんな抽出デザインにも適用できる。現実的には抽出デザインが複雑すぎて包含確率が解析的に求められないような場合に有用であろう。
　Fattorini (2006 Biometrika) は単位をシーケンシャルに選ぶ空間抽出でシミュレーションを使っている。PPSに代替がついた奴もシーケンシャルだとみなせるだろう[なるほど…]。なおFattoriniは本研究とは違う結果を示している。
　難関は計算時間である。無作為系統PPSはかんたんだからいいけど、同じPPS抽出でも他の方法だと時間的に厳しくなる。
　云々。

　… 仕事のなかでよく似たことをやっているので、気になってめくってみたのだが、特に驚くような話はなかった。直接の先行研究も前に読んだ Fattorini (2006) だけで、ちょっと安心した。
　あ、そうだ。シーケンシャルな抽出デザインのモンテカルロ・シミュレーションの場合、試行ごとに台帳の個体順序をランダマイズする？？するよね？？　というところを確認したかったんだった。どっちだったんだろう。この研究の場合も一種のシーケンシャルな抽出を問題にしているが、そうなる原因である回答拒否は台帳とは独立に起きるから、どっちでもいいのかもな。四の五の言わずに付録についてたRコードを読めって話ですね。