ngrad<-function(beta,x,y,cx=FALSE,cr=TRUE) {
  x<-as.matrix(x)
  if (cx) x<-apply(x,2,function(x){x-mean(x)})
  y<-as.matrix(y)
  beta<-as.matrix(beta)
  e<-y-x%*%beta
  r<-rank(e)
  n<-dim(x)[1]
  if (cr)
    ans<-sum((r-(n+1)/2)*x)
  else
    ans<-sum(r*x)
  ans
}