[R] help

[Greg Snow]

You should probably read the posting guide, there is a link at the
bottom of every e-mail on the list.

You provide a large amount of code but no description of what you are
trying to do with it and very little about what problem you are
having.  Most people (well me at least) are unwilling to run and look
through that much code without some more idea of what we are supposed
to be looking for.  Can you reduce the amount of code to a smaller
example that shows your problem?  can you give us more detail about
what you expect the cod to do and how what you are seeing differs from
what you are expecting?

On Thu, May 29, 2014 at 3:36 AM, thanoon younis
<thanoon.younis80 at gmail.com> wrote:
> 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 at 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.



-- 
Gregory (Greg) L. Snow Ph.D.
538280 at gmail.com