calc.FPR.chisq <- function(r,alpha=0.05){ # r=c(37,21,21,21) vs n=c(25,25,25,25)
m=length(r)
n=rep(mean(r),m)
p.val=chisq.test(rbind(r,n))$p.value
df=m-1
P0=n/sum(n)
P1=(r/n)/sum(r/n)
N=sum(r) # == sum(n)
ncp1=0
for(i in 1:m) ncp1=ncp1+N*(P1[i]-P0[i])^2/P0[i]
qcrit=qchisq(1-alpha,df,ncp=0)
curve(dchisq(x,df),0,20,xlab=quote(chi),ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
abline(v=qcrit,col='gray') ; text(qcrit,0,round(qcrit,2))
legend('topright',bty='n',legend=c('H0:Central','H1:Noncentral'),col=1:2,lty=1:2)
power=pchisq(qcrit,df,ncp1,lower=FALSE)
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR=y0/(y0+y1)
print(c(power=power,p.value=p.val,FPR=FPR),digits=3)
}

calc.FPR.chisq(c(10,16,34,9,10,26))
calc.FPR.chisq(c(1,0,9,2,1,1))