dear all members
i have some problems in R- code below, is there anyone helps me to correct
it.
many thanks in advance
thanoon
################## Ridge MM ##################
rm(list=ls())
library(MASS)
library(mvoutlier)
library(robustbase)
library(car)
library(quantreg)
n<-100
r<-0.95
R<-10
p<-3
B<-matrix(c(1,1,1,1),p+1,1)
n.out<-round(0.10*n)
set.seed(n)
b0=matrix(nrow=R,ncol=1)
b1=matrix(nrow=R,ncol=1)
b2=matrix(nrow=R,ncol=1)
b3=matrix(nrow=R,ncol=1)
b.MM=matrix(nrow=R,ncol=p+1)
b0.final=matrix(nrow=R,ncol=1)
b1.final=matrix(nrow=R,ncol=1)
b2.final=matrix(nrow=R,ncol=1)
b3.final=matrix(nrow=R,ncol=1)
MMR.BS=matrix(nrow=p+1,ncol=R)
WL=matrix(nrow=p+1,ncol=R)
var=matrix(nrow=p+1,ncol=R)
beta.error=matrix(nrow=p+1,ncol=R)
for (i in 1:R){
z11<-rnorm(n,0,1)
z22<-rnorm(n,0,1)
z33<-rnorm(n,0,1)
z44<-rnorm(n,0,1)
x1<-(1-r)*z11+sqrt(r)*z44
x2<-(1-r)*z22+sqrt(r)*z44
x3<-(1-r)*z33+sqrt(r)*z44
x<-as.matrix(data.frame(x1,x2,x3))
e<-rnorm(n,0,1)
e[1:n.out]<-100
y<-1+x1+x2+x3+e
stand<-function(x){(x-mean(x))/sd(x)*(n-1)^-0.5}
x11<-stand(x1)
x22<-stand(x2)
x33<-stand(x3)
X<-as.matrix(data.frame(1,x11,x22,x33))
y1<-stand(y)
data1<-data.frame(y1,x11,x22,x33)
b.MM<- rlm(y1~.,data1,psi=psi.bisquare,method="MM",k2=4.685)
beta=as.matrix(coef(b.MM))
b0[i]=beta[1]
b1[i]=beta[2]
b2[i]=beta[3]
b3[i]=beta[4]
eMM<-rlm(y1~.,data1,psi=psi.bisquare,method="MM",k2=4.685)$residuals
sMM<-sum(eMM^2)/(n-p-1)
b.MM.final<- rlm(y1~.,data1,psi=psi.bisquare,method="MM",k2=4.685)$coef
kMM<-4*sMM^2/t(b.MM.final)%*%b.MM.final
kMM<-c(kMM)
A<-t(X)%*%X
II<-diag(x = 1, nrow=p+1, ncol=p+1)
KMM<-II*kMM
MMR.BS[,i]<-solve(t(X)%*%X+KMM)%*%t(X)%*%X%*%b.MM.final
beta.error[,i]=(MMR.BS[,i]-c(1,1,1,1))^2
WL[i,]<-solve(II+kMM*solve(A))
MMR.var[i,]<-sMM*WL[,i]%*%solve(A)%*%t(WL[,i])
}
##############################################
#Bias,RMSE and SE for MMR
MMR.AVE.b1<-mean(MMR.BS[2,])
MMR.AVE.b2<-mean(MMR.BS[3,])
MMR.AVE.b3<-mean(MMR.BS[4,])
bias
MMR.bi.b1<-mean(MMR.BS[2,])-1
MMR.bi.b2<-mean(MMR.BS[3,])-1
MMR.bi.b3<-mean(MMR.BS[4,])-1
MSE
MSE.b1<- (beta.error[2,])+(MMR.var[2,])
MSE.b2<- (beta.error[3,])+(MMR.var[3,])
MSE.b3<- (beta.error[4,])+(MMR.var[4,])
RMSE
RMSE.b1<-sqrt((beta.error[2,])+(MMR.var[2,]))
RMSE.b2<-sqrt((beta.error[3,])+(MMR.var[3,]))
RMSE.b3<-sqrt((beta.error[4,])+(MMR.var[4,]))
sdd
sd.b1=sqrt(MSE.b1-(MMR.bi.b1)^2)
sd.b2=sqrt(MSE.b2-(MMR.bi.b2)^2)
sd.b3=sqrt(MSE.b3-(MMR.bi.b3)^2)
results
res.b1=cbind(MMR.AVE.b1,MMR.bi.b1,MSE.b1,RMSE.b1,sd.b1)
res.b2=cbind(MMR.AVE.b2,MMR.bi.b2,MSE.b2,RMSE.b2,sd.b2)
res.b3=cbind(MMR.AVE.b3,MMR.bi.b3,MSE.b3,RMSE.b3,sd.b3)
res.b1
res.b2
res.b3
[[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.
