ポアッソン分布の分解

• さて、炎症性関節総数がポアッソン分布に従うときに、右と左の炎症性関節数はどんな分布になるのだろうか
• それに進む前に、多変量ポアッソン分布についてリンクを張っておく。多変量正規分布があるように、複数の軸について、それぞれポアッソン分布に従うような多変量を扱うことができて、それが多変量ポアッソン分布(Mathematicaのページがこちら)
• では、単変量ポアッソン分布に従うような炎症関節数があるときに、左右の振り分けが独立であったとしたら、左右のその数は、どんな風になるだろうか。分布は、ポアッソン分布と似た感じの非対称な離散的１峰性分布になる。それを確かめるのが次のRのソース。ただし、k次元に一般化して処理している。次元数(k)を２にしておけば、左右の話しになる。

k<-2 # 次元 変数の数
n<-30 # 考慮する、次元ごとの値の最大値
npoly<-6 # 全次元について足し合わせ、単変量ポアッソン分布を仮定したときの平均値
ttt<-dpois(0:n,npoly) # 単変量ポアッソン分布の確率密度分布

aaa<-array(0,rep((n+1),k)) # 計算結果格納用アレイ
keta<-0:(k-1) # n進数計算のための道具
keta<-(n+1)^keta # n進数のアレイ内番地計算用の道具

counter<-rep(0,k) # 同上
counter[1]<--1
loop<-TRUE
while(loop){
	counter[1]<-counter[1]+1
	for(i in 1:(k-1)){
		if(counter[i]==n+2){
			counter[i]<-0
			counter[i+1]<-counter[i+1]+1
		}
	}
	banchi<-sum(keta*counter)+1
	vvv<-counter+1
	vvv2<-sum(vvv)-k+1
	if(vvv2<=n+1){
		aaa[banchi]<-ttt[vvv2]*exp(-sum(lgamma(vvv)))
	}
	if(sum(counter)==(n+1)*k){
		loop<-FALSE
	}
}
par(mfcol=c(2,2))
#image(aaa[,,1]) # k=3の場合
image(aaa)
plot(apply(aaa,1,sum),type="b",main="軸１の分布")
plot(apply(aaa,2,sum),type="b",main="軸２の分布")


ax<-apply(aaa,1,sum)
ax<-ax/sum(ax)
tmpmean<-sum((0:n)*ax)
plot(0:n,dpois(0:n,tmpmean),ylim=c(0,max(ax)))

par(new=TRUE)
plot(0:n,ax,type="l",col="red",ylim=c(0,max(ax)),main="1軸の分布に、その平均値をパラメタとして与えたポアッソン分布を重ねた")