[R] Problem to transfer Splus functions
Prof Brian Ripley
ripley at stats.ox.ac.uk
Mon Nov 5 08:26:28 CET 2001
On Mon, 5 Nov 2001, Michel ARNAUD wrote:
> Hello
> I would like to transfer some Splus functions in R.
> But I have a problem first about this assignation in Splus :
> xnom <- deparse(substitute(x))
You need to move such lines in R to the top of the function, before x
(etc) is altered. It's in the R FAQ.
> I am a bad programmer : I don't understand the R help
> How to modify these functions ?
> Thank you very much for your help
>
> Here are the four functions and a data test
> ####################################################################################
>
> acp <- function(x, wt = rep(1/nrow(x), nrow(x)), d = rep(1, ncol(x)),
> ctr = T,
> reduc = T, contav = F, method = "acp", tol = 10^(-8))
> {
> #============================================================
> # Methode d'analyse en composantes principales
> # avec des metriques diagonales:
> # wt vecteur des poids
> # d vecteurs des elements diagonaux de la metrique
> # si ctr=T l'acp est centree
> # si reduc=T l'acp est reduite
> # Method = "acp" or "afc" or "afcm"
> # method that produced the call to the function acp
> # tol values of the eigenvalues that are set to zero
> # to avoid zero "negative" values
> #============================================================
> x <- as.matrix(x)
> wt <- wt/sum(wt) #Normalisation des poids
> x <- crenom(x)
> nam <- dimnames(x)
> y <- centre(x, wt)
> moy <- (x - y)[1, ]
> names(moy) <- nam[[2]]
> if(ctr)
> x <- y
> y <- reduct(x, wt, ctr = ctr)
> sigm <- attr(y, "std")
> names(sigm) <- nam[[2]] # ..si reduc=T..alors
> if(reduc) x <- y
> #------------------------------controle prealable: impression des
> contributions (si contav=T)
> if(contav) {
> conta(x, d, wt)
> cat("\n")
> tex <- c("On continue ?", "Arret")
> ski <- menu(tex)
> switch(ski,
> ,
> stop())
> }
> y <- wt * x
> v <- t(x) %*% y # Calcul de la matrice a diagonaliser
> d12 <- sqrt(d)
> v <- d12 * v
> v <- t(d12 * t(v))
> res <- eigen(v, symmetric = T)
> # Recherche des elements propres de v
> #------------------------------------------------------------------------
>
> # Factors associated with very small eigenvalues removed
> #------------------------------------------------------------------------
>
> nf <- length(res$values)
> cond <- res$values < tol * sum(res$values)
> nfmax <- if(sum(cond) > 0) min((1:length(res$values))[cond]) -
> 1 else nf
> res$values <- res$values[1:nfmax]
> res$vectors <- res$vectors[, 1:nfmax]
> vecp <- res$vectors
> # Reorientation des deux premiers vecteurs propres pour avoir un max
> # de saturations positives
> v12 <- rep(1, ncol(x)) %*% (vecp[, 1:2]/abs(vecp[, 1:2]))
> v12 <- 2 * (as.numeric(v12 >= 0) - 0.5)
> vecp[, 1:2] <- t(v12 * t(vecp[, 1:2]))
> # calcul des vecteurs propres M normes
> res$vectors <- (1/sqrt(d)) * vecp
> # Generation des noms de facteurs
> nomf <- as.vector(outer("f", 1:nfmax, paste, sep = ""))
> dimnames(res$vectors) <- list(nam[[2]], nomf)
> #Affectation des noms
> # des vecteurs principaux
> res$cmpr <- x %*% (d12 * vecp)
> #Calcul des composantes principales
>
> xnom <- deparse(substitute(x))
> wtnom <- deparse(substitute(wt))
> wtequal <- F
> if(wtnom == paste("rep(1/nrow(", xnom, "), nrow(", xnom, "))",
> sep = ""))
> wtequal <- T
> dnom <- deparse(substitute(d))
> dusual <- F
> if(dnom == paste("rep(1, ncol(", xnom, "))", sep = ""))
> dusual <- T
> dimnames(res$cmpr)[[2]] <- nomf
> names(res$values) <- nomf
> res$d <- d
> res$pi <- wt
> res$reduc <- reduc
> res$ctr <- ctr
> res$moy <- moy
> res$sigma <- sigm
> attr(res, "class") <- "acp"
> res$xnom <- xnom
> res$wtequal <- wtequal
> if(!wtequal)
> res$wtnom <- wtnom
> res$dusual <- dusual
> if(!dusual)
> res$dnom <- dnom
> if(method == "acp")
> print(res)
> res
> }
>
> ####################################################################################
>
> crenom <- function(x, nl = "i", nv = "v")
> {
> #--------------------------------------------------------------------------
>
> # Cree des noms pour le tableau x si ces noms n'existent pas
> # Par defaut, les noms de ligne commencent par "i", ceux de colonnes
> # par "v"
> #--------------------------------------------------------------------------
>
> nomi <- if((!is.null(dimnames(x)) && !is.null(dimnames(x)[[1]])
> ) && length(dimnames(x)[[1]]) != 0) dimnames(x)[[1]]
> else (paste(nl, 1:nrow(x), sep = ""))
> nomv <- if(!is.null(dimnames(x)) && !is.null(dimnames(x)[[2]]) &&
> length(dimnames(x)[[2]]) != 0) dimnames(x)[[2]] else (
> paste(nv, 1:ncol(x), sep = ""))
> dimnames(x) <- list(nomi, nomv)
> x
> }
> ####################################################################################
>
> centre <- function(x, wt = rep(1, nrow(as.matrix(x))))
> {
> # centre le nuage des lignes de x ponderees par wt
> x <- as.matrix(x)
> g <- as.vector((wt %*% x)/sum(wt))
> t(t(x) - g)
> }
> ####################################################################################
>
> reduct <- function(x, wt = rep(1/nrow(x), nrow(x)), ctr = T)
> {
> # Calcul d'une matrice reduite (avec poids)
> #
> # wt vecteur des poids
> #============================================================
> nc <- ncol(x)
> nr <- nrow(x)
> nam <- dimnames(x)
> wt <- wt/sum(wt)
> mx <- wt %*% x
> xx <- t(t(x) - as.vector(mx))
> rwt <- sqrt(wt)
> y <- rwt * xx #.. integration des poids..
> stdx <- sqrt(diag(t(y) %*% y))
> if(ctr)
> x <- xx
> x <- t(t(x)/stdx)
> dimnames(x) <- nam
> attr(x, "std") <- stdx
> x
> }
> ca mg cl so4 hco3 co2 c te deb
> 62.8 3.5 1.75 27.5 173 2.93 300 10.40 0.092
> 62.8 4.0 1.65 27.0 177 2.92 300 10.40 0.490
> 64.5 3.6 1.65 28.5 179 3.10 303 10.40 0.560
> 64.0 5.0 1.70 26.0 184 3.33 308 10.40 0.605
> 61.6 5.0 1.90 28.0 178 3.08 303 10.40 0.690
> 63.5 4.4 1.75 27.5 179 2.76 302 10.50 0.740
> 66.8 4.4 1.65 31.5 183 2.82 311 10.50 0.760
> 69.3 4.4 1.85 37.0 188 3.02 330 10.45 0.810
> 66.3 5.5 1.80 35.5 186 2.92 327 10.45 0.830
> 66.5 4.7 1.75 32.5 182 2.86 317 10.45 0.860
> 66.1 3.3 1.80 28.5 181 2.78 307 10.45 0.890
> 62.8 3.2 1.70 23.0 177 2.79 293 10.40 0.950
> 62.3 3.8 1.60 23.0 177 1.92 294 10.40 1.050
> 69.0 3.9 1.85 35.5 185 3.44 324 10.35 1.190
> 71.5 5.0 1.95 38.5 195 3.94 344 10.35 1.220
> 73.2 5.5 2.25 40.5 208 4.38 352 10.35 1.220
> 71.7 5.7 2.20 35.0 204 4.94 347 10.35 1.240
> 73.0 4.9 2.30 36.5 200 4.64 348 10.35 1.260
> 69.3 4.5 2.25 29.5 193 4.38 333 10.35 1.350
> 64.5 3.3 2.00 18.5 187 4.25 306 10.35 1.630
> 60.8 3.9 2.00 14.0 183 4.08 292 10.30 1.830
> 62.5 3.8 2.45 12.5 190 4.81 296 10.10 1.580
> 61.7 3.4 2.25 11.0 187 4.77 288 10.05 2.210
> 62.9 3.0 2.20 10.5 189 4.61 291 10.05 1.930
> 66.0 2.7 2.30 9.5 198 3.49 301 10.10 1.200
> 66.5 2.3 2.25 9.5 198 3.65 303 10.05 1.050
> 65.5 3.0 2.40 11.0 198 3.82 304 10.05 0.905
> 65.5 2.5 2.30 11.0 195 3.61 298 10.05 0.490
> 64.7 3.0 2.25 9.5 195 3.69 297 10.05 0.495
> 66.5 3.4 2.25 13.0 198 3.91 308 10.10 0.830
> 66.9 3.5 2.50 14.0 199 4.21 309 10.05 1.010
> 64.5 3.4 2.20 11.5 192 3.79 298 10.00 0.920
> 63.7 3.3 2.20 11.5 189 3.92 297 10.00 0.650
> 63.7 3.4 2.30 13.0 190 3.95 297 10.10 0.560
>
>
> --
> Michel ARNAUD
> CIRAD
> TA60/15
> 73, av. Jean François Breton
> 34938 MONTPELLIER CEDEX 5
> tel 04 67 59 38 34 - Fax 04 67 59 38 27
>
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272860 (secr)
Oxford OX1 3TG, UK Fax: +44 1865 272595
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list