Hello,
Without trying it, I very much doubt that your code runs: First you use
'asim' and create it after.
So I've moved asim <- 1000 up some lines.
And it still doesn't work.
1. Run your own code before posting it.
2. Learn to use indentation.
3. Learn to use spaces to make your code easier to read.
Hope this helps,
Rui Barradas
P.S. At a glance, it seems you are dividing 'num' by J and only then
subtracting 1.
Missing parenthesis? Impossible to say. Follow the above first.
R.B.
Em 29-05-2012 07:06, umai88 escreveu:
Hello everyone, I want to calculate type I error rate for modified F
statistic for one way robust anova. I need to find the j group trimmed
mean and winsorized sum of squared deviations. Here I attached my code for
j=2 to make it simple. Originally I have j=4. Hope you can help. I need to
run it for 1000 times
My problem is:
i) the value of F-test obtain from my simulation below is in negative
value..There might be something wrong in my coding
ii) after obtain F value, how i can proceed to obtain type I error rate?
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=2
trim1<-rep(NA,asim)
trim2<-rep(NA,asim)
pw<-rep(NA,asim)
qw<-rep(NA,asim)
px<-rep(NA,asim)
qx<-rep(NA,asim)
ssd1<-rep(NA,asim)
ssd2<-rep(NA,asim)
madN1<-rep(NA,asim)
madN2<-rep(NA,asim)
lo.w<-rep(NA,asim)
up.w<-rep(NA,asim)
lo.x<-rep(NA,asim)
up.x<-rep(NA,asim)
hw<-rep(NA,asim)
hx<-rep(NA,asim)
xbar.t<-rep(NA,asim)
num<-rep(NA,asim)
df2<-rep(NA,asim)
denom<-rep(NA,asim)
F<-rep(NA,asim)
asim<-1000
for(j in 1:asim)
{
print(j)
set.seed(j)
a<-rnorm(15,0,1)
b<-rnorm(15,0,1)
w<-a*exp(h*a^2/2)
x<-b*exp(h*b^2/2)
matw<-sort(w)
matx<-sort(x)
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]))
madN1[j]<-mad(w)
madN2[j]<-mad(x)
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])))
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
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
hw[j]<-n-lo.w[j]-up.w[j]
hx[j]<-n-lo.x[j]-up.x[j]
H[j]=hw[j]+hx[j]
xbar.t[j]<-(hw[j]*trim1[j]+hx[j]*trim2[j])/H[j]
num[j]<-( (trim1[j]-xbar.t[j])^2 + (trim2[j]-xbar.t[j])^2)^2 )/J-1
df2[j]<-H[j]-J
denom[j]<-(ssd1[j]+ssd2[j])/df2[j]
F[j]<-num[j]/denom[j]
}
Graduate student,
Department of Mathematics& Statistics, UPM
--
View this message in context:
http://r.789695.n4.nabble.com/need-help-to-find-type-I-error-rate-for-modified-F-statistic-tp4631661.html
Sent from the R help mailing list archive at Nabble.com.
[[alternative HTML version deleted]]
______________________________________________
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.
______________________________________________
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.
