R/SErob.R

In musto101/wilcox_R: What the package does (short line)

SErob <-
function(Z,mu,Sigma,theta,d,r,tau){
	n=nrow(Z)
	p=ncol(Z)
	ps=p*(p+1)/2
	q=6
	Dup=Dp(p)
	DupPlus=solve(t(Dup)%*%Dup)%*%t(Dup)

   	InvSigma=solve(Sigma)
	sigma=vech(Sigma)
   	W=0.5*t(Dup)%*%(InvSigma%x%InvSigma)%*%Dup
   	
   	Zr=matrix(NA,n,p) # robustly transformed data
   	A=matrix(0,p+q,p+q)
   	B=matrix(0,p+q,p+q)
 	sumg=rep(0,p+q)

   	for (i in 1:n) {
    	zi=Z[i,]
    	zi0=zi-mu
    	di=d[i]

     	if (di<=r) {
     			u1i=1.0
     			u2i=1.0/tau
     			du1i=0
     			du2i=0
		}else {
     			u1i=r/di
     			u2i=u1i^2/tau
     			du1i=-r/di^2
     			du2i=-2*r^2/tau/di^3
 		}

      	#robust transformed data
    	Zr[i,]=sqrt(u2i)*t(zi0)
    			   	
    	####	gi

    	g1i=u1i*zi0	# defined in (24)
      	vTi=vech(zi0%*%t(zi0))
    	g2i=u2i*vTi-sigma	# defined in (25)
	gi=rbind(g1i,g2i)
    	sumg=gi+sumg

		B=B+gi%*%t(gi)
    	
    	####	gdoti

    	#	derivatives of di with respect to mu and sigma
    	ddmu=-1/di*t(zi0)%*%InvSigma
    	ddsigma=-t(vTi)%*%W/di

    	#	
    	dg1imu=-u1i*diag(p)+du1i*zi0%*%ddmu
    	dg1isigma=du1i*zi0%*%ddsigma
    	dg2imu=-u2i*DupPlus%*%(zi0%x%diag(p)+diag(p)%x%zi0)+du2i*vTi%*%ddmu
    	dg2isigma=du2i*vTi%*%ddsigma-diag(q)
    	
    	dgi=rbind(cbind(dg1imu,dg1isigma),cbind(dg2imu,dg2isigma))
    	A=A+dgi
   } # end of loop n

	A=-1*A/n
	B=B/n
	invA=solve(A)
	OmegaSW=invA%*%B%*%t(invA)
   	OmegaSW=OmegaSW[(p+1):(p+q),(p+1):(p+q)]
	
		
   SEsw=getSE(theta,OmegaSW,n)
	SEinf=SEML(Zr,theta)$inf
	
	Results=list(inf=SEinf,sand=SEsw,Zr=Zr)  	
   return(Results)
   	
}