[R] simulation of modified bartlett's test

[dila]

Hi, I run this code to get the power of the test for modified bartlett's
test..but I'm not really sure that my coding is right..

#normal distribution unequal variance
asim<-5000
pv<-rep(NA,asim)
for(i in 1:asim)
{print(i)
set.seed(i)
n1<-20
n2<-20
n3<-20
mu<-0
sd1<-sqrt(25)
sd2<-sqrt(50)
sd3<-sqrt(100)
g1<-rnorm(n1,mu,sd1)
g2<-rnorm(n2,mu,sd2)
g3<-rnorm(n3,mu,sd3)
x=c(g1,g2,g3)
group=c(rep(1,n1),rep(2,n2),rep(3,n3))
N=60
k=3
v1=var(g1)
v2=var(g2)
v3=var(g3)
#pooled variance
A=((n1-1)*v1+(n2-1)*v2+(n3-1)*v3)/(N-k)
#calculate B
B=((N-k)*(log(A)))-((n1-1)*log(v1)+(n2-1)*log(v2)+(n3-1)*log(v3))
#calculate C
C=1+(1/(3*(k-1))*(((1/(n1-1))+(1/(n2-1))+(1/(n3-1)))-(1/(N-k))))
#calculate layard estimator
xbar1=mean(g1)
xbar2=mean(g2)
xbar3=mean(g3)
sum1=sum((g1-xbar1)^4)
sum2=sum((g2-xbar2)^4)
sum3=sum((g3-xbar3)^4)
sum4=sum((g1-xbar1)^2)
sum5=sum((g2-xbar2)^2)
sum6=sum((g3-xbar3)^2)
y= (N*(sum1+sum2+sum3))/((sum4+sum5+sum6)^2)
#calculate bartlett modified statistic
bar2=B/(C*(1/2)*(y-1))
bar2
pv[i]<-pchisq(bar2,2,lower=FALSE)
}
mean(pv<0.01)
mean(pv<0.05)


--
View this message in context: http://r.789695.n4.nabble.com/simulation-of-modified-bartlett-s-test-tp4632305.html
Sent from the R help mailing list archive at Nabble.com.