Ntrial<-100
Ma<-30
as<-seq(from=0,to=Ma,length.out=Ntrial)
N<-10000
k<-sample(2:10,1)
mr<-runif(1)*5
xs<-RandomSphere(k,n=N)
C<-xs[,1]
S<-sqrt(1-C^2)
p1<-p2<-rep(0,Ntrial)
for(xx in 1:Ntrial){
a<-as[xx]
tmp<-a^2-mr^2*S^2
tmp[which(tmp<0)]<-0
tmp2<-mr*C
tmp3<-sqrt(tmp)
tmp4<-tmp2-tmp3
tmp5<-tmp2+tmp3
tmp4<-abs(tmp4)
tmp5<-abs(tmp5)
R12<-matrix(c(tmp4,tmp5),ncol=2)
R1<-apply(R12,1,min)
R2<-apply(R12,1,max)
plot(R1,R2)
if(a<mr){
Q2<-0.5 + (pchisq(R1^2,k)+pchisq(R2^2,k,lower.tail=FALSE))/2
}else{
Q2<-(pchisq(R1^2,k,lower.tail=FALSE)+pchisq(R2^2,k,lower.tail=FALSE))/2
}
p1[xx]<-mean(Q2)
p2[xx]<-pchisq(a^2,k,mr^2,lower.tail=FALSE)
}
plot(p1,p2)