TangoT.index<-function(Nrep,freq,p,mc ){ # # --------------------------------------------------------------- # An example of execution # # Nrep<-999 # f100<-c(0,1,0,1,2,0,0,0,0,1,1,1,6,4,5,1,1,2,0,2,4,3) # out<-TangoT.index(Nrep,f100) # # --------------------------------------------------------------- # # ARGUMENTS # # Nrep : The number of Monte Carlo replications, e.g., 999, 9999 # freq : a vector of observed frequencies # For example, you can set : freq<-c(4,3,4,4,4,7,2,3,11,9,8,3) # p : a vector of multinomial parameters under null hypothesis # Default values are "equal probabilities" # mc : a user defined matrix of closeness measure A # Defaut values are exp( -|i-j| ) # # VALUES # # Stnd.C : standardized test statistic Z # C.p.value : Prob{C > c} using approximation (4.22) # C.u : a vector of percent contribution of Ui to C # Freq : a vector of simulated test statistics under the null # ------------------------------------------------------------------ nn<-sum(freq) lenn<-length(freq) if (missing(p)) p<-rep(1/lenn,lenn) if (missing(mc)){ mc<-matrix(0,lenn,lenn) for (i in 1:lenn) { for (j in 1:lenn) mc[i,j]<- exp( -abs(i-j) ) } } ac<-mc w1<-ac w2<-ac pp<-matrix(p) w<-diag(p)-pp%*%t(pp) q<-freq/nn for (i in 1:lenn) { for (j in 1:lenn) { w1[i,j]<- ac[i,j]*q[i] w2[i,j]<- ac[i,j]*(q[i]-p[i]) } } g<- q%*%ac%*%q eg<- p%*%ac%*%p+sum(diag( ac%*%w ))/nn vg<-( 4* p%*%ac%*%w%*%ac%*%p + 2/nn*sum(diag( ac%*%w%*%ac%*%w )) )/nn skew<-8*( 3* p%*%ac%*%w%*%ac%*%w%*%ac%*%p + + (1/nn)*sum(diag(ac%*%w%*%ac%*%w%*%ac%*%w)) ) / sqrt(nn) / (nn*vg)**1.5 df<-8/skew**2 hh<-(g-eg)/sqrt(vg) tc<-hh[1,1] pval<- 1-pgamma((df+tc*sqrt(2*df))/2,df/2) u1<- w1%*%q / g[1,1] # py<-p vv<-py/sum(py) cvv<-cumsum(vv) for (ijk in 1:Nrep) { r1<-hist(runif(nn),breaks=c(0,cvv),plot=F) q<-r1$count/nn g<- q%*%ac%*%q eg<- p%*%ac%*%p+sum(diag( ac%*%w ))/nn vg<-( 4* p%*%ac%*%w%*%ac%*%p + 2/nn*sum(diag( ac%*%w%*%ac%*%w )) )/nn skew<-8*( 3* p%*%ac%*%w%*%ac%*%w%*%ac%*%p + + (1/nn)*sum(diag(ac%*%w%*%ac%*%w%*%ac%*%w)) ) / sqrt(nn) / (nn*vg)**1.5 df<-8/skew**2 hh1<-(g-eg)/sqrt(vg) hh<-c(hh,hh1) } MCp<-rank(-hh)[1]/(Nrep+1) list(Stnd.C=tc, C.p.value=pval, C.u= u1, Freq=hh, MC.p.value=MCp) }