#The chi-square test of goodness of fit for hypothesis on Benford's distribution. 
#Simulation evaluation of power (b) 
# on the basis of simulated distribution of the test statistic
# under fixed significance level (a) and sample size (n);
# ls - number of replication,
#start sample size:
n=141
ls=10000
a=0.05
#H0: p=p0; H1: p=p1
k=9
czas=proc.time()
#Benford's prob. function:
p0=matrix(0,k,1); 
for (i in 1:k) p0[i]=log10(1+1/i)
EX=matrix(1:9,1,9)%*%p0
#alternative distributions:
#p12=matrix(c(-.005,-.005,-.005,-.005,0,.005,.005,.005,.005),k,1)
#p12=matrix(c(-.01,-.01,-.01,-.01,0,.01,.01,.01,.01),k,1)
p12=matrix(c(-.04,-.03,-.02,-.01,0,.01,.02,.03,.04),k,1)
#p12=matrix(c(-.02,-.015,-.01,-.005,0,.005,.01,.015,.02),k,1)
p1=p0+p12
#p1=matrix(c(9/45,8/45,7/45,6/45,5/45,4/45,3/45,2/45,1/45),k,1)
#p1=matrix(1/k,k,1)
#p1=matrix(dbinom(0:8,9,0.13),k,1)
#p1=as.matrix(dexp(1:9,1/EX)/sum(dexp(1:9,1/EX)))

qs=matrix(0,ls,1)
ws=matrix(0,ls,1)

fitchi=function(n,p,w) 
{#evaluation of the test statistic
 n*sum((w-p)*(w-p)/p)}

bs=0; 
for (t in 1:ls) {ws=rmultinom(1,n,p0)/n; qs[t]=fitchi(n,p0,ws)}
qs=sort(qs)
wk=qs[floor((1-a)*ls)+1]
bs=0
for (t in 1:ls) 
  {ws=rmultinom(1,n,p1)/n;if(fitchi(n,p0,ws)>=wk) bs=bs+1}
bs=bs/ls
#sample size:
n
#power:
bs
#simulated critical value:
wk;
#significance level:
a
proc.time()-czas
alarm()