R/sonscabasic.R

sonscabasic <-
function (Xtable,mj) 
{
x <- Xtable/sum(Xtable)
rsums <- as.matrix(rowSums(x))
csums <- as.vector(colSums(x))
tauden <- 1 - sum(rsums^2)
drm1 <- diag( 1/( rsums + (rsums==0) ) * (1-(rsums==0)) )
dcm1 <- diag( 1/( csums + (csums==0) ) * (1-(csums==0)) )
drmh<-diag(rep(1,nrow(x))) #change the metric in NSCA
dcmh <- sqrt(dcm1)
dj <- diag(csums)
di <- diag(rsums)
uni <- matrix(1, 1, ncol(x))
uni1 <- rep(1, nrow(x))
Bpoly <- emerson.poly(mj, csums)$B
Bpoly2 <- sqrt(dj) %*% Bpoly
#pcc <- 1/sqrt(tauden)* ( x%*%dcm1 - rsums %*% (uni) ) 
pcc <- ( x%*%dcm1 - rsums %*% (uni) ) 
svdpccw<-svd(pcc%*%sqrt(dj))
u<-svdpccw$u
#mu<-svdpccw$d
#Z <- t(u)  %*%pcc %*% sqrt(dj)%*%Bpoly2
Z <- t(u)  %*%pcc %*% dj%*%Bpoly
ZtZ<-Z%*%t(Z)
tZZ<-t(Z)%*%Z
mu2<- diag(tZZ) #only the sum gives me the total inertia
mu<- diag(ZtZ) #these are coincident with each eigenvalue (mu^2)
#tau<-sum(mu)
#r<-rmax
#browser()
#sonsca <- new("cabasicresults",
 #         RX=pcc,CX=t(pcc),Rweights=dj,Cweights=diag(uni1),
  #        Raxes=Bpoly,Caxes=u,mu=mu,mu2=mu2,catype="SONSCA",tauDen=tauden,Z=Z,ZtZ=ZtZ,tZZ=tZZ)

ressonsca<-( list(RX=pcc,CX=t(pcc),Rweights=dj,Cweights=diag(uni1),
          Raxes=Bpoly,Caxes=u,mu=mu,mu2=mu2,tauDen=tauden,catype="SONSCA",Z=Z,ZtZ=ZtZ,tZZ=tZZ))

return(ressonsca)
}

Try the CAvariants package in your browser

Any scripts or data that you put into this service are public.

CAvariants documentation built on Oct. 20, 2023, 1:07 a.m.