折り返し

• ０で下限が切られている分布について、「正規分布」を折り返しましょうという話があった(こちら)
• 今日は「そもそもカイ分布は正規分布の折り返しなのだから、正規分布の折り返しとしてカイ分布が計算できる」という「そもそも」の話
• 正規分布の折り返しとして確率密度を計算する
• この場合はもちろん自由度が1

x <- seq(from=-10,to=10,by=0.01)
m <-1.3
y <- dnorm(x,m)
posi <- which(x>=0)
neg <- which(x<=0)
yposi <- y[posi]
yneg <- y[neg[length(neg):1]]
y1 <- yposi+yneg

plot(y1,type="l")

• カイ分布として同様の確率密度を計算する(Rではカイ二乗分布の関数があるので、それを使ってある)

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

y2 <- dchi(x[posi],df=1,m)
plot(y1,y2)

• 一致する
• 一致するけれど、期待値は異なります
• 折り返しを考慮してカーネル推定するときに(SASでの対応)、この折り返し関数(カイ分布)を積み重ねよう、そのときの分散(bw)をどうしましょうか、という話はあったようだが、そのときに「期待値が異なっている」部分をどうしましょう、という話はなかったようだが、はて？
• それはさておき、期待値をそろえた「折り返し」を作るのだけはやっておこう

m <-1.3
y <- dnorm(x,m)
posi <- which(x>=0)
neg <- which(x<=0)
yposi <- y[posi]
yneg <- y[neg[length(neg):1]]
y1 <- yposi+yneg

plot(y1,type="l")

y2 <- dchi(x[posi],df=1,m)
plot(y1,y2)


optim.out <- optim_nc_m(df=1,m)
y3 <- dchi(x[posi],df=1,optim.out$par)

plot(x[posi],y3,xlim=range(x[posi]),ylim=c(0,max(y1)))
par(new=TRUE)
plot(x[posi],y2,xlim=range(x[posi]),ylim=c(0,max(y1)),col=2)
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")
}

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
}