Hello,
I haven't found errors in your code. I implemented the test in the paper
(the first, fixed symetric mean) and it also gives me zero Type I
errors, when alpha = 0.05. Try to see the value of min(pv) or to plot
the histogram of 'pv', hist(pv) and you'll see that there are no
significant p-values, at that level.
Anyway I'll continue to look at it, but my first impression is that
there is something wrong somewhere and it doesn't seem to be in
following the formulae.
Rui Barradas
Em 07-07-2012 04:26, Agnes Ayang escreveu:
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 <ruipbarra...@sapo.pt
<mailto:ruipbarra...@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.ay...@gmail.com
<mailto:agnes.ay...@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)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
