Revision 4c9961af6c3e72b611474c0722c9262fb871ad46 authored by Matthias Templ on 05 July 2010, 18:25:41 UTC, committed by cran-robot on 05 July 2010, 18:25:41 UTC
1 parent b66ed0e
alrEM.R
alrEM <- function(x, pos=ncol(x), dl=rep(0.05, ncol(x)-1), eps=0.0001, maxit=50) {
stopifnot(all(x[,pos] != 0 & length(which(is.na(x[,pos]))) == 0))
m=matrix(rep(dl,each=nrow(x)),ncol=length(dl))
phi <- log(m/x[,pos])
xOrig <- x
if(length(colnames) < 1) names <- letters[1:ncol(x)] else names <- colnames(x)
nam <- colnames(x)[pos]
x <- cbind(x[, -pos], x[, pos])
w <- which(colnames(x) == "x[, pos]")
colnames(x)[w] <- nam
## close the X rows to 1:
xc <- constSum(x)
## convert zeros into missing values:
x[x==0] <- NA
## transform x into y=alr(x,pos):
# x is already in the correct order
y <- alr(x)$x.alr
y <- data.frame(y)
w <- is.na(y)
wr <- apply(y, 1, function(x) any(is.na(x)) )
wrr <- which(wr)
it <- 0
d <- 99999999
y <- as.matrix(y)
y[w] <- 0
it <- 0
amountMiss <- length(which(w))
while( d > eps & it <= maxit ){
it <- it + 1
yold <- y
for(i in 1:ncol(y)){
lm1 <- lm( y[,i,drop=FALSE] ~ y[,-i,drop=FALSE] )
yhat <- predict(lm1, new.data=y[,-i])
s <- sd(y[,i], na.rm=TRUE)
ex <- (phi[,i] - yhat)/s
yhat2 <- yhat - s*dnorm(ex)/pnorm(ex)
y[w[,i],i] <- ifelse(yhat2[w[,i]] >= phi[w[,i],i], phi[w[,i],i], yhat2[w[,i]])
}
d <- sum(abs(y)- sum(abs(yold)))/amountMiss
}
ximp <- suppressWarnings(invalr(y))
w <- which(colnames(ximp) == "rat")
colnames(ximp)[w] <- nam
colnames(ximp)[which(colnames(ximp) =="")] <- names[pos]
ximp <- ximp[, names]
res <- list(xOrig=xOrig, xImp=ximp, wind=NULL, iter=it, eps=eps)
invisible(res)
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...