[R] fixed trimmed mean for j-group

[umai88]

Hello...i want to find the empirical rate for type 1 error using fixed
trimmed mean.  To make it easy, i'm referring to journal given by this
website
http://www.academicjournals.org/ajmcsr/PDF/pdf2011/Yusof%20et%20al.pdf.
I already run the programme and there is no error in it but i got zero for
the empirical rate of type 1 error. The empirical rate for the type 1 error
given in the journal should lied between 0.025 and 0.07. Below is my
programme. Hope anyone can help me. Thank you.

## use for data transformation
g=0
h=0
  
  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)

g=0.5 and h=0
g=0.5 and h=0.5

w<-(exp(g*a)-1)/g*(exp((h*a^2)/2))  
x<-(exp(g*b)-1)/g*(exp((h*b^2)/2))
y<-(exp(g*c)-1)/g*(exp((h*c^2)/2))
z<-(exp(g*d)-1)/g*(exp((h*d^2)/2))

######################## FIXED SYMMETRIC TRIMMED MEAN 
###############################
asim<-5000
pv<-rep(NA, asim)
for(j in 1:asim)
{
print(j)
set.seed(j)
n1=15
n2=15
n3=15
n4=15
miu=0
sd1=1
sd2=1
sd3=1
sd4=1

a=rnorm(n1,miu,sd1)
b=rnorm(n2,miu,sd2)
c=rnorm(n3,miu,sd3)
d=rnorm(n4,miu,sd4)

## data transformation
g=0
h=0
  
  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)

mat1<-sort(w)
mat2<-sort(x)
mat3<-sort(y)
mat4<-sort(z)

alpha=0.15
k1=floor(alpha*n1)+1
k2=floor(alpha*n2)+1
k3=floor(alpha*n3)+1
k4=floor(alpha*n4)+1

r1=k1-(alpha*n1)
r2=k2-(alpha*n2)
r3=k3-(alpha*n3)
r4=k4-(alpha*n4)

## j-group trimmed mean

e1=k1+1
f1=n1-k1

e2=k2+1
f2=n2-k2

e3=k3+1
f3=n3-k3

e4=k4+1
f4=n4-k4

trim1=1/((1-2*alpha)*n1)*(sum(mat1[e1:f1]) + r1*(mat1[k1]+mat1[n1-k1+1]))
trim2=1/((1-2*alpha)*n2)*(sum(mat2[e2:f2]) + r2*(mat2[k2]+mat2[n2-k2+1]))
trim3=1/((1-2*alpha)*n3)*(sum(mat3[e3:f3]) + r3*(mat3[k3]+mat3[n3-k3+1]))
trim4=1/((1-2*alpha)*n4)*(sum(mat4[e4:f4]) + r4*(mat4[k4]+mat4[n4-k4+1]))

## sample winsorized mean
x1=(1-r1)*mat1[k1+1]+r1*mat1[k1]
x2=(1-r2)*mat2[k2+1]+r2*mat2[k2]
x3=(1-r3)*mat3[k3+1]+r3*mat3[k3]
x4=(1-r4)*mat4[k4+1]+r4*mat4[k4]

u1=(1-r1)*mat1[n1-k1]+r1*mat1[n1-k1+1]
u2=(1-r2)*mat2[n2-k2]+r2*mat2[n2-k2+1]
u3=(1-r3)*mat3[n3-k3]+r3*mat3[n3-k3+1]
u4=(1-r4)*mat4[n4-k4]+r4*mat4[n4-k4+1]

win1=1/n1* (sum(mat1[e1:f1])+k1*(x1+u1))
win2=1/n2* (sum(mat2[e2:f2])+k2*(x2+u2))
win3=1/n3* (sum(mat3[e3:f3])+k3*(x3+u3))
win4=1/n4* (sum(mat4[e4:f4])+k4*(x4+u4))

## g-winsorized sum of squared deviations

ssd1=sum((mat1[e1:f1]-win1)^2) + k1*((mat1[k1+1]-win1)^2 +
(mat1[n1-k1]-win1)^2)
ssd2=sum((mat2[e2:f2]-win2)^2) + k2*((mat2[k2+1]-win2)^2 +
(mat2[n2-k2]-win2)^2)
ssd3=sum((mat3[e3:f3]-win3)^2) + k3*((mat3[k3+1]-win3)^2 +
(mat3[n3-k3]-win3)^2)
ssd4=sum((mat4[e4:f4]-win4)^2) + k4*((mat4[k4+1]-win4)^2 +
(mat4[n4-k4]-win4)^2)

## calculate f statistic

J=4
h1=n1-2*floor(alpha*n1)
h2=n2-2*floor(alpha*n2)
h3=n3-2*floor(alpha*n3)
h4=n4-2*floor(alpha*n4)

H=h1+h2+h3+h4

xt=(h1*trim1+h2*trim2+h3*trim3+h4*trim4)/H

num= ((trim1-xt)^2+(trim2-xt)^2+(trim3-xt)^2+(trim4-xt)^2)/(J-1)
denom= (ssd1+ssd2+ssd3+ssd4)/(H-J)

ft=num/denom
pv[j]=1-pf(ft,(J-1),(H-J))
}
mean(pv<0.05)

##################################################################################

Graduate student of Universiti Putra Malaysia

--
View this message in context: http://r.789695.n4.nabble.com/fixed-trimmed-mean-for-j-group-tp4633327.html
Sent from the R help mailing list archive at Nabble.com.