Inuverse Sci. X Tech. Blog

← ブログ一覧

Landy-Szalay推定量②

#cosmology#statistics

はじめに

Landy-Szalay推定量についての計算です。小売の店舗の分布を解析するのに使いたいので計算したメモの第二弾です。 前回の記事Landy-Szalay推定量①も参考にしてください。

二点相関関数の計算③

ある領域Ω\Omegaに店舗がnn個あるとする。これを格子に分割し、小空間がKK個できたとしよう。十分にKKを大きくすれば、一つの格子にひとつの店舗を置く状況を実現できる。このとき、それぞれの格子の店舗の有無はBernoulli過程に従う1。 セルの中に店舗がある状態を11、ない状態を00として、ii番目のセルに対して次のような変数viv_iを導入しておく：

vi={1(店舗がセル内に存在する)0(店舗がセル内に存在しない)\begin{align} v_i = \begin{cases} 1 & (\text{店舗がセル内に存在する}) \\[8pt] 0 & (\text{店舗がセル内に存在しない}) \end{cases} \end{align}

いま、知りたいのはペアの数である。ビンr±dr/2r \pm dr/2に入る点の数である。そこで次のようなペア(i,j)(i, j)の両方のセルに店舗が存在し、かつビンr±dr/2r \pm dr/2の中に入るときに11、それ以外のとき00となる関数を

ϱijr={1(∥ri−rj∥=rij<r±dr/2)0otherwise\begin{align} \varrho^r_{ij} = \begin{cases} 1 & (\|r_i - r_j \| = r_{ij} < r\pm dr/2) \\[8pt] 0 & \text{otherwise} \end{cases} \end{align}

と置く。すると、データのペアカウントは

DD(r)=∑i<jvivjϱijr\begin{align} DD(r) &= \sum_{i < j} v_i v_j \varrho_{ij}^r \end{align}

と表すことができる。

セルがKK個のとき、nn店舗で固定された時のvivjv_i v_jの期待値を求めよう。vivjv_i v_jなので

E[vivj]=1⋅1⋅P(vi=1,vj=1)+0⋅(otherwise)=P(vi=1)⋅P(vj=1∣vi=1)=n(n−1)K(K−1)\begin{align} \mathbb{E}[v_i v_j] &= 1 \cdot 1 \cdot P(v_i = 1, v_j = 1) + 0 \cdot (\text{otherwise}) \notag \\[8pt] &= P(v_i = 1) \cdot P(v_j = 1 | v_i = 1) \notag \\[8pt] &= \frac{ n ( n - 1 ) }{ K ( K - 1 ) } \end{align}

である。よって、DD(r)DD(r)の期待値は

E[DD(r)]=E[∑i<jvivjϱijr]=∑i<jE[vivj]ϱijr=n(n−1)K(K−1)∑i<jϱijr\begin{align} \mathbb{E}[DD(r)] &= \mathbb{E} \left[ \sum_{i < j} v_i v_j \varrho_{ij}^r \right] \notag \\[8pt] &= \sum_{i < j} \mathbb{E}[ v_i v_j ] \varrho_{ij}^r \notag \\[8pt] &= \frac{ n ( n - 1 )}{ K ( K - 1 ) } \sum_{i < j} \varrho_{ij}^r \end{align}

となる。ここで、∑i<jϱijr\sum_{i < j} \varrho_{ij}^rはvi,vjv_i, v_jの値によらない幾何学的な関数である。幾何学的な関数ϱijr\varrho_{ij}^rはvi,vjv_i, v_jに依存するので、期待値の外側に出すことはできないのではないかと、読者は不安になるかもしれない。しかし実際は、ϱijr\varrho_{ij}^r自体はvi,vjv_i, v_jの値についてなんも関心をもたない。vi=1v_i = 1だろうがvj=0v_j = 0だろうがなんでもよく、セルが含まれる領域の形に依存するので、「幾何学的な関数」と呼ばれる。

分散を求めるためにはE[DD2]\mathbb{E}[DD^2]を計算する必要がある。単純に考えると、

E[DD2(r)]=(∑i<jE[vivj]ϱijr)2=∑i<j∑k<lE[vivjvkvl]ϱijrϱklr\begin{align} \mathbb{E}[DD^2(r)] &= \left( \sum_{i<j} \mathbb{E} [v_i v_j] \varrho_{ij}^r \right)^2 \notag \\[8pt] &= \sum_{i < j} \sum_{k < l} \mathbb{E} [v_i v_j v_k v_l] \varrho_{ij}^r \varrho_{kl}^r \end{align}

となる。これをこのまま計算するのは辛い。そこで、この式の意味を考える。これはi,ji, jのペアとk,lk, lのペアの二つに依存する。全てのペアの組み合わせの数は(K2)2=[K(K−1)/2]2\binom{K}{2}^2 = \left[ K(K-1)/2 \right]^2であるが、この内訳を考えると便利である。内訳は3つに分かれる：

  • (i). (i,j)(i, j)、(k,l)(k, l)の全ての点が違うときの二つのペアの組み合わせ
  • (ii). (i,j)(i, j)、(i,k)(i, k)のように、ある一点を共有する二つのペアの組み合わせ
  • (iii). (i,j)、(j,k)(i, j)、(j, k)が重なっているときの二つのペアの組み合わせ

(i), (ii), (iii)の分解を形式的に次のように書く：

∑i<j∑k<lϱijrϱklr=∑i,j,k,l∗ϱijrϱklr+∑i,j,k∗ϱijrϱikr+∑i<jϱij\begin{align} \sum_{ i < j } \sum_{ k < l } \varrho_{ ij }^r \varrho_{ kl }^r &= \sum_{ i, j, k, l }^\ast \varrho_{ ij }^r \varrho_{ kl }^r + \sum_{ i, j, k }^\ast \varrho_{ ij }^r \varrho_{ ik }^r + \sum_{ i < j } \varrho_{ ij } \end{align}

(i)の場合をまず考える。最初に選ぶペアの点の組み合わせの数は(K2)\binom{ K }{ 2 }である。 先に選んだ点と重複しないように注意すると、もう一つのペアの組み合わせは(K−22)\binom{K-2}{2}である。 よって、重複を許さない二つのペアの組み合わせの数は

(K2)(K−22)=K(K−1)(K−2)(K−3)4\begin{align} \binom{ K }{ 2 } \binom{ K - 2 }{ 2 } &= \frac{ K (K - 1) (K - 2) ( K - 3 ) }{ 4 } \end{align}

である。幾何学的な関数Gq(r)G_q(r)を導入すると、それを

∑i,j,k,l∗ϱijrϱklr=K(K−1)(K−2)(K−3)4Gq(r)\begin{align} \sum_{ i, j, k, l }^\ast \varrho_{ ij }^r \varrho_{ kl }^r &= \frac{ K (K - 1) (K - 2) ( K - 3 ) }{ 4 } G_q(r) \end{align}

と定義する。

(ii)を考えよう。(i)と変わらず、最初のペアの点の組み合わせの数は(K2)\binom{ K }{ 2 }である。 最初のペアともうひとつのペアはある一点で交わっているので、交わっていない点の選び方はK−2K - 2と分かる。さらに、交わる点の選び方は二つなので、(ii)のペアの組み合わせの数は

(K2)(K−2)×2=K(K−1)(K−2)\begin{align} \binom{K}{2} (K - 2) \times 2 &= K ( K - 1 ) ( K - 2 ) \end{align}

である。この場合の幾何学的な関数Gt(r)G_t(r)を導入し、

∑i,j,k∗ϱijrϱikr=K(K−1)(K−2)Gt(r)\begin{align} \sum_{ i, j, k }^\ast \varrho_{ ij }^r \varrho_{ ik }^r &= K ( K - 1 ) ( K - 2 ) G_t(r) \end{align}

と定義する。

最後に(iii)を対処する。二つのペアの両点が重なっているため、点の選び方はひとつのペアの点の選び方の組み合わせの数に等しい。よって、以前と同様に

∑i<jϱij=K(K−1)2Gp(r)\begin{align} \sum_{ i < j } \varrho_{ ij } &= \frac{ K ( K - 1 ) }{ 2 } G_p(r) \end{align}

である。

以上、(i)(ii)(iii)を統合して、上式を幾何学的な関数で表現すると

[K(K−2)2Gp(r)]2=K(K−1)(K−2)(K−3)4Gq(r)+K(K−1)(K−2)Gt(r)+K(K−1)2Gp(r)\begin{align} \left[ \frac{ K ( K - 2) }{ 2 } G_p ( r ) \right]^2 &= \frac{ K (K - 1) (K - 2) ( K - 3 ) }{ 4 } G_q(r) \notag \\ &\quad + K ( K - 1 ) ( K - 2 ) G_t(r) + \frac{ K ( K - 1 ) }{ 2 } G_p(r) \end{align}

となる。(i)(ii)(iii)で議論した組み合わせの数の議論から次の関係式も成り立つだろうと予想できる：

[K(K−2)2]2=K(K−1)(K−2)(K−3)4+K(K−1)(K−2)+K(K−1)2\begin{align} \left[ \frac{ K ( K - 2) }{ 2 } \right]^2 &= \frac{ K (K - 1) (K - 2) ( K - 3 ) }{ 4 } \notag \\ &\quad + K ( K - 1 ) ( K - 2 ) + \frac{ K ( K - 1 ) }{ 2 } \end{align}

この式は簡単に確かめることができる：

[K(K−2)2]2−K(K−1)(K−2)(K−3)4−K(K−1)(K−2)−K(K−1)2=K(K−1)4[K(K−1)−(K−2)(K−3)−4(K−2)−2]=K(K−1)4[K2−K−K2+5K−6−4K+8−2]=0\begin{align} &\left[ \frac{ K ( K - 2) }{ 2 } \right]^2 - \frac{ K (K - 1) (K - 2) ( K - 3 ) }{ 4 } \notag \\ &\quad - K ( K - 1 ) ( K - 2 ) - \frac{ K ( K - 1 ) }{ 2 } \notag \\ &= \frac{ K ( K - 1) }{ 4 } \left[ K ( K - 1) - (K - 2) ( K - 3 ) - 4 ( K - 2 ) - 2 \right] \notag \\ &= \frac{ K ( K - 1) }{ 4 } \left[ K^2 - K - K^2 + 5K - 6 - 4K + 8 -2 \right] = 0 \notag \\ \end{align}

やはり、(i)(ii)(iii)の分割の仕方は、「二つのペア」のペアが持つ点の選び方を網羅する。

データカウントの二次のモーメントは

E[DD2]=E[vivjvkvl]∑i,j,k,l∗ϱijrϱklr+E[vivjvk]∑i,j,k∗ϱijrϱikr+E[vivj]∑i<jϱij =n(n−1)(n−2)(n−3)K(K−1)(K−2)(K−3)∑i,j,k,l∗ϱijrϱklr+n(n−1)(n−2)K(K−1)(K−2)∑i,j,k∗ϱijrϱikr+n(n−1)K(K−1)∑i<jϱij=n(n−1)(n−2)(n−3)4Gq(r)+n(n−1)(n−2)Gt(r)+n(n−1)2Gp(r)\begin{align} \mathbb{E}[DD^2] &= \mathbb{E}[v_i v_j v_k v_l]\sum_{ i, j, k, l }^\ast \varrho_{ ij }^r \varrho_{ kl }^r + \mathbb{E}[v_i v_j v_k] \sum_{ i, j, k }^\ast \varrho_{ ij }^r \varrho_{ ik }^r + \mathbb{E}[v_i v_j]\sum_{ i < j } \varrho_{ ij } \notag \\\ &= \frac{ n ( n - 1 ) ( n - 2 ) ( n - 3 ) }{ K ( K - 1 ) ( K - 2 ) ( K - 3 ) } \sum_{ i, j, k, l }^\ast \varrho_{ ij }^r \varrho_{ kl }^r \notag \\ &\quad + \frac{ n ( n - 1 ) ( n - 2 ) }{ K ( K - 1 ) ( K - 2 ) } \sum_{ i, j, k }^\ast \varrho_{ ij }^r \varrho_{ ik }^r + \frac{ n ( n - 1 ) }{ K ( K - 1 ) } \sum_{ i < j } \varrho_{ ij } \notag \\ &= \frac{ n ( n - 1 ) ( n - 2 ) ( n - 3 ) }{ 4 } G_q(r) \notag \\ &\quad + n ( n - 1 ) ( n - 2 ) G_t(r) + \frac{ n ( n - 1 ) }{ 2 } G_p(r) \end{align}

セル数は一般的に非常に大きな数になる。セル数KKがK≫1K \gg 1であるとき、K4K^4で両辺を割ることで、

Gp(r)2≃Gq(r)\begin{align} G_p(r)^2 \simeq G_q(r) \end{align}

となる。

データカウントの期待値の分散を計算しよう。分散の定義より

Var[DD]=E[DD2]−(E[DD])2=n(n−1)(n−2)(n−3)4Gq(r)−[n(n−1)2Gp(r)]2+n(n−1)(n−2)Gt(r)+n(n−1)2Gp(r)\begin{align} {\rm Var}[DD] &= \mathbb{E}[DD^2] - (\mathbb{E}[DD])^2 \notag \\ &= \frac{ n ( n - 1 ) ( n - 2 ) ( n - 3 ) }{ 4 } G_q(r) - \left[ \frac{ n ( n - 1 ) }{ 2 } G_p(r) \right]^2 \notag \\ &\quad + n ( n - 1 ) ( n - 2 ) G_t(r) + \frac{ n ( n - 1 ) }{ 2 } G_p(r) \end{align}

となる。もし土地に散らばる店舗の数自体もまた期待値NNのPoisson分布

Po(n;N)=Nnn!e−N\begin{align} Po(n; N) = \frac{N^n}{n!} e^{-N} \end{align}

に従うのであれば、Var[DD]\mathrm{Var}[DD]のポアソンアンサンブル平均は次の結果に導ける。セル数が非常に大きい近似（K≫1K \gg 1）のもとで、

Epoisson[Var[DD]]≃N3Gt(r)+N22Gp(r)\begin{align} \mathbb{E}_{\rm poisson}[{\rm Var}[DD]] &\simeq N^3 G_t(r) + \frac{N^2}{2} G_p(r) \end{align}

となる2。

この式の導出のために、ポアソン分布の次の性質を用いた：

∑n=0n(n−1)⋯(n−k+1)⏟k個Po(n;N)=Nk\begin{align} \sum_{n = 0} \underbrace{n ( n - 1 ) \cdots (n - k + 1)}_{k\text{個}} Po(n;N) &= N^k \end{align}

参考文献

Footnotes

  1. Bernoulli過程では、複数の事象が連なって起こる時に、ある事象の確率変数はXiX_iは0,10, 1のみをとり、それぞれの事象が同じ確率ppで11の値をとる過程である。 ↩

  2. 宇宙論では銀河の観測数がポアソン分布に従う。しかし、店舗の解析でアンサンブル平均まで考えるのは冗長だと私は思う。なぜなら、現在観測している店舗の位置の統計的性質を知りたいからである。 ↩

← ブログ一覧