Berger, Y.G., Tille.Y. (2009) Sampling with Unequal Probabilities. Handbook of Statistics, Vol. 29A. Sample Surveys: Design, Methods and Applications. 39-54.
ポアソン抽出みたいな感じの標本抽出についてあれこれ調べていたんだけど、基礎知識が足りないのに論文ばかり読んでいても… と思って、試しに読んでみたやつ。
これもぼやきになっちゃいますけど、標本抽出論っていったいどこで学ぶんだ？ 皆目見当がつかない。
1. イントロダクション
標本抽出の研究の多くは不均一確率抽出についての研究だ。これまでに提案された抽出アルゴリズムは50を超える。解説書としてはBrewer & Hanif (1983), Tille (2006)がある。
サイズ\(N\)の有限母集団\(U\)がある。ユニットのラベルを\(k\)とする。
\(U\)の部分集合である標本\(s\)がある。抽出デザイン\(p(.)\)はすべての可能な標本上の確率測度であって、すべての\(s \in U\)について\(p(s) \geq 0\), かつ\(\sum_{s\in U} p(s) = 1\)である。\(s\)のサイズを\(n(s)\)、もしそれが定数ならば\(n\)と書く。
一次包含確率を\(\pi_k = p(k \in s)\)、同時包含確率を\(\pi_{kl} = p(k \in s \ \mathrm{and} \ l \in s)\)とする。
関心ある特性を\(y\), ユニット\(k\)の値を\(y_k\)とし、推定対象を母合計$$ Y = \sum_{k \in U} y_k$$ とする。
推定量のひとつが\(\pi\)推定量である。[Horvits-Thompson推定量のこと] $$ \hat{Y}_{\pi} = \sum_{k \in s} \frac{y_k}{\pi_k} $$ すべての\(k\)で\(\pi_k \gt 0\)ならデザイン不偏である。関心ある特性と一次包含確率の間に正の相関があるときは、均一確率抽出デザイン下よりも分散が小さくなる。相関が弱いとき向けの推定量は3節で示す。[Hajek推定のことだろう]
包含確率が既知の正のサイズ変数\(x\)と比例するという場面も多い[pps抽出のことだ]。包含確率は、すべての\(k\)で\(nx_k \leq X = \sum_{k \in U} x_k\)として $$ \pi_k = \frac{nx_k}{X} $$ となる。
不均一確率抽出デザインのもう一つの適用例は多段抽出だ。たとえば、PSUの選択確率をPSU内のSSUの数に比例させ、PSU内ではSRSにするとか。
分散推定においては\(\pi\)推定量の分散がとても大事だ。たいていの推定量は\(\pi\)推定量を含む形に線形化できるからである。
\(\hat{Y}_\pi\)の標本抽出分散は$$ var(\hat{Y}_\pi) = \sum_{k \in U} \sum_{l \in U} (\pi_{kl} – \pi_k \pi_l) \frac{y_k y_l}{\pi_k \pi_l} $$ その不偏推定量は、Horvitz-Thompsonいわく[原文では落ちているが\(var\)にハットを付けた]、$$ \hat{var}(\hat{Y}_{\pi}) = \sum_{k \in s} \sum_{l \in s} \frac{\pi_{kl} – \pi_k \pi_l}{\pi_{kl}} \frac{y_k y_l}{\pi_k \pi_l} $$ 標本サイズ固定ならSen-Yates-Grundyの推定量$$ \hat{var}(\hat{Y}_\pi) = \frac{1}{2} \sum_{k \in s} \sum_{l \in s} \frac{\pi_{kl} – \pi_k \pi_l}{\pi_{kl}} \left( \frac{y_k}{\pi_k} – \frac{y_l}{\pi_l} \right)^2 $$ も使える。すべての\(k, l \in U\)について\(\pi_{kl} > 0\)ならデザイン不偏だ。もっとも負の値をとりうるけれど。
この推定量はめったに使われない。同時包含確率は算出困難なことが多いし、総和記号がふたつあって計算が大変だからだ[いや、それはHorvitz-Thompsonのほうもそうでしょう。どちらもめったに使われないということかな]
2. 不均一確率抽出のいくつかの手法
2.1 ポアソン抽出
母集団ユニットを確率\(\pi_k\)で独立に選ぶこと。\(n(s)\)は確率変数になり、0となることもありうる。抽出デザインは$$ P(s) = \left[ \prod_{k \in s} \frac{\pi_k}{1-\pi_k} \right] \left[ \prod_{k \in U} (1-\pi_k) \right]$$ となる。
\(\pi\)推定量の分散とその不偏推定量は$$ var(\hat{Y}_\pi) = \sum_{k \in U} \frac{1}{\pi_k} (1-\pi_k) y^2_k$$ $$ \hat{var}(\hat{Y}_\pi) = \sum_{k \in s} (1-\pi_k) \frac{y^2_k}{\pi^2_k} $$ となる。
ポアソン抽出は、所与の包含確率の下で、エントロピー$$ I(p) = – \sum_{s \subset U} p(s) \log p(s)$$ を最大化するデザインである。その意味で、所与の包含確率を満たすもっともランダムな抽出だといえる。[へえー。ま、同時包含確率が\(\pi_{kl} = \pi_k \pi_l\)だとするのがいちばんランダムだということでしょうね。包含インジケータに共分散があるということはそのぶんランダムでないわけだから]
ポアソン抽出自体はめったに使われないが、無回答のモデルに使ったり、条件付きポアソン抽出の定義に使ったりする。
2.2 復元抽出
抽出確率\(p_k\)が正のサイズ変数\(x_k\)と比例しているとしよう。$$ p_k = \frac{x_k}{\sum_{l \in U} x_l} $$ 復元抽出の簡単な方法はこうだ。\(v_k = \sum_{l=1}^k p_l\)として(始点は\(v_0 = 0\)とする\)、\([0,1)\)の一様乱数\(u\)が\(v_{k-1} \leq u \lt v_k\)となるようなユニット\(k\)をとる。これを\(m\)回繰り返す。もっと効率的なアルゴリズムもある。[あるんかい]
\(Y\)のHansen-Hurwitz推定量は[中略]、これは不偏で[中略]、分散は[中略]、その不偏推定量は[中略]。
復元抽出より非復元抽出のほうが分散は減るだろう。包含確率\(\pi_k\)の非復元デザインは、もしそこでのHT推定量が\(p_k = \pi_k / n\)の非復元デザインにおけるHH推定量より常により正確ならば、良いデザインだということになる。これが満たされる条件というのがわかっていて、Rao-Sampfordデザインや固定標本サイズ最大エントロピーデザインはこれを満たす。
2.3 系統抽出
系統抽出には2種類ある。(a)母集団をなんらかの方法でソートしてから選ぶ方法。サイズ変数とかで。(b)母集団をランダム化してから選ぶ方法である。後者はランダム化系統デザインと呼ばれることもある。
系統標本はこうやって選ぶ。0と1の間の一様乱数\(u\)を生成する。$$ \pi^{(c)}_k = \sum_{j \in U; j \leq k} \pi_j $$ として[つまり包含確率の累積ね]、\(l = 1,\ldots, n\)について\(\pi^{(c)}_{k_l – 1} \lt u + l – 1 \leq \pi^{(c)}_{k_l} \)となる\(k_l\)を選ぶ。[文章で説明したほうがわかりやすくないですかね]
特に\(y\)になんらかトレンドがあるような場合、(b)より(a)のほうが正確になる。いっぽう、(a)だと標本抽出分散を不偏推定できないという欠点がある。(b)なら分散の不偏推定量がある。
2.4 Rao-Sampford抽出デザイン
まず最初のユニットを確率\(p_k = \pi_k / n\)で選ぶ。次に残りの\(n-1\)個のユニットを\(\pi_k / (\pi_k – 1)\)に比例した確率で復元抽出する。もし以上\(n\)個に重複がなかったら終了、あったら最初からやり直す。こうすると、包含確率は厳密に\(\pi_k\)になる。二次包含確率も導出できる。しかし、\(\pi_k\)が大きいとうまくいかない。改善したアルゴリズムも提案されている。
2.5 分割法による抽出
まず、\(\pi_k\)をふたつにわける。$$ \pi_k = \lambda \pi_k^{(1)} + (1-\lambda) \pi_k^{(2)} $$ \(\pi_k^{(1)}, \pi_k^{(2)} \)のどちらも0以上1以下、母集団を通じた合計を1とする。\(0 \lt \lambda \lt 1\)とする。で、確率\(\lambda\)で\(\pi_k^{(1)}\)を、確率\(1-\lambda\)で\(\pi_k^{(2)}\)を選ぶ。次に… [なんだかさっぱりわからないので断念。全然わからんが、とにかく標本サイズ固定で非復元な抽出ができるらしい]
2.6 Brewer抽出デザイン
もともとは\(n=2\)の抽出方法だったんだけど、その後任意の標本サイズに一般化された。一個づつの抽出を\(n\)ステップ繰り返すタイプの抽出方法。[めんどくさいので省略するけど、最終的に一次包含確率が狙った\(\pi_k\)になるように毎回抽出確率を変えていくのである。ふーん]
2.7 最大エントロピーないし条件付きポアソン抽出デザイン
ポアソン抽出について考えよう。包含確率を\(\tilde{\pi}_k\), 標本サイズ(確率変数)を\(\tilde{n}\)として、抽出デザインはこう書ける。$$ P(s) = \left[ \prod_{k \in s} \frac{\tilde{\pi}}{1-\tilde{\pi}} \right] \left[ \prod_{k \in U} (1-\tilde{\pi}) \right] $$ 条件付きポアソン抽出とはこうである。サイズ\(n\)であるすべての標本の集合\(\mathcal{S}_n\)に属する\(s\)について、$$ p(s) = P(s|\tilde{n}_s = n) = \frac{P(s)}{\sum_{s \in \mathcal{S}_n} P(s)} $$ そりゃそうだけどどうやって実装すんのかというと、Hajekは棄却抽出手続きを考えた。標本サイズが\(n\)になるまで何度も頑張るわけである。実際には、もっと効率的なアルゴリズムが提案されている。
難しいのは、アルゴリズム上での\(\tilde{\pi}_k\)が、デザインの包含確率\(\pi_k\)とは微妙に異なるという点である。\(\pi_k\)をどうやって求めるのかというと…[略]。さらに、\(\pi_k\)から\(\tilde{\pi}_k\)をどうやって逆算すんのかというと…[略]。
[いやー、山ほど研究があるんですね。こういうのはもはやソフトにお任せでいいのかも]
2.8 順序抽出
母集団の全ユニットについて正のサイズ変数\(x_k \gt 0\)が既知で、これに比例させて目標とする一次包含確率\(\pi_k\)を決めたとする。\(N\)個の\([0,1]\)一様乱数\(\omega_k\)を生成して(別に一様乱数でなくてもいいんだけど)、\(\omega_k / \pi_k\)が小さい順に\(n\)個を選ぶ。簡単でよろしいが、あいにく包含確率は厳密には\(\pi_k\)じゃない。
3. 非復元不均一確率抽出における点推定
[この節にいちばん関心があるので、少し細かくメモする]
Rao(1966) いわく、包含確率と関連していない特性の母合計について推定する場合、重みづけのない推定量$$ \hat{Y}_u = \frac{N}{n} \sum_{k \in s} y_k$$ を使ったほうがよい。デザインバイアスは$$ bias(\hat{Y}_u) = \frac{N}{n} \sum_{k \in U} y_k \pi_k – \sum_{k \in U} y_k = \frac{N^2}{n} \frac{1}{N} \sum_{k \in U} \left( y_k – \frac{\hat{Y}_u}{N} \right) \left( \pi_k – \frac{n}{N} \right)$$ つまり、\(y_k\)と\(\pi_k\)の共分散に比例するわけで、無相関ならバイアスがない。いっぽう分散について考えると、以下の超母集団モデルを\(\xi\)として$$ y_k = \mu + \epsilon_k $$ $$ E_\xi(\epsilon_k | \pi_k) = 0, E_\xi(\epsilon^2_k | \pi_k) = \sigma^2, E_\xi(\epsilon_k \epsilon_l | \pi_k) = 0$$ \(\xi\)の下での\(\hat{Y}_u, \hat{Y}_\pi\)の平均的な分散について考えてみると、前者のほうが小さい。
[標本平均とHT推定量の分散を比べる話、カジュアルにはなんどか目にしていたのだが、ちゃんと説明してあるのを探していたのである… ようやく見つけられてうれしいが、Rao(1966)というのがなんなのかがわからない。候補としては”On the variance of the ratio estimator for Midzuno-Sen sampling scheme”, Metrika. というのがあるんだけど、これかなあ。→ 探して読んでみたところ、どうやらこれではなさそう。
それにしても、超母集団モデルを持ち出さないと説明できない話なのだろうか。\(y_k\)と\(\pi_k\)の分散と共分散に応じて勝ち負けが決まるという話ではないのかしらん。]
Amahia et al.(1989) は次の推定量について検討している。\(y_k\)と\(\pi_k\)の観察された相関を\(\rho\)として$$ \hat{Y}_a = (1-\rho) \hat{Y}_u + \rho \hat{Y}_\pi$$ [へええ！ これ面白いアイデアだなあ。初めて聞いたけど、なんで使われてないんだろうか。おそらく、Amahia, Chaubey, Rao (1989) “Efficiency of a new PPS sampling for multiple characteristics” J.Stat.Plan.Inference.]
Hajek(1971)の推定量 $$ \hat{Y}_H = N \left( \sum_{k \in s} \frac{1}{\pi_k} \right)^{-1} \sum_{k \in s} \frac{y_k}{\pi_k} $$もよくつかわれている。近似的にデザイン不偏である。\(y_k\)と\(\pi_k\)の相関がないとき、上の超母集団モデルのもとで、分散が小さい。Sarndal 本のp.258をみよ。[ご親切にありがとう]
Hajek推定量が良く使われている理由は、重み付き平均になっているからで、とくにカウントの推定(合計が所与の定数にならないとおかしい)の際に便利である。
というわけで、\(y_k\)と\(\pi_k\)の間に相関があるときは\(\pi\)推定量が、相関がないときは\(\hat{Y}_u, \hat{Y}_H\)がよい。Basuの挙げた有名な事例を挙げると…[例のサーカスの象の話]。
ところで、Sen-Yates-Grundyの分散推定量$$ \hat{var}(\hat{Y}_\pi) = \frac{1}{2} \sum_{k \in s} \sum_{l \in s} \frac{\pi_{kl} – \pi_k \pi_l}{\pi_{kl}} \left( \frac{y_k}{\pi_k} – \frac{y_l}{\pi_l} \right)^2 $$は\(\hat{Y}_u, \hat{Y}_a, \hat{Y}_h\)のいずれにも使える。式のなかの\(y_k\)を\(y_k \pi_k N/n\)にすれば\(\pi_{kl} \gt 0\)のときの\(\hat{Y}_u\)の分散のデザイン不偏推定量になるし、\(y_k \pi_k /(n/N(1-\rho) + \rho \pi_k\)にすれば\(\pi_{kl} \gt 0\)のときの\(\hat{Y}_a\)の分散の近似的にデザイン不偏な推定量になるし、\(y_k – \hat{Y}_H\)にすれば\(\pi_{kl} \gt 0\)のときの\(\hat{Y}_H\)の分散の近似的にデザイン不偏な推定量になる。[へええ。\(n\)固定の時という限定がつくと思うけど、面白い豆知識だ]
サイズ変数は、関心ある変数とサイズ変数の相関が高くなるように選ぶべきである。もっとも実務的には、関心ある変数は複数あって、なかには相関がない変数もあるかもしれないということが多い。そういう場合には単純平均かHajek推定量がよろしかろう。
4. 同時包含確率を使わない分散推定
同時包含確率を求めるのは大変なことが多いし、Sen-Yates-Grundy推定量は計算が大変だし、同時包含確率をデータにつけるのも大変だ。
層内で確率不均一な一段層別抽出での、\(\pi\)推定量\(\hat{Y}_\pi\)の分散推定量について考えよう。層の数を\(H\), 層\(h\)への所属インジケータを\(z_{kh}\)とする。WLS回帰係数$$ \hat{B}_k = \sum_{k \in s} \lambda_k z_kh \frac{y_k}{\pi_k} / \left( \sum_{k \in s} \lambda_k z^2_{kh} \right)^{-1} $$ を求め、残差$$\hat{e}_k = \frac{y_k}{\pi_k} – \sum_{h=1}^H \hat{B}_h z_kh$$ を求めれば、$$ \hat{var}^*(\hat{Y}_\pi) = \sum_{k \in s} \alpha_k \hat{e}^2_k $$ という近似ができる。\(\alpha_k, \lambda_k\)をいい感じに決められればの話だが。
いくつかのアイデアを紹介しよう。
- \(\alpha_k = \lambda_k = 1\)。これは復元あり抽出の下での分散推定量に帰着する[上でメモを省略したけど、\(Y\)のHansen-Hurwitz推定量の分散]。抽出割合が大きいとき、分散の過大評価となる。
- \(\alpha_k = 1- \frac{\pi_k(n_h-1)}{n_h}, \lambda_k = 1\)。Hartley-Rao(1962)の分散推定量というのに帰着する。ランダム化系統抽出のもとでの一致推定量である。[…中略…] 抽出割合が小さいか、逆にすごく大きいときにおすすめ。
- \(\alpha_k = \lambda_k = \frac{(1-\pi_k)n_h}{n_h-1}\)。Hajek(1964)の分散推定量というのに帰着する。[…中略…] \(\sum_{l \in U_h} \pi_l(1-\pi_l) \rightarrow \infty\)という仮定の下で最大エントロピー抽出の分散の近似となる。
- […後略…]
[結局最後まで層別抽出の話だった。包含確率が個体ごとにちがうんだけど同時包含確率がわからないという場合はどうすればいいのだろうか。というか、ソフトウェアはどうしているのだろうか。だって確率不均一抽出のデータってそういうのがほとんどですよね、ローデータにウェイトはついているけど抽出デザイン情報はわからないという]
5. 平均の関数の分散推定
関心あるパラメータ\(\theta\)が、\(Q\)個の調査変数の有限母集団平均の関数になっていることがある。\(\theta = g(\mu_1, \ldots, \mu_Q)\), ただし\(g(\cdot)\)は平滑で微分可能、という場合である。比とか下位母集団平均とか相関とか回帰とかがそうですね。いっぽうL統計量[なにそれ]とかロジスティック回帰とかはここにははいらない。
その推定量は、Hajek推定量を\(\hat{\mu}_{qH}\)として[Hはインデクスじゃないぞ、気をつけろ]、\(\hat{\theta} = g(\hat{\mu}_{1H}, \ldots, \hat{\mu}_{QH})\)である。[いや、そうですか？ それは\(g(\cdot)\)によっては使い物にならないのではないか]
\(\hat{\theta}\)の分散推定量としては、まず線形化分散推定量がある。[…中略…]
慣習的なジャックナイフ分散推定は、非復元確率不均一抽出のもとではかならずしも一致性がない。一般化ジャックナイフ分散推定量というのがあって[…中略…]。
項目欠損の単一代入法の場合、Rao-Shao法というのがあって[…後略…]。
6. 釣り合い抽出
6.1 定義
補足変数たちの\(\pi\)推定量がそれらの既知の母合計と等しいとき、釣り合いが取れたデザイン balanced design であるという。抽出デザインのなかにカリブレーションが埋め込まれているといってもいい。[なるほど]
釣り合い抽出の一般的な方法にcube法がある。分散推定量もある。SASやRに載っている。
6.2 釣り合い抽出とcube法
補足変数を\(x_1, \ldots, x_p\)とし、ユニット\(k\)の値のベクトルを\(\mathbf{x}_k\)とする。所与の包含確率\(\pi_k\)について、デザイン\(p(.)\)が下式を満たすことを釣り合いという。$$ \sum_{k \in s} \frac{\mathbf{x}_k}{\pi_k} = \sum_{k \in U} \mathbf{x}_k$$ \(\mathbf{x}_k = \pi_k\)ならば固定サイズ制約のことである。\(\mathbf{x}_k\)が層のインジケータなら層別抽出のことである。
実現不可能な場合もある。たとえば上の式の右辺が整数だったりするとそうなりやすい。
cube法について。
ある標本を包含インジケータのベクトル \(\mathbf{s} = (s_1, \ldots, s_N) \in \mathbb{R}^N\)で表現すると、これは\(\mathbb{R}^N\)における\(N\)次元立方体の\(2^N\)個の頂点とみることができる。あるデザインとは、それぞれの頂点に確率\(p(.)\)を付与して、標本の期待値が一次包含確率のベクトル\(\mathbf{\pi}\)になるようにすることだといえる。
[面白いことを考える人もいるものだ… めんどくさいから省略するけど、cube法というのはflightフェイズとlandingフェイズに分かれるのだそうである。前者では立方体のなかをランダムウォークしてどこかの頂点で止まる。landingフェイズではそこから制約された下位空間に着地する。それぞれ複数のアルゴリズムがあるのだそうな。
まあよほど暇になったら勉強することにして、ユーザの観点から見ると、有限母集団の全ユニットについて希望の包含確率があって、かつ標本における共変量の分布が母集団と一致していてほしいときに、何だか知らんがうまいことやってくれる方法がある、というわけだ。具体的な使い道が思いつかないけれど、面白いなあ]
————-
勉強になりましたですー。
あれこれ調べ物をする前に、最初にこういう解説を読むべきなのかもしれないけど、でもほんとに最初に読んじゃったらありがたみがわかないんですよね。