2010-12-01

相関のある２つの確率変数の同時分布と和の分布

R 分布確率密度分布

こちらの記事では、その冒頭に「怪しい」と書いた通り、怪しい話だった。
とは言え、いくつか、このことを考えている
考えを進める上で、「和」の分布について考えることは、「同時分布」について考えることの、一部であること
一部ではあるけれども、「同時分布」について考えるのと、同じ道筋で「和」の分布を考えなくてはならないかどうかは、別であること(「和」の分布だけなら、何かしら、特殊事情がありそうに思われる)
同時分布・和の分布のいずれにしろ、個々の確率変数について平均と分散によって、標準化してから、とりかかるのがよさそうなこと
個々の変数が、平均に関して対称であるかどうかも大きな要素であること
正規分布で当てはまる話(軸を直交軸から斜交軸に変える)は、非対称分布には(少なくともそのままでは)当てはめられないこと
同じく、正規分布で当てはまる話(直交・斜交軸)は、対称分布であっても(少なくともそのままでは)当てはめられないこと
さらに、この当てはめられなさは、対称性な安定分布の場合も「当てはめられないまま」であること
確率密度関数・特性関数が安定分布のそれであり、かつ、安定な分布の関数が平均と分散とでパラメタ化されている場合についてであるが、2変数が独立であれば、和の分布の平均は、２変数の平均の和、和の分布の分散は、２変数の分散の和であるから、密度関数・特性関数はすぐにわかること
もし、２変数が非独立であるときには、２変数の共分散と相関係数とは１対１対応であるから、和の分布の分散(２変数の分散と共分散の和)もわかり、平均は相変わらず、２変数の平均の和である。もし、２変数が非独立な場合にも、安定な分布の非独立な２変数の和が、同一の形であるのであれば(ここがわからない)、これまた、平均と分散とでパラメタ化されている場合においては、確率密度関数・特性関数はすぐにわかる
同時分布の条件を整理すると
- すべてのｘの値、すべてのｙの値について、ｙ軸・ｘ軸に関して積分した値は、２変数の相関によらず変わらないこと
  - それを満足させるには２変数がそろって対称性分布であるならば、(x,y),(-x,-y)についてある同じ値を増やし、(x,-y),(-x,y)について同じ値を減らすという操作は可能なこと
  - それは、「分割表の周辺度数を変えずにセルの値を変える操作」と同じことであること
この先、考えを進めるためには、確率密度関数を三角関数に分解したことを想いだし、三角関数は、対称性がよいので、三角関数の相関がうまく表せれば、うまくいくのではないか、と思ったりする

minXY<--10
maxXY<-10
kizami<-0.05
xlim<-ylim<-c(-2*minXY,2*maxXY)
x<-seq(from=-maxXY,to=maxXY,by=kizami)
y<-x

Xpr<-dnorm(x,mean=0,sd=1)
#Xpr<-dexp(abs(x),0.01)
Ypr<-dnorm(y,mean=0,sd=1)
#Ypr<-dexp(abs(y),0.01)
Xpr<-cos(x)
Ypr<-cos(y)
Zpr0<-Xpr %*% t(Ypr)
Zpr<-Zpr0
image(Zpr)

# 以下、おまけ


Niter<-100000
#r<-min(Zpr)*0.1
r<-0.1
for(xxx in 1:Niter){
#for(i in 1:(length(x)/2)){
	#for(j in (i+1):(length(y)/2)){
	tmp<-sample(1:(length(x)/2),2)
	i<-min(tmp)
	j<-max(tmp)
	
	theta<-atan(((length(x)+1)/2-j)/((length(x)+1)/2-i))*2
	#print(i)
	#print(j)
	
	#i2<-tmp[1]+1
	#j2<-tmp[2]+1
	#tmpx<-sample(1:length(x),2)
	#tmpy<-sample(1:length(y),2)
	#tmpx[1]<-tmp[1]
	#tmpx[2]<-(length(x)+1)-tmp[1]
	#tmpy[1]<-tmp[2]
	#tmpy[2]<-(length(x)+1)-tmp[2]
	#tmpy[1]<-tmpx[1]
	#tmpy[2]<-tmpx[2]

	tmpi<-length(x)-i
	tmpj<-length(y)-j
	
	xs<-c(i,i,tmpi,tmpi,j,j,tmpj,tmpj)
	ys<-c(j,tmpj,j,tmpj,i,tmpi,i,tmpi)
	pms<-c(1,-1,-1,1,1,-1,-1,1)
	sumPr<-0
	for(k in 1:length(xs)){
		sumPr<-sumPr+Zpr[xs[k],ys[k]]
	}
	delta<-sumPr*r/length(xs)*cos(theta)
	for(k in 1:length(xs)){
		Zpr[xs[k],ys[k]]<-Zpr[xs[k],ys[k]]+delta*pms[k]
	}
	#tmpx<-c(i,j)
	#tmpy<-c(i,j)
	#delta<-min(Zpr[tmpx,tmpy])*r
	#if(tmpx[1]*tmpy[1]+tmpx[2]*tmpy[2]-tmpx[1]*tmpy[2]-tmpx[2]*tmpy[1]<0){
		#r<--r
	#}
	#Zpr[tmpx[1],tmpy[1]]<-Zpr[tmpx[1],tmpy[1]]+delta
	#Zpr[tmpx[2],tmpy[2]]<-Zpr[tmpx[2],tmpy[2]]+delta
	#Zpr[tmpx[1],tmpy[2]]<-Zpr[tmpx[1],tmpy[2]]-delta
	#Zpr[tmpx[2],tmpy[1]]<-Zpr[tmpx[2],tmpy[1]]-delta
	#}
	image(Zpr)

}


par(mfcol=c(2,2))
image(Zpr0)

plot(Xpr,apply(Zpr,1,sum))
plot(Ypr,apply(Zpr,2,sum))

image(Zpr)

par(mfcol=c(1,1))

vv<-vv0<-0
for(i in 1:length(x)){
	for(j in 1:length(y)){
		vv<-vv+Zpr[i,j]*x[i]*y[j]
		vv0<-vv0+Zpr0[i,j]*x[i]*y[j]
	}
}

vv
vv0