[R] fixed trimmed mean for group
Rui Barradas
ruipbarradas at sapo.pt
Sat Jul 7 15:46:24 CEST 2012
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 <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)
>
>
>
More information about the R-help
mailing list