non-central parameter

• p値がある
• ncp=0なるカイ自乗分布に照らしたときの、カイ自乗値を出せる
• これが、実はncpが0ではない非心カイ自乗分布から計算されたものとして、その本当のp値が与えられたとする
• ncpを逆算しよう
• optim()を使ってやろう

# df=1 ncp=0にのっとるp値を作る
N<-10000
df<-1
Ps<-runif(N)

# 一部のp値をずらす
P0<-10^(-2)
a<-1.5

P1s<-Ps
these<-which(P1s<P0)
P1s[these]<-exp(a*(log(Ps[these])-log(P0))+log(P0))

plot(log(ppoints(N,a=0),10),log(sort(P1s),10),cex=0.1)
abline(a=0,b=1,col=2)

plot(Ps,P1s)

plot(sort(log(Ps)),sort(log(P1s)))

# Pの値に相当するncpを求めたい
Ncps<-rep(0,N)
	fr <- function(x) {   
    (pchisq(Chisq,df=df,lower.tail=FALSE,ncp=x)-P)^2
	}

for(i in 1:length(these)){
	Chisq<-qchisq(P1s[these[i]],df=df,lower.tail=FALSE)
	P<-Ps[these[i]]
	Ncps[these[i]]<-optim(c(2), fr)$par

}

ord<-order(Ps)
plot(Ps[ord],Ncps[ord])