Actually this r code is for my final year project paper..if this code
cannot give the correct answer which is between 0.025 and 0.075, i canot
obtain the result for my project paper. I have tried 3 times to sent
this programme to the r help but till now there is no answer whether the
r help already received my programe or not.
On Sat, Jul 7, 2012 at 12:34 AM, Rui Barradas <ruipbarradas at sapo.pt
<mailto:ruipbarradas at sapo.pt>> wrote:
Hello,
Why me? You haven't written to R-help, to one (or more?) R-helper, why?
I'm going to look at it with more attention today and during the
weekend because:
1. The paper you include a link to is very readable;
2. Your code works without changes, a very, very strong point to you;
3. It doesn't seem like homework, but if it is, you have made a
visible effort.
4. Finally, your code is also very readable, well organized.
Starting with this last point, some immediate notes.
1. Put spaces before and after =, +, *, etc.
2. To assign a value use '<-' not '=', the equal sign is reserved to
pass parameters to functions.
Anyway, I'll look at it, and I'm curious, why me?
Rui Barradas
Em 06-07-2012 08:41, agnes.ayang at gmail.com
<mailto:agnes.ayang at gmail.com> escreveu:
rDear Mr Rui Baradas,
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
<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 you 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)