2012-07-19

独立な確率変数が２つ

２つのp値を併せることと統計量の和を取ること
- ２つのp値があったときに、それを総合的に評価したいことがある
- やり方はいろいろあるだろうけれども、わかりやすいのは、２つのp値の積をとること
- なぜならば、珍しさ $p_1$ と珍しさ $p_2$ とが観察されたとき、両方の観察の珍しさは(両者が独立であるならば)、 $p_1$ と $p_2$ の積だから
- $p.prod = p_1 \times p_2$
- この $p.prod$ の大きさは意味があるものの、このような値を取る統計量として考えたときには、この値の分布に照らして、小さい値は珍しく、大きい値はありきたり、という値に変えてやりたい。 $p.prod$ の値の分布に照らして、その累積確率としてのp値(これは一様)に対応付けたい、ということである
- ちなみに $p.prod$ の分布の様子は以下のように、一様分布とは全く違う

p1.null <- runif(10000)
p2.null <- runif(10000)

p.prod <- p1.null * p2.null

hist(p.prod)

- p値の積が使いたい、統計量の和として表したい、という２つの要請がある
- $p_1 \times p_2$ を何かしらの和にするには対数を取ればよく、そうすると $\log(p_1) + \log(p_2)$ となる
- ここで指数分布の特徴を思い出そう
  - 確率密度関数は $\lambda e^{-\lambda x}$
  - 累積分布関数は $1-e^{-\lambda x}$
  - p値は、1-累積分布なので、指数分布からのp値は $p = e^{-\lambda x}$ となる
  - 両辺の対数を取れば $\log(p) = -\lambda x$ ( $\lambda=1$ として簡単にしてしまおう)
- したがって、 $-(\log(p_1) + \log(p_2))$ は指数分布に従うような統計量 $x$ の和として表されることがわかる
- 言い換えよう
- $p_1 \times p_2$ が適当な順序をもたらすようなとき、指数分布を考えて、 $p_i$ に相当する指数分布の値の和をとるとよい
- 実際 $-2(\log(p_1)+\log(p_2)$ は２つの項のそれぞれが指数分布(自由度２のカイ自乗分布でもある)であり、この２つの自由度２のカイ自乗分布の和は自由度４のカイ自乗分布になって、一様分布であるp値を求めることは容易である

p1.null <- runif(10000)
p2.null <- runif(10000)
X <- 2*(-log(p1.null)  -log(p2.null))
hist(X)
hist(pchisq(X,4))

- - ちなみに、対数を取っても、自由度２のカイ自乗分布を用いても、p値から得られる統計量が等しいことをRで確認しておくと以下の通り

p1.null <- runif(10000)

df <- 2
st1 <- qchisq(p1.null,df,lower.tail=FALSE)
plot(st1,-2*log(p1.null))

２つのp値の併せ方を自由度で調整する
- 上記では、 $p_1 \times p_2$ の大小を守るために自由度２のカイ自乗分布を介してp値を統計量に換算し、その和を取った
- 異なる自由度のカイ自乗分布を介してp値を統計量に換算してやるとどうなるかを見てみよう
  - ２次元正規分布があり、そこに直交する自由度１のテストを実施し、その２つのp値を併せることを考える。
  - 一つは、p値を自由度２のカイ自乗分布で統計量化した上でその和をとり、それを一様分布pにするべく、自由度４のカイ自乗分布に照らしてp値にする( $p_1\times p_2$ の大小を守る方法)
  - もう一つは、p値を自由度１のカイ自乗分布で統計量化した上でその和をとり、それを一様分布にするべく、自由度２のカイ自乗分布に照らしてp値にする

X <- matrix(rnorm(10000*2),ncol=2)

phi <- pi/2
v1 <- c(cos(-phi/2),sin(-phi/2))
v2 <- c(cos(phi/2),sin(phi/2))

chi1 <- (X %*% v1)^2
chi2 <- (X %*% v2)^2

p1.null <- pchisq(chi1,1,lower.tail=FALSE)
p2.null <- pchisq(chi2,1,lower.tail=FALSE)

df <- 2
st1 <- qchisq(p1.null,df,lower.tail=FALSE)
st2 <- qchisq(p2.null,df,lower.tail=FALSE)

X2 <- st1 + st2

p.df2<-pchisq(X2,df*2,lower.tail=FALSE)
hist(p.df2)

df <- 1
st1 <- qchisq(p1.null,df,lower.tail=FALSE)
st2 <- qchisq(p2.null,df,lower.tail=FALSE)

X3 <- st1 + st2
p.df1<-pchisq(X3,df*2,lower.tail=FALSE)
hist(p.df1)
plot(p.df1,p.df2)

plot(X,col=X2)
plot(X,col=X3)

- - どちらの方法でやっても、最終的なp値の分布は一様分布になっている

- - ２つの方法で得られる統計量は少しずつ異なる

- - ２つの方法のそれぞれで得られる統計量が２次元正規分布の形とどのような感じになっているのかを描く
    - 自由度１を介する方法では、２次元正規分布の中心からの逸脱を評価している。方向は(たまたま２つの直交ベクトルを使って評価したが、それは方便であって、取り出したベクトルの向きそのものに意味はない)関係ない。他方、自由度２を介する方法は、２つのベクトルの間の方向を強調している

- - ちなみに、自由度１でやると円、自由度２でやると、対角線に尖りが出る。自由度をさらに大きくすれば、尖りが先鋭化し、自由度を１より小さくすれば、軸方向へ先鋭化し、対角方向がへこんでくる
独立な確率変数の和の平均と分散の期待値
- p値から出発して、それを併せて評価するために、なにがしかの統計量に変換してその和を取る話をしてきた
- 二つのp値(そして二つの統計量)は相互に独立であることｗ前提にして話をしてきた
- 一般に、独立な二つの確率変数の和では、その期待値と分散は、元の２つの確率変数の期待値の和と分散の和とになる(こちらを参照)
- このことは、独立な確率変数の和を考えているときに、期待値と分散を知ることは簡単であることを意味し、また、期待値と分散とがわかると一意に決まるような分布を取ることがわかっているときには、確率変数の和の分布そのものを知ることも容易であることを意味している
独立な確率変数の和とガンマ分布
- ２変数で決まる分布にガンマ分布がある
- ２変数で決まるので、期待値と分散が決まることで分布が確定する
  - その関係は $\theta k = \text{mean},\theta^2 k = \text{var}$
  - その確率密度関数は $\frac{x^{k-1}exp^{-\frac{x}{\theta}}}{\Gamma{k}\theta^k}$
- ２つのガンマ分布が独立であるとき、二つの確率変数があって、それぞれがこの分布のいずれかに従っているとき、確率変数の和の分布は、ガンマ分布であって(ガンマ分布とガンマ分布の和がガンマ分布になることを示さないといけないが…省略して)、その２変数は
  - $\theta_0 = \frac{\theta_1^2 k1 + \theta_2^2 k_2}{\theta_1 k_1 + \theta_2 k_2}$ :これは $\theta_1,\theta_2$ の重み付き平均
  - $k_0 = \frac{(\theta_1 k_1 + \theta_2 k_2)^2}{\theta_1^2 k_1 + \theta_2^2 k_2}$
  - ガンマ分布の和とパラメタの関係の確認

gammaTK <- function(a,b){
	t1 <- a[1]
	k1 <- a[2]
	t2 <- b[1]
	k2 <- b[2]
	t.sum <- (k1*t1^2+k2*t2^2)/(k1*t1+k2*t2)
	k.sum <- (k1*t1+k2*t2)^2/(k1*t1^2+k2*t2^2)
	return(c(t = t.sum,k=k.sum))
}
a <- runif(2)*4
b <- runif(2)*6

d <- gammaTK(a,b)

R1 <- rgamma(10000,shape = a[2],scale = a[1])
R2 <- rgamma(10000,shape = b[2],scale = b[1])
R3 <- rgamma(10000,shape = d[2],scale = d[1])

plot(sort(R1+R2),sort(R3))
abline(0,1)

ガンマ分布とカイ自乗分布と指数分布との関係
- カイ自乗分布は(でが半整数の)ガンマ分布
  - したがって、カイ自乗分布(に従う確率変数)の和は、そのようなガンマ分布(に従く確率変数)の和だから、 $(\theta_1 = 2, k_1)$ と $(\theta_2 = 2, k_2)$ というパラメタで表されるガンマ分布(であるところの２つのカイ自乗分布(自由度が $\frac{k_i}{2}$ のカイ自乗分布)(に従う確率変数)の和なので、 $\theta_0 = \frac{2^2 k_1 + 2^2 k_2}{2 k_1 + 2 k_2}=2$ 、 $k_0 = \frac{(2 k_1 + 2 k_2)^2}{2^2 k_1 + 2^2 k_2}=k_1+k_2$ となり、結局、２つのカイ自乗分布の自由度の和を自由度とするカイ自乗分布になることがわかる
- 指数分布は、 $k=1$ のガンマ分布
- 相等しい二つの指数分布は
  - 相等しい指数分布に従う確率変数の和が従う分布は $\theta_0 = \frac{\theta_1^2 \times 1 + \theta_1^2 \times 1}{\theta_1 \times 1 + \theta_1 \times 1}=\theta_1$ , $k_0 = \frac{(\theta_1 \times 1 + \theta_1 \times 1)^2}{\theta_1^2 \times 1 + \theta_1^2 \times 1}=2$ となり、結局、 $\theta$ は変えないで、 $k$ が２倍のガンマ分布になった
- 独立な確率変数の和は、「独立」なので、自由度が増える方向( $k$ が増える方向)で $\theta$ を変えない方向への変化になっている
- 一般にガンマ分布では、を同じくする分布２つについての確率変数の和は、を変えずにが２つのの和となるようなガンマ分布になる
  - Rで確かめる

n.iter <- 100
as <- bs <- ds <- matrix(0,n.iter,2)
for(i in 1:n.iter){
	a <- runif(2)*4
	b <- c(a[1],runif(1)*10)
	d <- gammaTK(a,b)
	as[i,] <- a
	bs[i,] <- b
	ds[i,] <- d
}

plot(data.frame(as,bs,ds,as[,2]+bs[,2]))