i want to find type I error rates using modified F statistic..below is my R
code..
h=0
g=0
n=15
alpha=0.1
k=floor(alpha*n)+1
r=k-(alpha*n)
i=k+1
m=n-k
J=4
trim1<-rep(NA,asim)
trim2<-rep(NA,asim)
trim3<-rep(NA,asim)
trim4<-rep(NA,asim)
pw<-rep(NA,asim)
qw<-rep(NA,asim)
px<-rep(NA,asim)
qx<-rep(NA,asim)
py<-rep(NA,asim)
qy<-rep(NA,asim)
pz<-rep(NA,asim)
qz<-rep(NA,asim)
ssd1<-rep(NA,asim)
ssd2<-rep(NA,asim)
ssd3<-rep(NA,asim)
ssd4<-rep(NA,asim)
madN1<-rep(NA,asim)
madN2<-rep(NA,asim)
madN3<-rep(NA,asim)
madN4<-rep(NA,asim)
lo.w<-rep(NA,asim)
up.w<-rep(NA,asim)
lo.x<-rep(NA,asim)
up.x<-rep(NA,asim)
lo.y<-rep(NA,asim)
up.y<-rep(NA,asim)
lo.z<-rep(NA,asim)
up.z<-rep(NA,asim)
hw<-rep(NA,asim)
hx<-rep(NA,asim)
hy<-rep(NA,asim)
hz<-rep(NA,asim)
H<-rep(NA,asim)
xbar.t<-rep(NA,asim)
num<-rep(NA,asim)
df2<-rep(NA,asim)
denom<-rep(NA,asim)
F<-rep(NA,asim)
F.table<-rep(NA,asim)
asim<-10
for(j in 1:asim)
{
print(j)
set.seed(j)
a<-rnorm(15,0,1)
b<-rnorm(15,0,1)
c<-rnorm(15,0,1)
d<-rnorm(15,0,1)
w<-a*exp(h*a^2/2)
x<-b*exp(h*b^2/2)
y<-c*exp(h*c^2/2)
z<-d*exp(h*d^2/2)
matw<-sort(w)
matx<-sort(x)
maty<-sort(y)
matz<-sort(z)
trim1[j]=1/((1-2*alpha)*n)*(sum(matw[i:m])+r*(matw[k]+matw[n-k+1]))
trim2[j]=1/((1-2*alpha)*n)*(sum(matx[i:m])+r*(matx[k]+matx[n-k+1]))
trim3[j]=1/((1-2*alpha)*n)*(sum(maty[i:m])+r*(maty[k]+maty[n-k+1]))
trim4[j]=1/((1-2*alpha)*n)*(sum(matz[i:m])+r*(matz[k]+matz[n-k+1]))
madN1[j]<-mad(w)
madN2[j]<-mad(x)
madN3[j]<-mad(y)
madN4[j]<-mad(z)
lo.w[j]<-length(which(w-median(w)<(-2.24*madN1[j])))
up.w[j]<-length(which(w-median(w)>(2.24*madN1[j])))
lo.x[j]<-length(which(x-median(x)<(-2.24*madN2[j])))
up.x[j]<-length(which(x-median(x)>(2.24*madN2[j])))
lo.y[j]<-length(which(y-median(y)<(-2.24*madN3[j])))
up.y[j]<-length(which(y-median(y)>(2.24*madN3[j])))
lo.z[j]<-length(which(z-median(z)<(-2.24*madN4[j])))
up.z[j]<-length(which(z-median(z)>(2.24*madN4[j])))
pw[j]<-up.w[j]+2
qw[j]<-n-up.w[j]-1
px[j]<-up.x[j]+2
qx[j]<-n-up.x[j]-1
py[j]<-up.y[j]+2
qy[j]<-n-up.y[j]-1
pz[j]<-up.z[j]+2
qz[j]<-n-up.z[j]-1
ssd1[j]<-lo.w[j]*(matw[lo.w[j]+1]-trim1[j])^2 +
(matw[pw[j]:qw[j]]-trim1[j])^2 + up.w[j]*(matw[n-up.w]-trim1[j])^2 -
((lo.w)*(matw[lo.w+1]-trim1[j])+ (up.w)*(matw[n-up.w]-trim1[j]))^2
ssd2[j]<-lo.x[j]*(matx[lo.x[j]+1]-trim2[j])^2 +
(matx[px[j]:qx[j]]-trim2[j])^2 + up.x[j]*(matx[n-up.x]-trim2[j])^2 -
((lo.x)*(matw[lo.x+1]-trim2[j])+ (up.x)*(matx[n-up.x]-trim2[j]))^2
ssd3[j]<-lo.y[j]*(maty[lo.y[j]+1]-trim3[j])^2 +
(maty[py[j]:qy[j]]-trim3[j])^2 + up.y[j]*(maty[n-up.y]-trim3[j])^2 -
((lo.y)*(matw[lo.y+1]-trim3[j])+ (up.y)*(maty[n-up.y]-trim3[j]))^2
ssd4[j]<-lo.z[j]*(matz[lo.z[j]+1]-trim4[j])^2 +
(matz[pz[j]:qz[j]]-trim4[j])^2 + up.z[j]*(matz[n-up.z]-trim4[j])^2 -
((lo.z)*(matw[lo.z+1]-trim4[j])+ (up.z)*(matz[n-up.z]-trim4[j]))^2
hw[j]<-n-lo.w[j]-up.w[j]
hx[j]<-n-lo.x[j]-up.x[j]
hy[j]<-n-lo.y[j]-up.y[j]
hz[j]<-n-lo.z[j]-up.z[j]
H[j]=hw[j]+hx[j]+hy[j]+hz[j]
xbar.t[j]<-(hw[j]*trim1[j]+hx[j]*trim2[j]+hy[j]*trim3[j]+hz[j]*trim4[j])/H[j]
num[j]<-( (trim1[j]-xbar.t[j])^2 + (trim2[j]-xbar.t[j])^2 +
(trim3[j]-xbar.t[j])^2+ (trim4[j]-xbar.t[j])^2 )/J-1
df2[j]<-H[j]-J
denom[j]<-(ssd1[j]+ssd2[j]+ssd3[j]+ssd4[j])/df2[j]
F[j]<-num[j]/denom[j]
}
--
View this message in context:
http://r.789695.n4.nabble.com/obtain-type-I-error-from-modified-F-statistic-tp4631660.html
Sent from the R help mailing list archive at Nabble.com.