lmc <- function(CV,fun,a,tol=0.001,mode='aut')
{
#Valores iniciales
n <- length(CV$data[,1]) #Número de elementos medidos
k <- length(CV[[1]][,1]) #Número de lags
m <- length(fun)+1 #Número de funciones + Nugget
p <- CV$variables #Número de variables medidas
pin <- (p*(1+p)/2)
#Paquetes necesarios
require(matrixcalc)
#Ajustar b y nug de manera natural
#RESULT <- GmodCross2(CV,fun,a,mode)
RESULT <- GmodCross3(CV,fun,a)
#Matrices A y Dist
Dist <- as.matrix(dist(CV$data[,1:2]))
A <- Amatrix(CV$data[,1:2],n=length(CV[[1]][,1]))
#Calcular TRAZA
#Matriz kxk trazas inicial (depende de fun, a y A(h))
##Combinación funciones y h
CombFUN <- t(cbind(rbind(1:m,1:m),combn(1:m,2))) #Combinaciones entre funciones
CombH <- t(cbind(rbind(1:k,1:k),combn(1:k,2))) #Combinaciones entre lags
CombFUNH <- expand.grid(1:m,1:k) #Mezcla
names(CombFUNH) <- c('FUN','H')
##Matrices nxn de las funciones "tipo"
###OJO! Si hay más de una función repetida hay que tener en cuenta el rango de la misma
#Cambios en las funciones el 20/05/2016 En el github sigue la versión antigua.
CorrSph <- function(D,as){Fsph <- (ifelse(D>as,Fsph <- 0,Fsph <- (1-1.5*(D/as)+0.5*(D/as)^3)))}#MIRAR SI ES 1-1.5...
CorrExp <- function(D,as){Fexp <- (exp(-D/as))}
CorrGau <- function(D,as){Fgau <- (exp(-(D/as)^2))}
CorrLin <- function(D,as){Flin <- (ifelse(D>as,Flin <- 0,Flin <- (1-(D/as))))} #lineal
CorrPow <- function(D,as){Fpow <- (ifelse(as<=2 & as>0,Fsph <- (1-D^as),stop('el rango no se encuentra entre 0 y 2')))}
##FUN
FUN <- list(diag(1,n,n))
for (i in 1:length(fun))
{if(fun[i]=='Exp'){FUN[[i+1]]<-CorrExp(Dist,a[i])}
else{if(fun[i]=='Sph'){FUN[[i+1]]<-CorrSph(Dist,a[i])}else{FUN[[i+1]]<-CorrGau(Dist,a[i])}}}
names(FUN) <- c('Nug',fun)
prevtraza <- function(v1,v2)
{list(FUN[[v1]]%*%A[[v2]])}
FUNH <- mapply(prevtraza,CombFUNH[,1],CombFUNH[,2])
traza <- function(w1)
{
ind1 <- rep(which(CombFUNH$H==CombH[w1,1]), each=length(fun))
ind2 <- rep(which(CombFUNH$H==CombH[w1,2]),each=length(fun))
r <- matrix.trace(FUNH[[ind1[1]]]%*%FUNH[[ind2[1]]])
for (i in 2:length(ind1))
{r <- r+matrix.trace(FUNH[[ind1[i]]]%*%FUNH[[ind2[i]]])}
t <- r
}
CombH <- cbind(CombH,unlist(lapply(1:length(CombH[,1]),traza)))
TRAZAS <- matrix(0,k,k)
TRAZAS[lower.tri(TRAZAS,diag=TRUE)] <- CombH[,3]
TRAZAS[upper.tri(TRAZAS)] = t(TRAZAS)[upper.tri(TRAZAS)] #Matriz kxk TRAZAS
#Calcular la primeras matrices de B0 a partir de RESULT
comb <- t(cbind(rbind(1:p,1:p),combn(1:p, 2)))
comb <- cbind(comb,rep(0,(p*(1+p)/2)))
Bnew <- RESULT$B
Bnewlast <- Bnew
bcalc <- function(B)
{sum((Reduce('+',B))^2)} #Cálculo de la constante b que multiplica a la matriz de trazas
b <- bcalc(Bnew)
WSS <- function(index)
{bg <- Bnew[index,1]*RESULT$FUN[,1]
for (i in 2:m)
{bg <- bg+Bnew[index,i]*RESULT$FUN[,i]}
wss <- t(CV[[index]][,2]-bg)%*%solve(b*TRAZAS)%*%(CV[[index]][,2]-bg)#ANTES TENIA bg[[index]]
}
SUMWSS <- function(comb,wss)
{
s <- matrix(0,p,p)
comb[,3] <- wss
s[comb[,1:2]] <- comb[,3]
s[lower.tri(s)] = t(s)[lower.tri(s)]
sum(s)
}
wss0 <- unlist(lapply(1:pin,WSS))
WSS1 <- SUMWSS(comb,wss0)
WSSInic <- WSS1
WSS0 <- 0
getbg <- function(index)
{
bg <- Bnewlast[index,1]*RESULT$FUN[,1]
for (i in 2:m)
{bg <- rbind(bg,Bnewlast[index,i]*RESULT$FUN[,i])}
rownames(bg) <- c('Nug',fun)
t(bg)
}
getgammas0 <- function(index, funi)
{bg <- lapply(1:pin,getbg)
bgd <- CV[[index]][,2]-(rowSums(bg[[index]])-bg[[index]][,funi])
bgd}
fitmodel <- function(index,funi)
{
gammas0 <- mapply(getgammas0,1:pin,funi)
bnew <- solve(t(RESULT$FUN[,funi])%*%solve(b*TRAZAS)%*%RESULT$FUN[,funi])%*%t(RESULT$FUN[,funi])%*%solve(b*TRAZAS)%*%gammas0[,index]
}
Bt <- function(comb,bfit)
{
comb[,3] <- bfit
s <- matrix(0,p,p)
s[comb[,1:2]] <- comb[,3]
s[lower.tri(s)] = t(s)[lower.tri(s)]
s
}
DVS <- function(BN)
{PC <- eigen(BN)
PC$values[PC$values<0] <- 0
ord <- order(PC$values,decreasing = T)
PC$values <- PC$values[ord]
PC$vectors <- PC$vectors[,ord]
dG <- PC$vectors%*%diag(PC$values)%*%t(PC$vectors)}
#Inversa generalizada
Bvector <- function(B)
{bp <- c(diag(B),B[lower.tri(B)])
bp}
it <- 0
while (abs(WSS1-WSS0)>tol)
{
s <- 1
while (s<=m)
{
bfit <- mapply(fitmodel,1:pin,s)
Bfit <- Bt(comb,bfit)
Bdvs <- DVS(Bfit)
Bnew[,s] <- Bvector(Bdvs)
b <- bcalc(Bnew)
s <- s+1
}
wss1 <- unlist(lapply(1:(p*(1+p)/2),WSS))
WSS2 <- SUMWSS(comb,wss1)
if (abs(WSS1-WSS0) < abs(WSS2-WSS1)) break
if (abs(WSS2)>abs(WSS1)) break
WSS0 <- WSS1 #EL FALLO ESTÁ EN EL CONTROL DE WSS OJO!
WSS1 <- WSS2
Bnewlast <- Bnew
it <- it+1
#El problema es que WSS1-WSS0 puede ser mayor que el anterior, entonces no para y se aleja
}
orderBfinal <- function(comb,BT)
{
B <- list(0)
s <- matrix(0,p,p)
for (i in 1:m)
{
comb[,3] <- BT[,i]
s[comb[,1:2]] <- comb[,3]
s[lower.tri(s)] = t(s)[lower.tri(s)]
B[[i]] <- s
}
names(B) <- c('Nug',fun)
B
}
Bfinal <- orderBfinal(comb,Bnewlast)
if(it==0)
{
for (i in 1:m)
{
Bfinal[[i]] <- DVS(Bfinal[[i]])
}
for (i in 1:m)
{
Bnewlast[,i] <- c(diag(Bfinal[[i]]),Bfinal[[i]][lower.tri(Bfinal[[i]])])
}
}
NEWRESULT <- list(var=CV$variables,rango=a,B=Bfinal,Bvector=Bnewlast,FUN=RESULT$FUN,WSSini=WSSInic,WSSfinal=WSS1,iter=it,data=CV$data)
#if( NEWRESULT$WSSini<NEWRESULT$WSSfinal ) stop('Valores iniciales mal elegidos. Prueba con nuevas combinaciones')
NEWRESULT
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.