Nothing
R/orthobasis.R
In adegenet: Exploratory Analysis of Genetic and Genomic Data
Defines functions
.orthobasis.listw
## THIS CODE HAS BEEN TAKEN FROM ADE4
## IT HAS BEEN MOVED TO ADESPATIAL
## THIS IS MERELY A TEMPORARY FIX.
## only orthobasis.listw is used.
## renamed into .orthobasis.listw to avoid having to document it.
## - Thibaut Jombart, April 2015
#######################################################
.orthobasis.listw <- function( listw) {
appel = match.call()
if(!inherits(listw,"listw")) stop ("object of class 'listw' expected")
if(listw$style!="W") stop ("object of class 'listw' with style 'W' expected")
n = length(listw$weights)
fun <- function (x) {
num = listw$neighbours[[x]]
wei = listw$weights[[x]]
res = rep(0,n)
res[num] = wei
return (res)
}
b0 <- matrix(unlist(lapply(1:n,fun)),n,n)
b0=(t(b0)+b0)/2
b0=bicenter.wt(b0)
a0 <- eigen(b0, symmetric = TRUE)
#barplot(a0$values)
a0 <- a0$vectors
a0 <- cbind(rep(1,n),a0)
a0 <- qr.Q(qr(a0))
a0 <- as.data.frame(a0[,-1])*sqrt(n)
row.names(a0) <- attr(listw,"region.id")
names(a0) <- paste("VP", 1:(n-1), sep = "")
z <- apply(a0,2,function(x) sum((t(b0*x)*x))/n)
attr(a0,"values") <- z
attr(a0,"weights") <- rep(1/n,n)
attr(a0,"call") <- appel
attr(a0,"class") <- c("orthobasis","data.frame")
return(a0)
}
## print.orthobasis <- function(x,...) {
## if (!inherits(x,"orthobasis")) stop ("for 'orthobasis' object")
## cat("Orthonormal basis: ")
## n <- nrow(x)
## p <- ncol(x)
## if (n!=(p+1)) stop ("Non convenient dimension: author's error")
## cat("data.frame with",n,"rows and",ncol(x),"columns\n")
## cat("--------------------------------------\n")
## cat("Columns are an orthonormal basis of 1n-orthogonal for\n")
## cat("the inner product defined by the weights attribute\n")
## cat("---------------------------------------\n")
## w <- attributes(x)
## if (!is.null(w$"names")) cat("names =", w$names[1],"...",w$names[p],"\n")
## if (!is.null(w$"row.names")) cat("row.names =", w$row.names[1],"...",w$row.names[n],"\n")
## if (!is.null(w$"weights")) cat("weights =", w$weights[1],"...",w$weights[n],"\n")
## if (!is.null(w$"values")) cat("values =", w$values[1],"...",w$values[p],"\n")
## if (!is.null(w$"class")) cat("class =", w$class,"\n")
## if (!is.null(w$"call")) {
## cat("call =")
## print(w$"call")
## }
## }
## orthobasis.mat <- function(mat, cnw=TRUE) {
## if (!is.matrix(mat)) stop ("matrix expected")
## if (any(mat<0)) stop ("negative value in 'mat'")
## if (nrow(mat)!=ncol(mat)) stop ("squared matrix expected")
## mat <- (mat+t(mat))/2
## nlig <- nrow(mat)
## if (is.null(dimnames(mat))) {
## w <- paste("P",1:nrow(mat),sep="")
## dimnames(mat) <- list(w,w)
## }
## labels <- dimnames(mat)[[1]]
## if (cnw) {
## margi <- apply(mat,1,sum)
## margi <- max(margi)-margi
## mat <- mat+diag(margi)
## }
## mat <- mat/sum(mat)
## wt <- rep ((1/nlig),nlig)
## # calculs extensibles a une ponderation quelconque
## wt <- wt/sum(wt)
## # si mat wt est la ponderation marginale associee a mat
## # tot = sum(mat)
## # mat = mat-matrix(wt,nlig,nlig,byrow=TRUE)*wt*tot
## # encore plus particulier mat = mat-1/nlig/nlig
## # en general les precedents sont des cas particuliers
## U <- matrix(1,nlig,nlig)
## U <- diag(1,nlig)-U*wt
## mat <- U%*%mat%*%t(U)
## wt <- sqrt(wt)
## mat <- t(t(mat)/wt)
## mat <- mat/wt
## eig <- eigen(mat,sym=TRUE)
## w0 <- abs(eig$values)/max(abs(eig$values))
## tol <- 1e-07
## w0 <- which(w0<tol)
## if (length(w0)==0) stop ("abnormal output : no null eigenvalue")
## else if (length(w0)==1) w0 <- (1:nlig)[-w0]
## else if (length(w0)>1) {
## # on ajoute le vecteur derive de 1n
## w <- cbind(wt,eig$vectors[,w0])
## # on orthonormalise l'ensemble
## w <- qr.Q(qr(w))
## # on met les valeurs propres a 0
## eig$values[w0] <- 0
## # on remplace les vecteurs du noyau par une base orthonormee contenant
## # en premiere position le parasite
## eig$vectors[,w0] <- w[,-ncol(w)]
## # on enleve la position du parasite
## w0 <- (1:nlig)[-w0[1]]
## }
## mat <- eig$vectors[,w0]/wt
## mat <- data.frame(mat)
## row.names(mat) <- labels
## names(mat) <- paste("S",1:(nlig-1),sep="")
## attr(mat,"values") <- eig$values[w0]
## attr(mat,"weights") <- rep(1/nlig,nlig)
## attr(mat,"call") <- match.call()
## attr(mat,"class") <- c("orthobasis","data.frame")
## return(mat)
## }
## "orthobasis.haar" <- function(n) {
## # on definit deux fonctions :
## appel = match.call()
## a <- log(n)/log(2)
## b <- floor(a)
## if ((a-b)^2>1e-10) stop ("Haar is not a power of 2")
## # la premiere est ecrite par Daniel et elle donne la demonstration (par analogie avec la fonction qui construit la base Bscores)
## # que la base Bscores est exactement la base de Haar quand on prend une phylogenie reguliere resolue.
## "haar.basis.1" <- function (n) {
## pari <- matrix(c(1,n),1)
## "div2" <- function (mat) {
## res <- NULL
## for (k in 1 : nrow(mat)) {
## n1 <- mat[k,1]
## n2 <- mat[k,2]
## diff <- n2-n1
## if (diff <=0) break
## n3 <- floor((n1+n2)/2)
## res <- rbind(res,c(n1,n3),c(n3+1,n2))
## }
## if (!is.null(res)) pari <<- rbind(pari,res)
## return(res)
## }
## mat <- div2(pari)
## while (!is.null(mat)) mat <- div2(mat)
## res <- NULL
## for (k in 1:nrow(pari)) {
## x<-rep(0,n)
## x[(pari[k,1]):(pari[k,2])] <- 1
## res <-c(res,x)
## }
## res = matrix(res,n)
## res <- qr.Q(qr(res))
## res <- res[, -1] * sqrt(n)
## res <- data.frame(res)
## row.names(res) <- paste("u",1:n,sep="")
## names(res) <- paste("B",1:(n-1),sep="")
## return(res)
## }
## # la seconde exploite les potentialites de la librairie waveslim, en remarquant qu'il existe un lien etroit entre la definition des filtres et la definition
## # des bases. Cette strategie permettra a l'avenir de definir les bases associees a d'autres famille de fonctions.
## "haar.basis.2" <- function (n) {
## if (!require(waveslim)) stop ("Please install waveslim")
## J <- a #nombre de niveau
## res <- matrix(0, nrow = n,ncol = n-1)
## filter.seq <- "H" #filtre correspondant au niveau 1
## h <- waveslim::wavelet.filter(wf.name = "haar", filter.seq = filter.seq) #parametre du filtre au niveau 1
## k <- 0
## for(i in 1:J){
## z <- rep(h,2**(J-i))
## x <- 1:n
## y <- rep((n-1-k):(n-2**(J-i)-k),rep(2**i,2**(J-i)))
## for(j in 1:n) res[x[j],y[j]] <- z[j]
## k <- k+2**(J-i)
## filter.seq <- paste(filter.seq, "L", sep = "")
## h <- waveslim::wavelet.filter(wf.name = "haar", filter.seq = filter.seq)
## }
## res <- res*sqrt(n)
## res <- data.frame(res)
## row.names(res) <- paste("u", 1:n, sep = "")
## names(res) <- paste("B", 1:(n-1), sep = "")
## return(res)
## }
## # suivant que n est grand (n > 257) ou non, on choisit l'une des deux strategies :
## if (n < 257)
## res <- haar.basis.1(n)
## else
## res <- haar.basis.2(n)
## attr(res,"values") <- NULL
## attr(res,"weights") <- rep(1/n,n)
## attr(res,"call") <- appel
## attr(res,"class") <- c("orthobasis","data.frame")
## return(res)
## }
## "orthobasis.line" <- function (n) {
## appel <- match.call()
## # solution de Cornillon p. 12
## res <- NULL
## for (k in 1:(n-1)) {
## x <- cos(k*pi*(2*(1:n)-1)/2/n)
## x <- sqrt(n)*x/sqrt(sum(x*x))
## res <-c(res,x)
## }
## res=matrix(res,n)
## res <- data.frame(res)
## row.names(res) <- paste("u",1:n,sep="")
## names(res) <- paste("B",1:(n-1),sep="")
## w <- (1:(n-1))*pi/2/n
## valpro <- 4*(sin(w)^2)/n
## poivoisi <- c(1,rep(2,n-2),1)
## poivoisi <- poivoisi/sum(poivoisi)
## norm <- unlist(apply(res, 2, function(a) sum(a*a*poivoisi)))
## y <- valpro*n*n/2/(n-1)
## val <- norm - y
## attr(res,"values") <- val
## attr(res,"weights") <- rep(1/n,n)
## attr(res,"call") <- appel
## attr(res,"class") <- c("orthobasis","data.frame")
## # verification locale. Ce paragraphe verifie que les vecteurs et les valeurs
## # proposee par Cornillon p. 12 sont bien les vecteurs propres de l'operateur de voisinage
## # rangee dans la solution analytique par variance locale croissante
## # l'article de Meot est errone et a donne le graphe circulaire pour le graphe lineaire
## # d0=neig2mat(neig(n.lin=n))
## # d0 = d0/n
## # d1=apply(d0,1,sum)
## # d0=diag(d1)-d0
## # fun2 <- function(x) {
## # z <- sum(t(d0*x)*x)/n
## # z <- z/sum(x*x)
## # return(z)
## # }
## # lambda <- unlist(apply(res,2,fun2))
## # print(lambda)
## # print(attr(res,"values"))
## # plot(lambda,attr(res,"values"))
## # abline(lm(attr(res,"values")~lambda))
## # print(coefficients(lm(attr(res,"values")~lambda)))
## # verification que les valeurs derivees des valeurs propres sont exactement des indices de Moran
## # d = neig2mat(neig(n.lin=n))
## # d = d/sum(d) # Moran type W
## # moran <- unlist(lapply(res,function(x) sum(t(d*x)*x)))
## # print(moran)
## # plot(moran,attr(res,"values"))
## # abline(lm(attr(res,"values")~moran))
## # print(summary(lm(attr(res,"values")~moran)))
## return(res)
## }
## "orthobasis.circ" <- function (n) {
## appel = match.call()
## if (n<3) stop ("'n' too small")
## "vecprosin" <- function(k) {
## x <- sin(2*k*pi*(1:n)/n)
## x <- x/sqrt(sum(x*x))
## }
## "vecprocos" <- function(k) {
## x <- cos(2*k*pi*(1:n)/n)
## x <- x/sqrt(sum(x*x))
## }
## "valpro" <- function(k,bis=TRUE) {
## x <- (4/n)*((sin(k*pi/n))^2)
## if (bis) x <- c(x,x)
## return(x)
## }
## k <- floor(n/2)
## if (k==n/2) {
## #n est pair
## w1 <- matrix(unlist(lapply(1:k,vecprocos)),n,k)
## w2 <- matrix(unlist(lapply(1:(k-1),vecprosin)),n,k-1)
## res <- cbind(w1,w2)
## res[,seq(1,2*k-1,by=2)]<-w1
## res[,seq(2,2*k-2,by=2)]<-w2
## vp <- unlist(lapply(1:(k-1),valpro))
## vp <- c(vp, valpro(k,FALSE))
## } else {
## # n est impair
## w1 <- matrix(unlist(lapply(1:k,vecprocos)),n,k)
## w2 <- matrix(unlist(lapply(1:k,vecprosin)),n,k)
## res <- cbind(w1,w2)
## res[,seq(1,2*k-1,by=2)]<-w1
## res[,seq(2,2*k,by=2)]<-w2
## vp <- unlist(lapply(1:k,valpro))
## }
## res=sqrt(n)*res
## res <- as.data.frame(res)
## row.names(res) <- paste("u",1:n,sep="")
## names(res) <- paste("B",1:(n-1),sep="")
## attr(res,"values") <- 1 - n*vp/2
## attr(res,"weights") <- rep(1/n,n)
## attr(res,"call") <- appel
## attr(res,"class") <- c("orthobasis","data.frame")
## # verification qu'on a exactement des indices de Moran a partie des valeurs propres
## # d = neig2mat(neig(n.cir=n))
## # d = d/sum(d) # Moran type W
## # moran <- unlist(lapply(res,function(x) sum(t(d*x)*x)))
## # print(moran)
## # plot(moran,attr(res,"values"))
## # abline(lm(attr(res,"values")~moran))
## # print(summary(lm(attr(res,"values")~moran)))
## return(res)
## }
## "orthobasis.neig" <- function( neig) {
## appel = match.call()
## if(!inherits(neig,"neig")) stop ("object of class 'neig' expected")
## n <- length(attr(neig,"degree"))
## m <- sum(attr(neig,"degree"))
## poivoisi <- attr(neig,"degree")/m
## if (is.null(names(poivoisi))) names(poivoisi) <- as.character(1:n)
## d0 = neig2mat(neig)
## d0 = diag(poivoisi)-d0/m
## eig <- eigen(d0, sym = TRUE)
## ########
## tol <- 1e-07
## w0 <- abs(eig$values)/max(abs(eig$values))
## w0 <- which(w0<tol)
## if (length(w0)==0) stop ("abnormal output : no null eigenvalue")
## else if (length(w0)==1) w0 <- (1:n)[-w0]
## else if (length(w0)>1) {
## # on ajoute le vecteur derive de 1n
## wt <- rep(1,n)
## w <- cbind(wt,eig$vectors[,w0])
## # on orthonormalise l'ensemble
## w <- qr.Q(qr(w))
## # on met les valeurs propres a 0
## eig$values[w0] <- 0
## # on remplace les vecteurs du noyau par une base orthonormee contenant
## # en premiere position le parasite
## eig$vectors[,w0] <- w[,-ncol(w)]
## # on enleve la position du parasite
## w0 <- (1:n)[-w0[1]]
## }
## w0 <- rev(w0)
## valpro <- eig$values[w0]
## eig <- eig$vectors[,w0]
## eig <- as.data.frame(eig)*sqrt(n)
## z <- apply(eig,2,function(x) sum(x*x*poivoisi))
## z <- z - valpro*n
## w <- rev(order(z))
## z <- z[w]
## eig <- eig[,w]
## row.names(eig) <- names(poivoisi)
## names(eig) <- paste("VP", 1:(n-1), sep = "")
## attr(eig,"values") <- z
## attr(eig,"weights") <- rep(1/n,n)
## attr(eig,"call") <- appel
## attr(eig,"class") <- c("orthobasis","data.frame")
## return(eig)
## }
Try the adegenet package in your browser
Any scripts or data that you put into this service are public.
adegenet documentation built on Feb. 16, 2023, 6 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .orthobasis.listw
## THIS CODE HAS BEEN TAKEN FROM ADE4
## IT HAS BEEN MOVED TO ADESPATIAL
## THIS IS MERELY A TEMPORARY FIX.
## only orthobasis.listw is used.
## renamed into .orthobasis.listw to avoid having to document it.
## - Thibaut Jombart, April 2015
#######################################################
.orthobasis.listw <- function( listw) {
appel = match.call()
if(!inherits(listw,"listw")) stop ("object of class 'listw' expected")
if(listw$style!="W") stop ("object of class 'listw' with style 'W' expected")
n = length(listw$weights)
fun <- function (x) {
num = listw$neighbours[[x]]
wei = listw$weights[[x]]
res = rep(0,n)
res[num] = wei
return (res)
}
b0 <- matrix(unlist(lapply(1:n,fun)),n,n)
b0=(t(b0)+b0)/2
b0=bicenter.wt(b0)
a0 <- eigen(b0, symmetric = TRUE)
#barplot(a0$values)
a0 <- a0$vectors
a0 <- cbind(rep(1,n),a0)
a0 <- qr.Q(qr(a0))
a0 <- as.data.frame(a0[,-1])*sqrt(n)
row.names(a0) <- attr(listw,"region.id")
names(a0) <- paste("VP", 1:(n-1), sep = "")
z <- apply(a0,2,function(x) sum((t(b0*x)*x))/n)
attr(a0,"values") <- z
attr(a0,"weights") <- rep(1/n,n)
attr(a0,"call") <- appel
attr(a0,"class") <- c("orthobasis","data.frame")
return(a0)
}
## print.orthobasis <- function(x,...) {
## if (!inherits(x,"orthobasis")) stop ("for 'orthobasis' object")
## cat("Orthonormal basis: ")
## n <- nrow(x)
## p <- ncol(x)
## if (n!=(p+1)) stop ("Non convenient dimension: author's error")
## cat("data.frame with",n,"rows and",ncol(x),"columns\n")
## cat("--------------------------------------\n")
## cat("Columns are an orthonormal basis of 1n-orthogonal for\n")
## cat("the inner product defined by the weights attribute\n")
## cat("---------------------------------------\n")
## w <- attributes(x)
## if (!is.null(w$"names")) cat("names =", w$names[1],"...",w$names[p],"\n")
## if (!is.null(w$"row.names")) cat("row.names =", w$row.names[1],"...",w$row.names[n],"\n")
## if (!is.null(w$"weights")) cat("weights =", w$weights[1],"...",w$weights[n],"\n")
## if (!is.null(w$"values")) cat("values =", w$values[1],"...",w$values[p],"\n")
## if (!is.null(w$"class")) cat("class =", w$class,"\n")
## if (!is.null(w$"call")) {
## cat("call =")
## print(w$"call")
## }
## }
## orthobasis.mat <- function(mat, cnw=TRUE) {
## if (!is.matrix(mat)) stop ("matrix expected")
## if (any(mat<0)) stop ("negative value in 'mat'")
## if (nrow(mat)!=ncol(mat)) stop ("squared matrix expected")
## mat <- (mat+t(mat))/2
## nlig <- nrow(mat)
## if (is.null(dimnames(mat))) {
## w <- paste("P",1:nrow(mat),sep="")
## dimnames(mat) <- list(w,w)
## }
## labels <- dimnames(mat)[[1]]
## if (cnw) {
## margi <- apply(mat,1,sum)
## margi <- max(margi)-margi
## mat <- mat+diag(margi)
## }
## mat <- mat/sum(mat)
## wt <- rep ((1/nlig),nlig)
## # calculs extensibles a une ponderation quelconque
## wt <- wt/sum(wt)
## # si mat wt est la ponderation marginale associee a mat
## # tot = sum(mat)
## # mat = mat-matrix(wt,nlig,nlig,byrow=TRUE)*wt*tot
## # encore plus particulier mat = mat-1/nlig/nlig
## # en general les precedents sont des cas particuliers
## U <- matrix(1,nlig,nlig)
## U <- diag(1,nlig)-U*wt
## mat <- U%*%mat%*%t(U)
## wt <- sqrt(wt)
## mat <- t(t(mat)/wt)
## mat <- mat/wt
## eig <- eigen(mat,sym=TRUE)
## w0 <- abs(eig$values)/max(abs(eig$values))
## tol <- 1e-07
## w0 <- which(w0<tol)
## if (length(w0)==0) stop ("abnormal output : no null eigenvalue")
## else if (length(w0)==1) w0 <- (1:nlig)[-w0]
## else if (length(w0)>1) {
## # on ajoute le vecteur derive de 1n
## w <- cbind(wt,eig$vectors[,w0])
## # on orthonormalise l'ensemble
## w <- qr.Q(qr(w))
## # on met les valeurs propres a 0
## eig$values[w0] <- 0
## # on remplace les vecteurs du noyau par une base orthonormee contenant
## # en premiere position le parasite
## eig$vectors[,w0] <- w[,-ncol(w)]
## # on enleve la position du parasite
## w0 <- (1:nlig)[-w0[1]]
## }
## mat <- eig$vectors[,w0]/wt
## mat <- data.frame(mat)
## row.names(mat) <- labels
## names(mat) <- paste("S",1:(nlig-1),sep="")
## attr(mat,"values") <- eig$values[w0]
## attr(mat,"weights") <- rep(1/nlig,nlig)
## attr(mat,"call") <- match.call()
## attr(mat,"class") <- c("orthobasis","data.frame")
## return(mat)
## }
## "orthobasis.haar" <- function(n) {
## # on definit deux fonctions :
## appel = match.call()
## a <- log(n)/log(2)
## b <- floor(a)
## if ((a-b)^2>1e-10) stop ("Haar is not a power of 2")
## # la premiere est ecrite par Daniel et elle donne la demonstration (par analogie avec la fonction qui construit la base Bscores)
## # que la base Bscores est exactement la base de Haar quand on prend une phylogenie reguliere resolue.
## "haar.basis.1" <- function (n) {
## pari <- matrix(c(1,n),1)
## "div2" <- function (mat) {
## res <- NULL
## for (k in 1 : nrow(mat)) {
## n1 <- mat[k,1]
## n2 <- mat[k,2]
## diff <- n2-n1
## if (diff <=0) break
## n3 <- floor((n1+n2)/2)
## res <- rbind(res,c(n1,n3),c(n3+1,n2))
## }
## if (!is.null(res)) pari <<- rbind(pari,res)
## return(res)
## }
## mat <- div2(pari)
## while (!is.null(mat)) mat <- div2(mat)
## res <- NULL
## for (k in 1:nrow(pari)) {
## x<-rep(0,n)
## x[(pari[k,1]):(pari[k,2])] <- 1
## res <-c(res,x)
## }
## res = matrix(res,n)
## res <- qr.Q(qr(res))
## res <- res[, -1] * sqrt(n)
## res <- data.frame(res)
## row.names(res) <- paste("u",1:n,sep="")
## names(res) <- paste("B",1:(n-1),sep="")
## return(res)
## }
## # la seconde exploite les potentialites de la librairie waveslim, en remarquant qu'il existe un lien etroit entre la definition des filtres et la definition
## # des bases. Cette strategie permettra a l'avenir de definir les bases associees a d'autres famille de fonctions.
## "haar.basis.2" <- function (n) {
## if (!require(waveslim)) stop ("Please install waveslim")
## J <- a #nombre de niveau
## res <- matrix(0, nrow = n,ncol = n-1)
## filter.seq <- "H" #filtre correspondant au niveau 1
## h <- waveslim::wavelet.filter(wf.name = "haar", filter.seq = filter.seq) #parametre du filtre au niveau 1
## k <- 0
## for(i in 1:J){
## z <- rep(h,2**(J-i))
## x <- 1:n
## y <- rep((n-1-k):(n-2**(J-i)-k),rep(2**i,2**(J-i)))
## for(j in 1:n) res[x[j],y[j]] <- z[j]
## k <- k+2**(J-i)
## filter.seq <- paste(filter.seq, "L", sep = "")
## h <- waveslim::wavelet.filter(wf.name = "haar", filter.seq = filter.seq)
## }
## res <- res*sqrt(n)
## res <- data.frame(res)
## row.names(res) <- paste("u", 1:n, sep = "")
## names(res) <- paste("B", 1:(n-1), sep = "")
## return(res)
## }
## # suivant que n est grand (n > 257) ou non, on choisit l'une des deux strategies :
## if (n < 257)
## res <- haar.basis.1(n)
## else
## res <- haar.basis.2(n)
## attr(res,"values") <- NULL
## attr(res,"weights") <- rep(1/n,n)
## attr(res,"call") <- appel
## attr(res,"class") <- c("orthobasis","data.frame")
## return(res)
## }
## "orthobasis.line" <- function (n) {
## appel <- match.call()
## # solution de Cornillon p. 12
## res <- NULL
## for (k in 1:(n-1)) {
## x <- cos(k*pi*(2*(1:n)-1)/2/n)
## x <- sqrt(n)*x/sqrt(sum(x*x))
## res <-c(res,x)
## }
## res=matrix(res,n)
## res <- data.frame(res)
## row.names(res) <- paste("u",1:n,sep="")
## names(res) <- paste("B",1:(n-1),sep="")
## w <- (1:(n-1))*pi/2/n
## valpro <- 4*(sin(w)^2)/n
## poivoisi <- c(1,rep(2,n-2),1)
## poivoisi <- poivoisi/sum(poivoisi)
## norm <- unlist(apply(res, 2, function(a) sum(a*a*poivoisi)))
## y <- valpro*n*n/2/(n-1)
## val <- norm - y
## attr(res,"values") <- val
## attr(res,"weights") <- rep(1/n,n)
## attr(res,"call") <- appel
## attr(res,"class") <- c("orthobasis","data.frame")
## # verification locale. Ce paragraphe verifie que les vecteurs et les valeurs
## # proposee par Cornillon p. 12 sont bien les vecteurs propres de l'operateur de voisinage
## # rangee dans la solution analytique par variance locale croissante
## # l'article de Meot est errone et a donne le graphe circulaire pour le graphe lineaire
## # d0=neig2mat(neig(n.lin=n))
## # d0 = d0/n
## # d1=apply(d0,1,sum)
## # d0=diag(d1)-d0
## # fun2 <- function(x) {
## # z <- sum(t(d0*x)*x)/n
## # z <- z/sum(x*x)
## # return(z)
## # }
## # lambda <- unlist(apply(res,2,fun2))
## # print(lambda)
## # print(attr(res,"values"))
## # plot(lambda,attr(res,"values"))
## # abline(lm(attr(res,"values")~lambda))
## # print(coefficients(lm(attr(res,"values")~lambda)))
## # verification que les valeurs derivees des valeurs propres sont exactement des indices de Moran
## # d = neig2mat(neig(n.lin=n))
## # d = d/sum(d) # Moran type W
## # moran <- unlist(lapply(res,function(x) sum(t(d*x)*x)))
## # print(moran)
## # plot(moran,attr(res,"values"))
## # abline(lm(attr(res,"values")~moran))
## # print(summary(lm(attr(res,"values")~moran)))
## return(res)
## }
## "orthobasis.circ" <- function (n) {
## appel = match.call()
## if (n<3) stop ("'n' too small")
## "vecprosin" <- function(k) {
## x <- sin(2*k*pi*(1:n)/n)
## x <- x/sqrt(sum(x*x))
## }
## "vecprocos" <- function(k) {
## x <- cos(2*k*pi*(1:n)/n)
## x <- x/sqrt(sum(x*x))
## }
## "valpro" <- function(k,bis=TRUE) {
## x <- (4/n)*((sin(k*pi/n))^2)
## if (bis) x <- c(x,x)
## return(x)
## }
## k <- floor(n/2)
## if (k==n/2) {
## #n est pair
## w1 <- matrix(unlist(lapply(1:k,vecprocos)),n,k)
## w2 <- matrix(unlist(lapply(1:(k-1),vecprosin)),n,k-1)
## res <- cbind(w1,w2)
## res[,seq(1,2*k-1,by=2)]<-w1
## res[,seq(2,2*k-2,by=2)]<-w2
## vp <- unlist(lapply(1:(k-1),valpro))
## vp <- c(vp, valpro(k,FALSE))
## } else {
## # n est impair
## w1 <- matrix(unlist(lapply(1:k,vecprocos)),n,k)
## w2 <- matrix(unlist(lapply(1:k,vecprosin)),n,k)
## res <- cbind(w1,w2)
## res[,seq(1,2*k-1,by=2)]<-w1
## res[,seq(2,2*k,by=2)]<-w2
## vp <- unlist(lapply(1:k,valpro))
## }
## res=sqrt(n)*res
## res <- as.data.frame(res)
## row.names(res) <- paste("u",1:n,sep="")
## names(res) <- paste("B",1:(n-1),sep="")
## attr(res,"values") <- 1 - n*vp/2
## attr(res,"weights") <- rep(1/n,n)
## attr(res,"call") <- appel
## attr(res,"class") <- c("orthobasis","data.frame")
## # verification qu'on a exactement des indices de Moran a partie des valeurs propres
## # d = neig2mat(neig(n.cir=n))
## # d = d/sum(d) # Moran type W
## # moran <- unlist(lapply(res,function(x) sum(t(d*x)*x)))
## # print(moran)
## # plot(moran,attr(res,"values"))
## # abline(lm(attr(res,"values")~moran))
## # print(summary(lm(attr(res,"values")~moran)))
## return(res)
## }
## "orthobasis.neig" <- function( neig) {
## appel = match.call()
## if(!inherits(neig,"neig")) stop ("object of class 'neig' expected")
## n <- length(attr(neig,"degree"))
## m <- sum(attr(neig,"degree"))
## poivoisi <- attr(neig,"degree")/m
## if (is.null(names(poivoisi))) names(poivoisi) <- as.character(1:n)
## d0 = neig2mat(neig)
## d0 = diag(poivoisi)-d0/m
## eig <- eigen(d0, sym = TRUE)
## ########
## tol <- 1e-07
## w0 <- abs(eig$values)/max(abs(eig$values))
## w0 <- which(w0<tol)
## if (length(w0)==0) stop ("abnormal output : no null eigenvalue")
## else if (length(w0)==1) w0 <- (1:n)[-w0]
## else if (length(w0)>1) {
## # on ajoute le vecteur derive de 1n
## wt <- rep(1,n)
## w <- cbind(wt,eig$vectors[,w0])
## # on orthonormalise l'ensemble
## w <- qr.Q(qr(w))
## # on met les valeurs propres a 0
## eig$values[w0] <- 0
## # on remplace les vecteurs du noyau par une base orthonormee contenant
## # en premiere position le parasite
## eig$vectors[,w0] <- w[,-ncol(w)]
## # on enleve la position du parasite
## w0 <- (1:n)[-w0[1]]
## }
## w0 <- rev(w0)
## valpro <- eig$values[w0]
## eig <- eig$vectors[,w0]
## eig <- as.data.frame(eig)*sqrt(n)
## z <- apply(eig,2,function(x) sum(x*x*poivoisi))
## z <- z - valpro*n
## w <- rev(order(z))
## z <- z[w]
## eig <- eig[,w]
## row.names(eig) <- names(poivoisi)
## names(eig) <- paste("VP", 1:(n-1), sep = "")
## attr(eig,"values") <- z
## attr(eig,"weights") <- rep(1/n,n)
## attr(eig,"call") <- appel
## attr(eig,"class") <- c("orthobasis","data.frame")
## return(eig)
## }
Try the adegenet package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.