2013-03-03

下限のある分布の比較

決断 Decision Theory カイ分布カイ二乗分布超幾何関数ベッセル関数距離 R

決断
選択肢があるときに、それぞれの選択肢を選んだ際にどんな帰結があったかのデータを読み、それぞれの選択肢の背後にある分布を予想し、その予想分布を比較することで、選択肢を選ぶ、ということをやっている
カテゴリカルな場合については決断ルールが階層化した状態も含めて検討した(こちら)
量的帰結については、選択肢がもたらす帰結の分布を推定するところまでやった(そのあと、決断のために分布を比較して「決断」を下すルールが未決)(こちら)
下限のある分布に伴う決断は、「生命予後(あとどれくらい生きられるか)」の予測に基づく決断を含むので、やりたい
下限のある分布で、下限が０であれば、「距離」に関する分布となる。距離といえば、カイ分布・カイ二乗分布、それらの非心な分布がある
それについての基礎的な確認はこちらの３日間(その１、その２、その３)
量的帰結の場合には、正規分布・オプションでカーネル分布推定をした
- 観測値を期待値として持つ正規分布の合わせたものとして、全観測値からの分布を推定するものだった。個々の観測値に対応付ける正規分布の分散をどうするかが問題になるわけだが、それについては、リスク関数を設定しそれが最小になるようにしてあった
下限あり・距離としての帰結の場合には、「下限あり・距離」の分布を合わせたものにすることができそうだ
- 観測値を期待値とする分布で「下限あり・距離」の分布を持ってくる
- 分布を決めるパラメタで期待値以外のそれは決まらないから、それを調整するために、リスク関数を設定してそれが最小になるようにすればよいだろう
- 「下限あり・距離」の分布として次のようなものがすぐに思いつく
  - 非心カイ二乗分布(期待値は分布を定める２つのパラメタの和 $df+M$ )
  - 非心カイ分布(期待値は分布を定める２つのパラメタを使って複雑に決まる $\sqrt{\frac{\pi}{2}} L_{\frac{1}{2}}^{(\frac{df}{2}-1)}(-\frac{m^2}{2})$
  - 非心カイ分布の期待値を計算する関数mean_chi_nc(df,m): df は自由度、mは非心カイ分布を定めるもう一つのパラメタ(非心カイ分布が対応する正規分布の期待値)

mean_chi_nc <- function(df,m,N=100){
	sqrt(pi/2)*Laguerre_poly(1/2,df/2-1,-m^2/2,N)
}
Laguerre_poly <- function(n,alpha,x,N=100){
	a <- -n
	b <- alpha+1
	tmp <- chokika_series_zs(a,b,x,N)
	tmp <- tmp * exp(lgamma(n+alpha+1)-lgamma(n+1)-lgamma(alpha+1))
	tmp
}
pochhammer <- function(x,n){
	if(n==0){
		return(1)
	}else{
		return(prod(x+(0:(n-1))))
	}
}
chokika_single <- function(a,b,z,n){
	a.p <- rep(NA,length(a))
	b.p <- rep(NA,length(b))
	for(i in 1:length(a)){
		a.p[i] <- pochhammer(a[i],n)
	}
	for(i in 1:length(b)){
		b.p[i] <- pochhammer(b[i],n)
	}
	signz <- sign(z)^n
	tmp <- -lgamma(n+1)+n*log(abs(z))
	if(length(a) > 0){
		tmp <- tmp + sum(log(abs(a.p)))
	}
	if(length(b) > 0){
		tmp <- tmp - sum(log(abs(b.p)))
	}
	return(signz*prod(sign(a.p))*prod(sign(b.p))*exp(tmp))
}

chokika_series <- function(a,b,z,N = 100){
	ret <- 0
	for(i in 0:N){
		ret <- ret + chokika_single(a,b,z,i)
	}
	ret
}
chokika_series_zs <- function(a,b,zs,N=100){
	ret <- rep(NA,length(zs))
	for(i in 1:length(zs)){
		ret[i] <- chokika_series(a,b,zs[i],N)
	}
	ret
}

カイ分布とカイ二乗分布はどちらも２パラメタの分布。その２パラメタごとに期待値を計算するとこんな等高線が描ける

ここでやっている非心カイ分布を描いておこう
- 自由度1でパラメタmを0から20まで細かく振ってある

dchi <- function(x,df,m){
	y <- x^2
	M <- m^2
	G <- dchisq(y,df,M)
	G*2*x
}

df <- 1

ms <- seq(from=0,to=2,length=20)
chi.out <- matrix(NA,length(ms),length(x))

for(i in 1:length(ms)){
	chi.out[i,] <- dchi(x,df,ms[i])
}

matplot(x,t(chi.out),type="l")

期待値が同じで自由度が異なるカイ分布を描く(対応するカイ分布パラメタを推定してやる)

m <- 2
dfs <- seq(from=1,to=10,length=10)
chi.out <- matrix(NA,length(dfs),length(x))

for(i in 1:length(dfs)){
	out.optim <- optim_nc_m(dfs[i],m)
	# 推定されたm(非心カイ分布のパラメタ)の値
	print(out.optim$par)

	print(mean_chi_nc(dfs[i],out.optim$par))
	chi.out[i,] <- dchi(x,dfs[i],out.optim$par)
}

matplot(x,t(chi.out),type="l")

ある観測値が与えられたとき、期待値がその観測値となるような下限のある分布をとることを考えよう。カイ分布、カイ二乗分布のそれぞれについて、その観測値に相当する等高線上の自由度とMとの組み合わせで決まる、いろいろな分布がこの条件を満たすことがわかる
そのあまたある分布のどれがもっともよいかをリスク関数の最小化処理で選びたい
まず、カイ分布の方は期待値からパラメタを求めないといけないけれど、式が(描けてはいるけれど)おばけみたいな式なのでoptim()関数を使って指定の自由度に対するもう一つのパラメタ値を推定することにしよう。

df <- seq(from=0,to=5,length=20)
df<-df[-1]
M <- df

dfM1 <- dfM2 <- matrix(0,length(df),length(M))
for(i in 1:length(dfM1[,1])){
	for(j in 1:length(dfM1[1,])){
		dfM2[i,j] <- df[i]+M[j]
		dfM1[i,j] <- mean_chi_nc(df[i],M[j])
	}
}
par(mfcol=c(1,2))
contour(dfM1,xlab="df",ylab="M",main="カイ分布期待値")

contour(dfM2,xlab="df",ylab="M",main="カイ二乗分布期待値")

期待値と自由度を指定して、非心カイ分布のパラメタを推定しよう

optim_nc_m <- function(df,M,init.m=1,N=100){
	
	fn <- function(x){
		(sqrt(pi/2)*Laguerre_poly(1/2,df/2-1,-x^2/2,N)-M)^2
	}
	optim(init.m,lower=0,fn,method= "L-BFGS-B")
}
# df=1,期待値=2
df <- 1
M <- 2
# 期待値がMになるように推定する
out.optim <- optim_nc_m(df,M)
# 推定されたm(非心カイ分布のパラメタ)の値
out.optim$par
# 推定された非心カイ分布のパラメタ値を用いて、
#非心カイ分布の期待値を計算すると
#確かにMになっている
mean_chi_nc(df,out.optim$par)