nrbf.R

nrbf <- function(X, Y, X_test = X, k = ncol(X), gamma = 1.0, seed = 123, plot = TRUE){
  library(corpcor)
  library(neuralnet)
  set.seed(seed)
  #### Definição dos argumentos:
  # X: Matriz de input 
  # Y: Matriz de output
  # k: número de centros (polos)
  # gamma: parametro de aprendizado
  
  N <- nrow(X) # numero de observacoes
  # criar funcao de ativacao gaussiana
  ativ_gaussiana <- function(gamma, x, y){
    # x e y sao vetores numericos
    
    v <- as.matrix(x - y)
    modulo <- norm(v, "F")
    exp(-gamma * modulo^2)
  }
  
    
  # criar clusteres em matriz de input normalizada
  cl <- kmeans(scale(X), centers = k)
  cl_centros <- cl$centers # (Matriz de dimensões k x ncol(X))
  # criar matriz vazia Phi para preencher  posteriormente
  matriz_phi <- matrix(NA_real_, nrow = N, ncol = k + 1)
  
  #browser()
  # iterar em cada linha
  for (l in 1:N){
    # Preencher a primeira coluna com 1 (coluna de bias)
    matriz_phi[l, 1] <- 1
    # iterar em cada coluna 
    for (c in 1:k){
      # calcular modulo do vetor das diferencas
      # preencher celula (com excecao da primeira coluna, que eh a de bias)
      matriz_phi[l, c + 1] <- ativ_gaussiana(gamma = gamma, x = X[l, ], y = cl_centros[c, ])
    }
  }
  
  w <- corpcor::pseudoinverse(t(matriz_phi) %*% matriz_phi) %*% t(matriz_phi) %*% Y 
  
  # obter previsoes para Y
  # inicializar vetor de previsoes a partir do bias do vetor de pesos
  prev <- rep(w[1], N)
  for (n in 1:N){
    for (j in 1:nrow(cl_centros)){
      prev[n] <- prev[n] + w[j + 1] * ativ_gaussiana(gamma = gamma, x = X_test[n, ], y = cl_centros[j, ])
    }
  }
  
  # plotar rede neural
  if (plot){
    # combinar input e target em um dataframe
    df_nnet <- data.frame(input, target = target)
    model_formula <- paste0("target ~ ", paste0(colnames(input), collapse = " + "))
    nn <- neuralnet(model_formula, data = df_nnet, hidden = k)
    plot(nn, information = FALSE, show.weights = FALSE)
  }
  
  list(weights=w, centers=cl_centros, gamma=gamma, prev = prev)  
}



#Funçao de base radial
funcao_base <- function(x, sigma) {exp(1) ^ (-0.5 * sigma * x)} 


funcao_radial <- function(calculo, x, sigma, medioides) {
  
  n <- nrow(calculo)
  res <- vector("numeric", length = n)
  for (i in 1:(n-1)){
    res[i] <- calculo[i,1] * funcao_base((sum((x - medioides[i,]) ^ 2) ^ 0.5) ^ 2,
                                         sigma)
  }
  res[n] <- calculo[n, 1]
  sum(res)
  
}

# carregar codigo dentro de loop
extrair_r2 <- function(seed, k){
  #k <- ncol(dados) - 1
  
  k <- k
  # set.seed(seed)
  # medioides <- kmeans(dados,k, algorithm = "MacQueen")
  # 
  # medioides <- medioides$centers
  # 
  # medioides <- as.matrix(medioides)
  # 
  # polo1 <- rep(0, ncol(dados))
  # 
  # medioides <- rbind(polo1, medioides)
  
  
  #Defini?ao de medioides atrav?s de aleat?rios
  
  
  medioides <- matrix(data=runif((ncol(dados)*k), 0, 1), nrow=k, ncol = ncol(dados)) 
  
  polo1 <- rep(0, ncol(dados))
  
  medioides <- rbind(polo1, medioides)
  
  ##############################################
  
  #Cálculo do Sigma
  
  sigma    = matrix(data=NA, nrow=nrow(medioides), ncol = 1) #criaçao da matriz
  
  for(i in 2:nrow((medioides - 1) ^ 2 + (medioides - 1) / 2))  { 
    for(j in 2:nrow((medioides - 1) ^ 2 + (medioides - 1) / 2)) {
      
      sigma[i] <- (sum ((medioides[i,] - medioides[j-1,]) ^ 2) ^ 0.5) ^ 2 
      
    }}
  
  sigma <- max(sigma, na.rm=TRUE) #escolha do máximo entre todas as distancias calculadas
  #Nota: talvez seja melhor calular um sigma para cada RBF, ao invés de usar o máximo para todas
  
  sigma_2 <- sigma ^ 2 #cálculo feito com um único sigma para todas as RBF
  ##############################################
  
  
  # Criaçao da matriz de distancias utilizando norma euclidiana às quais a funçao será aplicada
  
  base_de_dados     = matrix(data=NA, nrow=nrow(dados), ncol= nrow(medioides))
  
  for(i in 1:nrow(dados)) { 
    for(j in 1:nrow(medioides)) {
      
      base_de_dados[i,j] <- (sum ((dados[i,] - medioides[j,]) ^ 2) ^ 0.5) ^ 2 
      
    }}
  
  
  
  
  #aplicaçao da funçao
  nova_base <- apply(base_de_dados, c(1,2), funcao_base, sigma = sigma_2)
  ##############################################
  
  #adi?ao dos coefiecientes
  
  coeficiente <- rep(1, nrow(dados))
  nova_base <- as.matrix(cbind(nova_base, coeficiente), nrow = nrow(dados), ncol = ncol(medioides + 1))
  
  
  #multiplicaçao da matriz pela sua transposta, calculo da inversa e multiplicaçao pela transposta
  # (pseudo-inversa)
  transposta <- t(nova_base) %*% nova_base
  #inversa <- solve(transposta)
  #pseudo <- inversa %*% t(nova_base)
  pseudo <- pseudoinverse(transposta) %*% t(nova_base)
  
  
  #Resultado - multiplicaçao da pseudo inversa pela resultados observados
  
  
  calculo <- pseudo %*% resultado
  ##############################################
  
  #Funcao resultado 

  
  valores_estimados <- c()
  
  for(i in 1:nrow(dados)) {
    
    valores_estimados[i] <- funcao_radial(calculo, dados[i,], sigma_2, medioides) 
    
  }
  
  valores_estimados <- as.matrix(valores_estimados)
  
  media_resultados <- colMeans(resultado)
  
  residuo <- valores_estimados - media_resultados
  
  residuo <- residuo ^ 2
  
  soma_residuo <- colSums(residuo)
  
  
  #Valores observados
  
  residuo_obs <- resultado - media_resultados
  
  residuo_obs_2 <- residuo_obs ^ 2
  
  soma_residuo1 <- colSums(residuo_obs_2)
  
  r_2 <- soma_residuo / soma_residuo1
  
  # retornar r2 obtido
  #unname(r_2)
  out <- list(
    coef = calculo,
    r_2 = unname(r_2),
    sigma_2 = unname(sigma_2),
    medioides = medioides
  )
  out
}


# carregar codigo dentro de loop
extrair_r2b <- function(X, y, seed, k){
  #k <- ncol(dados) - 1
  
  k <- k
  # set.seed(seed)
  # medioides <- kmeans(dados,k, algorithm = "MacQueen")
  # 
  # medioides <- medioides$centers
  # 
  # medioides <- as.matrix(medioides)
  # 
  # polo1 <- rep(0, ncol(dados))
  # 
  # medioides <- rbind(polo1, medioides)
  
  
  #Defini?ao de medioides atrav?s de aleat?rios
  
  
  medioides <- matrix(data=runif((ncol(X)*k), 0, 1), nrow=k, ncol = ncol(X)) 
  
  polo1 <- rep(0, ncol(X))
  
  medioides <- rbind(polo1, medioides)
  
  ##############################################
  
  #Cálculo do Sigma
  
  sigma    = matrix(data=NA, nrow=nrow(medioides), ncol = 1) #criaçao da matriz
  
  for(i in 2:nrow((medioides - 1) ^ 2 + (medioides - 1) / 2))  { 
    for(j in 2:nrow((medioides - 1) ^ 2 + (medioides - 1) / 2)) {
      
      sigma[i] <- (sum ((medioides[i,] - medioides[j-1,]) ^ 2) ^ 0.5) ^ 2 
      
    }}
  
  sigma <- max(sigma, na.rm=TRUE) #escolha do máximo entre todas as distancias calculadas
  #Nota: talvez seja melhor calular um sigma para cada RBF, ao invés de usar o máximo para todas
  
  sigma_2 <- sigma ^ 2 #cálculo feito com um único sigma para todas as RBF
  ##############################################
  
  
  # Criaçao da matriz de distancias utilizando norma euclidiana às quais a funçao será aplicada
  
  base_de_dados     = matrix(data=NA, nrow=nrow(X), ncol= nrow(medioides))
  
  for(i in 1:nrow(X)) { 
    for(j in 1:nrow(medioides)) {
      
      base_de_dados[i,j] <- (sum ((X[i,] - medioides[j,]) ^ 2) ^ 0.5) ^ 2 
      
    }}
  
  
  
  
  #aplicaçao da funçao
  nova_base <- apply(base_de_dados, c(1,2), funcao_base, sigma = sigma_2)
  ##############################################
  
  #adi?ao dos coefiecientes
  
  coeficiente <- rep(1, nrow(X))
  nova_base <- as.matrix(cbind(nova_base, coeficiente), nrow = nrow(X), ncol = ncol(medioides + 1))
  
  
  #multiplicaçao da matriz pela sua transposta, calculo da inversa e multiplicaçao pela transposta
  # (pseudo-inversa)
  transposta <- t(nova_base) %*% nova_base
  #inversa <- solve(transposta)
  #pseudo <- inversa %*% t(nova_base)
  pseudo <- corpcor::pseudoinverse(transposta) %*% t(nova_base)
  
  
  #Resultado - multiplicaçao da pseudo inversa pela resultados observados
  
  # browser()
  calculo <- pseudo %*% y
  ##############################################
  
  #Funcao resultado 
  
  
  valores_estimados <- c()
  
  for(i in 1:nrow(X)) {
    
    valores_estimados[i] <- funcao_radial(calculo, X[i,], sigma_2, medioides) 
    
  }
  
  valores_estimados <- as.matrix(valores_estimados)
  
  media_resultados <- colMeans(y)
  
  residuo <- valores_estimados - media_resultados
  
  residuo <- residuo ^ 2
  
  soma_residuo <- colSums(residuo)
  
  
  #Valores observados
  
  residuo_obs <- y - media_resultados
  
  residuo_obs_2 <- residuo_obs ^ 2
  
  soma_residuo1 <- colSums(residuo_obs_2)
  
  r_2 <- soma_residuo / soma_residuo1
  
  # retornar r2 obtido
  #unname(r_2)
  out <- list(
    coef = calculo,
    r_2 = unname(r_2),
    sigma_2 = unname(sigma_2),
    medioides = medioides
  )
  out
}