Nothing
#-------------------------------------------------------
#
# Created : 29/10/02
# Last Modified : Time-stamp: <2018-05-11 20:49:49 (antoine)>
#
# Description : Principal component analysis
#
# Author : Antoine Lucas
# lucas@toulouse.inra.fr
#
# Licence : GPL
#
#-------------------------------------------------------
acp <- function(x,center=TRUE,reduce=TRUE,wI=rep(1,nrow(x)),wV=rep(1,ncol(x)))
{
x <- as.matrix(x)
if(center)
x <- t(t(x) - as.vector(( wI %*% x)/sum(wI)))
## x <- scale(x ,center = center, scale = FALSE)
if (reduce)
x <- apply(x,2,function(u) { u/stats::sd(u)})
## Di.X'.Dv.X
EIG <- eigen( (t(x)* wI) %*% (x * wV) ,symmetric=FALSE)
V <- EIG$vector # ou bien: V=svd(x)$v
EIG$values <- Re(EIG$values)
V <- V %*% diag(sign(EIG$values))
val <- sqrt(abs(EIG$values))
scores <- x %*% V
V <- as.matrix(Re(V))
scores <- as.matrix(Re(scores))
dimnames(V)[[2]] <- paste("Comp",1:dim(x)[2])
if(!is.null( dimnames(x)[[2]] ))
dimnames(V)[[1]] <- dimnames(x)[[2]]
if(!is.null(dimnames(x)[[1]]))
dimnames(scores)[[1]] <- dimnames(x)[[1]]
dimnames(scores)[[2]] <- paste("Comp",1:dim(x)[2])
##cmpr <- x %*% (sqrt(wV) * as.matrix(V))
sdev <- apply(scores,2,stats::sd)
res <- list(eig=val,sdev=sdev,scores=scores,loadings=V)
class(res) <- "acp"
res
}
pca <- acp
print.acp <- function(x, ...)
{
#cat("Call:\n"); dput(x$call)
cat("\nStandard deviations:\n")
print(x$sdev, ...)
cat("\nEigen values:\n")
print(x$eig, ...)
invisible(x)
}
#
# SECTION GRAPHIQUES
#
plot.acp <- function(x,i=1,j=2,text=TRUE,label='Composants',col='darkblue',main='Individuals PCA',variables=TRUE,individual.label=NULL,...)
{
U <- x$scores
XLAB <- paste(label,i)
YLAB <- paste(label,j)
plot.new()
plot.window(range(U[,i]),range(U[,j]))
axis(1,labels=TRUE,tick=TRUE)
axis(2,labels=TRUE,tick=TRUE)
box()
title(xlab=XLAB,ylab=YLAB,main=main)
if(text){
if(is.null(individual.label))
{
individual.label=dimnames(x$scores)[[1]]
}
text(labels=individual.label,U[,i],U[,j],col=col,...)
}
else{
points(U[,i],U[,j],col=col,...)
}
if(variables)
{
par(new=TRUE)
biplot.acp(x,circle=FALSE,label="",main="")
}
}
biplot.acp <- function(x,i=1,j=2,label='Composants',col='darkblue',length=0.1,main='Variables PCA',circle=TRUE,...)
{
U <- x$loadings
LIM <- c(-1.3,1.3)
XLAB <- paste(label,i)
YLAB <- paste(label,j)
# PLOT DES AXES
graphics::plot.new()
graphics::plot.window(LIM,LIM)
graphics::axis(1,labels=TRUE,tick=TRUE)
graphics::axis(2,labels=TRUE,tick=TRUE)
graphics::box()
graphics::title(xlab=XLAB,ylab=YLAB,main=main)
# PLOT DU NOM DES FLECHES
graphics::text(x=U[,i]*1.3,y=U[,j]*1.3,labels=dimnames(U)[[1]],col=col)
# PLOT DES FLECHES
graphics::arrows(0,0,U[,i],U[,j],length = length,col=col)
# CERCLE
if(circle)
{
t2p <- 2 * pi * seq(0,1, length = 200)
xc <- cos(t2p)
yc <- sin(t2p)
graphics::lines(xc,yc,col='darkblue')
}
}
# Graphique: Eboulis des valeurs propres
plot2 <- function(x,pourcent=FALSE,eigen=TRUE,label='Comp.',col='lightgrey',main='Scree Graph',ylab='Eigen Values')
{
if(eigen){ U <- x$eig }
else { U <- x$sdev }
if(pourcent){U <- U/sum(U) }
n <- length(U)
names <- paste(label,1:n)
graphics::barplot(U,main=main,ylab=ylab,col=col,names.arg=names)
}
plotAll <- function(x)
{
graphics::par(mfrow=c(2,2))
plot2(x)
## boxplot(as.list(as.data.frame(x$cmpr)))
graphics::plot(x,variables=FALSE)
stats::biplot(x)
graphics::plot(x,main="Both",variables=TRUE)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.