Inuverse Sci. X Tech. Blog

← ブログ一覧

Landy-Szalay推定量①

#cosmology#statistics#retail

はじめに

Landy-Szalay推定量についての計算です。小売の店舗の分布を解析するのに使いたいので計算したメモの第一弾です。

二点相関関数

二点相関関数の定義を述べる前に、前準備しておこう。

ある小さな領域ΔS\Delta Sに店舗が存在する確率がポアソン分布で記述できるとする。 平均店舗数密度をnnとすれば、ΔS\Delta Sにおける店舗数の期待値はnΔSn\Delta Sである。 ΔS\Delta Sに存在する店舗数をN(ΔS)N(\Delta S)書くとき、店舗数の確率分布は

P(N(ΔS)=k)=(nΔS)kk!e−nΔS\begin{align} P(N(\Delta S) = k) = \frac{(n \Delta S)^k}{k!} e^{- n\Delta S} \end{align}

で与えられる。ΔS\Delta Sに店舗が１つしかないとき、

P(N(ΔS)=1)=nΔS⋅e−nΔS\begin{align} P(N(\Delta S) = 1) = n \Delta S \cdot e^{- n \Delta S} \end{align}

となるので、ΔS≪1\Delta S \ll 1となるようにして、ΔS\Delta Sの１次までをとれば

P(N(ΔS)=1)=nΔS+O(ΔS2)\begin{align} P(N(\Delta S) = 1) = n \Delta S + \mathcal{O}(\Delta S^2) \end{align}

である。あるΔS1\Delta S_1に１店舗、あるΔS2\Delta S_2に1店舗ある場合ある場合、それらの同時確率は

ΔP=n2ΔS1ΔS2\begin{align} \Delta P = n^2 \Delta S_1 \Delta S_2 \end{align}

である。もし店舗が存在する位置同士に相関があって、ポアソンよりも店舗が見つかりやすくなる（見つかりにくくなる）効果としてξ(r)\xi(r)を考えると、

ΔP=n2(1+ξ(r))ΔS1ΔS2\begin{align} \Delta P = n^2 (1 + \xi(r)) \Delta S_1 \Delta S_2 \end{align}

となり、このξ(r)\xi(r)を二点相関関数という。ξ(r)>0\xi(r) > 0のとき、店舗同士はクラスターを形成しやすい傾向にあり、逆にξ(r)<0\xi(r) < 0ならば店舗通しは離れやすい性質を持つ。ξ(r)=0\xi(r) = 0ならばポアソン分布に従い店舗が配置される。

二点相関関数の推定量

二点相関関数の推定量は昔から調べられている。店舗間の距離rrからr+Δrr + \Delta rのビンの中に存在する店舗間距離のペアの数をDD(r)DD(r)とし、ランダムサンプルから得られるデータ点に対して、同様にペアの数をカウントしたものをRR(r)RR(r)とする。いま見ている相関とはランダムだった時に比べて、どのスケールでどれくらい集まりやすいか（離れやすいか）を表しているので、相関関数の推定量は単純には

ξ^(r)=DD(r)RR(r)−1\begin{align} \hat{\xi}(r) = \frac{DD(r)}{RR(r)} - 1 \end{align}

と表される。しかしながら、通常はポアソン分布の分散の最小値に近くなLandy-Szalay推定量

ξ^LS(r)=DD(r)−2DR(r)+RR(r)RR(r)\begin{align} \hat{\xi}_{\rm LS}(r) = \frac{DD(r) - 2DR(r) + RR(r)}{RR(r)} \end{align}

が使われる。

Landy-Szalay推定量についての計算

DD(r)DD(r)は観測することで得た標本のデータ点のペアカウントであるので確率変数である。RR(r)RR(r)はランダムサンプル同士のペアカウントで、DR(r)DR(r)は実データとランダムサンプル同士のペアカウントである。よく使われる推定量は

1+ξ^1(r)=DD(r)RR(r)nr(nr−1)n(n−1) ,1+ξ^2(r)=DD(r)DR(r)2nrn−1\begin{align} 1 + \hat{\xi}_1(r) &= \frac{DD(r)}{RR(r)} \frac{ n_r (n_r - 1) }{ n(n - 1) } \,, \\ 1+ \hat{\xi}_2(r) &= \frac{DD(r)}{DR(r)} \frac{ 2 n_r }{ n - 1 } \end{align}

となる。 観測されるデータ点の数がnnで、ランダムサンプルの数がnrn_rとするとき、DDDDのペアカウントの総数はn(n−1)/2n(n - 1) / 2であり、同様にRRRRのペアカウントの総数はnr(nr−1)/2n_r (n_r - 1) / 2である。 一方、DRDRのペアカウントはnnrn n_rである。 右辺のn,nrn, n_rで書かれる因子は規格化因子である。

データがnn点取得できて、かつランダムサンプルがnrn_rの時において、ペアカウントがDD(r)DD(r)となる条件付き確率をP(DD(r)∣n,nr)P(DD(r)|n, n_r)とすると、期待値からのゆらぎをα,β,γ\alpha, \beta, \gammaとして、それぞれ

DD(r)=E[DD(r)](1+α(r)) ,DR(r)=E[DR(r)](1+β(r)) ,RR(r)=E[RR(r)](1+γ(r))\begin{align} DD(r) &= \mathbb{E}[DD(r)] (1 + \alpha(r))\,, \\ DR(r) &= \mathbb{E}[DR(r)] (1 + \beta(r))\,, \\ RR(r) &= \mathbb{E}[RR(r)] (1 + \gamma(r)) \end{align}

と定義する。この定義より自動的にE[α]=E[β]=E[γ]=0\mathbb{E}[\alpha] = \mathbb{E}[\beta] = \mathbb{E}[\gamma] = 0である。またランダムサンプルの数は原理的には無限に増やすことができるので、任意の精度までγ\gammaを小さくできる。

例えば、ξ^2(r)\hat{ \xi }_2(r)の期待値と分散を求める。

ξ^2(r)=E[DD]E[DR]2nrn−11+α1+β−1\begin{align} \hat{\xi}_2(r) &= \frac{ \mathbb{E}[DD] }{ \mathbb{E}[DR] } \frac{ 2 n_r }{ n - 1 } \frac{ 1 + \alpha }{ 1 + \beta } - 1 \end{align}

であり、ここで微少量によらない部分の係数をA=E[DD]E[DR]2nrn−1A = \frac{ \mathbb{E}[DD] }{ \mathbb{E}[DR] } \frac{ 2 n_r }{ n - 1 } とする。(1+β)−1=1−β+β2+O(β3)(1 + \beta)^{-1} = 1 - \beta + \beta^2 + \mathcal{O}(\beta^3)であるので、ゆらぎの２次まで保つと

1+α1+β=(1+α)(1−β+β2+O(β3))=1+α−β−αβ+β2+O(δ3)\begin{align} \frac{ 1 + \alpha }{ 1 + \beta } &= (1 + \alpha) (1 - \beta + \beta^2 + \mathcal{O}(\beta^3)) \notag \\ &= 1 + \alpha - \beta - \alpha \beta + \beta^2 + \mathcal{O}(\delta_3) \end{align}

となる。ここでδ3\delta^3はゆらぎの３次という意味であり、これにあわせてδ1=α−β,δ2=−αβ+β2\delta_1 = \alpha - \beta, \delta_2 = - \alpha \beta + \beta^2とおくと、相関関数ξ^2(r)\hat{\xi}_2(r)は

ξ^2(r)=(A−1)+Aδ1+δ2+O(δ3)\begin{align} \hat{\xi}_2(r) &= (A - 1) + A \delta_1 + \delta_2 + \mathcal{O}(\delta_3) \end{align}

と書くことができる。これの期待値をとると

E[ξ^2(r)]=(A−1)+AE[δ1]+AE[δ2]+O(δ3)\begin{align} \mathbb{E}[\hat{\xi}_2(r)] &= (A - 1) + A \mathbb{E}[\delta_1] + A \mathbb{E}[\delta_2] + \mathcal{O}(\delta_3) \end{align}

であるが、１次のゆらぎの期待値は00なので、結局

1+E[ξ^2(r)]=A(1+E[δ2])+O(δ3)=E[DD]E[DR]2nrn−1(1−E[αβ]+E[β2])+O(δ3)\begin{align} 1 + \mathbb{E}[\hat{\xi}_2(r)] &= A(1 + \mathbb{E}[\delta_2]) + \mathcal{O}(\delta_3) \\ &= \frac{ \mathbb{E}[DD] }{ \mathbb{E}[DR] } \frac{ 2 n_r }{ n - 1 } (1 - \mathbb{E}[\alpha \beta] + \mathbb{E}[\beta^2]) + \mathcal{O}(\delta_3) \end{align}

を得る。

一方、分散を求めるためにξ^2(r)\hat{\xi}_2(r)の２次のモーメントを計算すると、

(ξ^2(r))2=((A−1)+Aδ1+Aδ2+O(δ3))2=(A−1)2+2A(A−1)δ1+A2δ12+2A(A−1)δ2+O(δ3)\begin{align} (\hat{\xi}_2(r))^2 &= \big( (A - 1) + A \delta_1 + A \delta_2 + \mathcal{O}(\delta_3) \big)^2 \\ &= (A - 1)^2 + 2 A (A - 1) \delta_1 + A^2 \delta_1^2 + 2 A (A - 1) \delta_2 + \mathcal{O}(\delta_3) \end{align}

なので

E[(ξ^2(r))2]=(A−1)2+A2E[δ12]+2A(A−1)E[δ2]+O(δ3)\begin{align} \mathbb{E}[(\hat{\xi}_2(r))^2] &= (A - 1)^2 + A^2 \mathbb{E}[\delta_1^2] + 2 A (A - 1) \mathbb{E}[\delta_2] + \mathcal{O}(\delta_3) \end{align}

である。ようやく分散が計算できて

Var[ξ^2(r)]=E[(ξ^2(r))2]−(E[ξ^2(r)])2=(A−1)2+A2E[δ12]+2A(A−1)E[δ2]−((A−1)+AE[δ2])2+O(δ3)=(A−1)2+A2E[δ12]+2A(A−1)E[δ2]−(A−1)2−2A(A−1)E[δ2]+O(δ3)=A2E[δ12]+O(δ3)=(E[DD]E[DR]2nrn−1)2(E[α2]+E[β2]−2E[αβ])+O(δ3)\begin{align} {\rm Var}[\hat{\xi}_2(r)] &= \mathbb{E}[(\hat{\xi}_2(r))^2] - (\mathbb{E}[\hat{\xi}_2(r)])^2 \notag \\ &= (A - 1)^2 + A^2 \mathbb{E}[\delta_1^2] + 2 A (A - 1) \mathbb{E}[\delta_2] - \big( (A - 1) + A \mathbb{E}[\delta_2] \big)^2 + \mathcal{O}(\delta_3) \notag \\ &= (A - 1)^2 + A^2 \mathbb{E}[\delta_1^2] + 2 A (A - 1) \mathbb{E}[\delta_2] \notag \\ &\qquad - (A - 1)^2 - 2A(A - 1)\mathbb{E}[\delta_2] + \mathcal{O}(\delta_3) \notag \\ &= A^2 \mathbb{E}[\delta_1^2] + \mathcal{O}(\delta_3) \notag \\ &= \left( \frac{ \mathbb{E}[DD] }{ \mathbb{E}[DR] } \frac{ 2 n_r }{ n - 1 } \right)^2 ( \mathbb{E}[\alpha^2] + \mathbb{E}[\beta^2] - 2 \mathbb{E}[\alpha\beta] ) + \mathcal{O}(\delta_3) \end{align}

である。

Landy-Szalay推定量についての計算２：

定義より、

α=DD−E[DD]E[DD] ,β=DR−E[DR]E[DR] ,\begin{align} \alpha &= \frac{DD - \mathbb{E}[DD]}{\mathbb{E}[DD]}\,, \\ \beta &= \frac{DR - \mathbb{E}[DR]}{\mathbb{E}[DR]}\,, \\ \end{align}

であるから、２次のモーメントがDD,DR,RRDD, DR, RRのペアカウントのみで表現できて

E[α2]=E[DD2]−(E[DD])2(E[DD])2=Var[DD](E[DD])2 ,E[β2]=E[RR2]−(E[RR])2(E[RR])2=Var[RR](E[RR])2 ,E[αβ]=E[DD⋅DR]−E[DD]⋅E[DR]E[DD]⋅E[DR]=Cov[DD,DR](E[DD])2 ,\begin{align} \mathbb{E}[\alpha^2] &= \frac{ \mathbb{E}[DD^2] - (\mathbb{E}[DD])^2 }{ (\mathbb{E}[DD])^2 } = \frac{ {\rm Var}[DD] }{ (\mathbb{E}[DD])^2 }\,, \\ \mathbb{E}[\beta^2] &= \frac{ \mathbb{E}[RR^2] - (\mathbb{E}[RR])^2 }{ (\mathbb{E}[RR])^2 } = \frac{ {\rm Var}[RR] }{ (\mathbb{E}[RR])^2 }\,, \\ \mathbb{E}[\alpha \beta] &= \frac{ \mathbb{E}[DD \cdot DR] - \mathbb{E}[DD] \cdot \mathbb{E}[DR] }{ \mathbb{E}[DD] \cdot \mathbb{E}[DR] } = \frac{ {\rm Cov}[DD, DR] }{ (\mathbb{E}[DD])^2 }\,, \end{align}

と表される。

参考文献

実は…

最近オーストラリアに一週間ほど滞在しました。

オーストラリアには二つの小売寡占企業がありまして、その名前を「Coles」「Woolworths」と言います。私もとある小売系のシステムに携わっているので、一応関係者の端くれとしてお店を観察したり、シェアハウスの人に聞き込みをしてみたのです。

それによれば、

  • オーストラリアの商品の入れ代わり頻度は日本よりもずっと小さい
  • オーストラリアの方がDC（配送センター）から店舗までの距離が大きい

らしいのです。私はこのような推測を実際調べたくなる性で、また調べている途中に気になるテーマが変わってしまう性でもあります。今回も上の二つを調べていたところ、日本とオーストラリアの店舗の分布が統計的に同質のものなのか、そうではないものなのかが気になってしまったのです。

こういうわけで気づけば宇宙論の論文を開いて計算を追っています。いつ計算が終わることやら…。

← ブログ一覧