2012-07-17

ガンマ分布のパラメタと非独立の程度の関係

こちらとその前の日にとある論文での「２つの非独立統計量の和をガンマ分布でｐ値化する話」にまつわることをやっている
自由度２の観察空間を考える
２つの自由度１の検定を単位ベクトルで表現する
２つの検定の非独立の程度は、単位ベクトルのなす角 $\phi$ で表すことにする
この観察空間では、期待点を原点として、自由度２の正規分布が観察される
２つの自由度１の検定をカイ二乗統計量で行えば、その統計量の期待値は自由度に等しいから、1
２つの自由度１の検定統計量の和を、「合わせた統計量」とすれば、その期待値は、２つの自由度１検定が独立であっても非独立であっても、1+1=2
他方、２つの検定がまったく同一であれば、その自由度は１、まったく独立であれば、自由度は２
２つの統計量の和である「合わせた統計量」の等高線が楕円になることは、昨日示した通り
また、「合わせた統計量」の分布がガンマ分布に似ている(合致するのかもしれない(するような気がする)かどうかは引用文献のそのまた引用文献に戻らなくてはならないが、読んでないので、ここでは、どちらかだとしよう)
今、ガンマ分布は、期待値と分散とを定めると、ガンマ分布を決める２つのパラメタが決まってしまうような分布(パラメタの個数が２個の分布)
期待値は２つの自由度１検定の非独立度によらず２なのだから、ガンマ分布の２つのパラメタを決める作業は「１変数推定」に過ぎない…
推定してみると、 $\theta = 3 + \cos(2t)$ のような関係なのでは…と思われるが

# ２つの自由度１のカイ二乗統計量の和であるような統計量について

# いくつかの非依存度で統計量を取ってみよう
ts <- seq(from=0,to=1,length=10)*pi/2
# ガンマ分布の２つのパラメタを推定してみよう
thetas <- ks <- rep(0,length(ts))
# ガンマ分布の２つのパラメタを数式的に与えてみよう
thetas.2 <- 3+cos(2*ts)
# k * theta = 2 (今の設定でのガンマ分布の期待値)
ks.2 <- 2/thetas.2

# シミュレーションする観測点の数
n.pt <- 100000
# 統計量空間の次元
df <- 2
# df次元正規乱数
X <- matrix(rnorm(df*n.pt),ncol=df)
# この空間に複数の自由度１テストを表す単位ベクトルを置く
# テストは相互に非独立

# いろいろなtの値で実施
for(tt in 1:length(ts)){
	tti <- ts[tt]
	# ２つの検定ベクトル
	v1 <- c(1,0)
	v2 <- c(cos(tti),sin(tti))
	# この空間での乱点から、各テストベクトルへの射影の２乗はカイ二乗値(自由度１)
	chi2.1 <- (X%*%v1)^2
	chi2.2 <- (X%*%v2)^2

	# 自由度１でｐ値化する
	p.1 <- pchisq(chi2.1,df=1,lower.tail=FALSE)
	p.2 <- pchisq(chi2.2,df=1,lower.tail=FALSE)

	# Fisher'sの方法では自由度２で統計量に戻したが
	# df.2 として任意の自由度で統計量に戻し、その上で足し合わせよう
	df.2 <- 1
	x <- qchisq(p.1,df.2,lower.tail=FALSE)+qchisq(p.2,df.2,lower.tail=FALSE)

	# 統計量の標本平均と標本分散を出す
	m.x <- mean(x)
	v.x <- var(x)
	# ガンマ分布の平均はtheta k、分散はtheta^2 k なので
	# 標本平均・標本分散とから、theta,kの推定値を出す
	#theta <- v.x/m.x
	# 別のやり方としては、k * theta = 2 が「正しい」と考えて
	theta <- v.x/2
	k <- 2/theta
	#k <- m.x/theta
	# 値を格納しよう
	thetas[tt]<-theta
	ks[tt]<-k
	# 推定したガンマ分布パラメタを使って、ガンマ分布乱数を作ってみよう
	y <- rgamma(n.pt,scale=theta,shape=k)

	# 確かに、標本統計量との関係がとれている
	xlim <- ylim <- range(c(x,y))
	plot(sort(x),sort(y),xlim=xlim,ylim=ylim)
	abline(0,1,col=2)


	#plot(data.frame(p.1,p.2,p.3))
	plot(data.frame(p.1,p.2))
}

plot(data.frame(ts,thetas,ks))

# 推定値と、式で計算した値との関係は？
plot(thetas, thetas.2)