This is my current code: library(readr) library(ggplot2) library(dplyr) library(GeoModels) library(gridExtra) library(geoR) library(geosphere) library(viridis) library(fields) setwd("/data/fcp-X79-03/Proyecto_Helena/") #CARGAR DATOS. datos_st28 <- read.csv("dataStadium28.csv") head(datos_st28) #VARIABLES. datos_base1 <- subset(datos_st28, BS == 1) datos_base2 <- subset(datos_st28, BS == 2) dist3B1 <- datos_base1$dEuc dist3B2 <- datos_base2$dEuc dist2B1 <- sqrt((datos_base1$x)^2 + (datos_base1$y)^2) dist2B2 <- sqrt((datos_base2$x)^2 + (datos_base2$y)^2) datos_base1_coords <- data.frame(x = datos_base1$x, y = datos_base1$y) datos_base2_coords <- data.frame(x = datos_base2$x, y = datos_base2$y) pdf("plx.pdf") plot(datos_base2$PLOmni,datos_base2$x, xlab = "x", ylab = "Path loss") dev.off() pdf("ply.pdf") plot(datos_base2$PLOmni,datos_base2$y, xlab = "y", ylab = "Path loss") dev.off() media_muestral = mean(datos_base2$PLOmni) #REGRESIÓN lm.fit <- lm(datos_base2$PLOmni ~ dist2B2) residuos <- lm.fit$residuals pdf("hisres.pdf") hist(residuos, main = "", xlab = "Residuos centrados", ylab = "Frecuencia") dev.off() pdf("qqres.pdf") qqnorm(residuos, xlab = "Cuantiles teóricos", ylab = "Cuantiles muestrales", main = "") qqline(residuos, col = "red", lwd = 2) dev.off() #VARIOGRAMA EMPÍRICO variograma_empirico <- variog(coords = datos_base2_coords, data = residuos, max.dist = 75) pdf("varioresiduos.pdf") plot(variograma_empirico, xlab = "h", ylab = expression(hat(gamma)(h))) dev.off() geo_residuos <- as.geodata(cbind(datos_base2_coords, residuos), coords.col = c("x", "y"), data.col = 3) vg_direccional_res <- variog4(geo_residuos, max.dist = 75) plot(vg_direccional_res, xlab = "h", ylab = expression(hat(gamma)(h))) #ANISOTROPÍA MUESTRA geo_data2 <- as.geodata(datos_base2, coords.col = c("x", "y"), data.col = "PLOmni") vg_direccional_B2 <- variog4(geo_data2, max.dist = 75) # Reducir el máximo pdf("anisgeomuestra.pdf") plot(vg_direccional_B2, xlab = "h", ylab = expression(hat(gamma)(h))) dev.off() #ANISOTROPÍA RESIDUOS datos_base2$resid <- residuos geo_data2 <- as.geodata(datos_base2, coords.col = c("x", "y"), data.col = "resid") vg_direccional_B2_res <- variog4(geo_data2, max.dist = 75) # Reducir el máximo pdf("anisgeo.pdf") plot(vg_direccional_B2_res, xlab = "h", ylab = expression(hat(gamma)(h))) dev.off() #MATRIZ DE DISEÑO DE LA REGRESIÓN matriz_diseño <- model.matrix(lm.fit) beta_hat <- coef(lm.fit) param_media <- matriz_diseño %*% beta_hat mean=param_media variograma_matern <- function(h, nugget, sigma2, phi, nu) { #Función de la covarianza de Matern. matern_corr <- function(h, phi, nu) { h[h == 0] <- 1e-10 # Caso h = 0. p1 <- (2^(1 - nu)) / gamma(nu) p2 <- (h / phi)^nu p3 <- besselK(h / phi, nu) return(p1 * p2 * p3) } # Calculamos el variograma. rho_h <- matern_corr(h, phi, nu) gamma_h <- nugget + (sigma2-nugget) * (1 - rho_h) return(gamma_h) } #PRIMERA ESTIMACIÓN: USAMOS \beta_0 y \beta_1 para estimar sill, nugget y escala. sill =1; nugget = 0.6 corrmodel = "Matern" scale = 45; smooth =0.5 param = list( nugget = nugget, scale =scale , sill = sill , smooth = smooth, mean=mean ) optimizer ="nlminb" I= Inf fixed = list(smooth = smooth, mean = 0) start = list(scale = scale, sill = var(residuos), nugget=nugget) lower = list(scale =0 , sill =0 , nugget=0) upper = list(scale =I, sill =I, nugget = I) fitML1 <- GeoFit(data = residuos , coordx =datos_base2_coords , corrmodel = corrmodel , model ="Gaussian", optimizer = optimizer , lower =lower , upper =upper , likelihood ="Full", type ="Standard", start =start , fixed = fixed) fitted1 <- matriz_diseño%*%beta_hat plot(matriz_diseño[,2], datos_base2$PLOmni, asp = 1) points(matriz_diseño[,2], fitted1,col = "red", pch = 16) my.res1 <- datos_base2$PLOmni - fitted1 res1 = GeoResiduals(fitML1) GeoQQ(res1) vario1 = GeoVariogram(data = my.res1 , coordx = datos_base2_coords , maxdist =75) plot(vario1$centers, vario1$variograms, ylim = c(0, 3), type = "b", xlab = "h", ylab = expression(hat(gamma)(h))) h0 <- seq(0,160, l = 10000) toplot1 <- GeoCorrFct(h0, corrmodel = "Matern", model = "Gaussian", param = c(fitML1$param, fixed), variogram = T) lines(h0, toplot1$corr, col = "red", type = "l") gamma_1 <- variograma_matern(h0, fitML1$param$nugget,fitML1$param$sill,fitML1$param$scale,smooth) lines(h, gamma_1, lwd = 2, col = "purple") #ESTIMACIÓN MEDIA NO CONSTANTE: Estimamos todos. sill =1; nugget =0.6 corrmodel = "Matern" mean = 87.83642 mean1 = 0.1007216 scale =45; smooth =0.5 param = list( nugget = nugget, scale =scale , sill = sill , smooth = smooth , mean= mean, mean1=mean1) optimizer ="nlminb" I= Inf fixed = list(smooth = smooth) start = list(mean = mean, mean1 = mean1, scale = scale, sill = sill, nugget = nugget) lower = list(mean = 0, mean1 = 0, scale =0 , sill =0 , nugget=0) upper = list(mean = I, mean1 = I,scale =I, sill =I, nugget = I) fitML3 <- GeoFit(data = datos_base2$PLOmni, coordx =datos_base2_coords , corrmodel = corrmodel , model ="Gaussian", optimizer = optimizer , lower =lower , upper =upper , likelihood ="Full", type ="Standard", start =start , fixed = fixed, X=matriz_diseño) fitted3 <- matriz_diseño%*%c(fitML3$param$mean, fitML3$param$mean1) plot(matriz_diseño[,2], datos_base2$PLOmni, asp = 1) points(matriz_diseño[,2], fitted3,col = "red", pch = 16) my.res3 <- datos_base2$PLOmni - fitted3 res3 = GeoResiduals(fitML3) GeoQQ(res3) vario3 = GeoVariogram(data = my.res3 , coordx = datos_base2_coords , maxdist =75) pdf("medianoconstante.pdf") plot(vario3$centers, vario3$variograms, ylim = c(0, 3), type = "b", xlab = "h", ylab = expression(hat(gamma)(h))) h0 <- seq(0,160, l = 10000) toplot3 <- GeoCorrFct(h0, corrmodel = "Matern", model = "Gaussian", param = c(fitML3$param, fixed), variogram = T ) #lines(h0, toplot3$corr, col = "red", type = "l") gamma_3 <- variograma_matern(h0, fitML3$param$nugget,fitML3$param$sill,fitML3$param$scale,smooth ) lines(h, gamma_3, lwd = 2, col = "red") dev.off() #ESTIMACIÓN MEDIA CONSTANTE: Estimamos todos. sill =1; nugget =0.6 corrmodel = "Matern" mean = mean(datos_base2$PLOmni) scale = 45; smooth =0.5 param = list( nugget = nugget, scale =scale , sill = sill, mean= mean,smooth=smooth) optimizer ="nlminb" I= Inf #fixed = list(smooth = smooth) start = list(mean = mean, scale = scale, sill = var(datos_base2$PLOmni), nugget = nugget,smooth=smooth) lower = list(mean = 0, scale =0 , sill =0 , nugget=0,smooth=0) upper = list(mean = I,scale =I, sill =I, nugget = I,smooth=I) fitML5 <- GeoFit(data = datos_base2$PLOmni , coordx =datos_base2_coords , corrmodel = corrmodel , model ="Gaussian", optimizer = optimizer , lower =lower , upper =upper , likelihood ="Full", type ="Standard", start =start) res5 = GeoResiduals(fitML5) GeoQQ(res5) vario5 = GeoVariogram(data = datos_base2$PLOmni , coordx = datos_base2_coords , maxdist =150, numbins = 10) pdf("mediaconstante.pdf") plot(vario5$centers, vario5$variograms, #ylim = c(0, 1.5), type = "b", xlab = "h", ylab = expression(hat(gamma)(h)), ylim = c(0, 15)) h0 <- seq(0,160, l = 10000) toplot5 <- GeoCorrFct(h0, corrmodel = "Matern", model = "Gaussian", param = c(fitML5$param), variogram = T) #lines(h0, fitML5$param$sill*toplot5$corr, col = "red", type = "l") gamma_5 <- variograma_matern(h, fitML5$param$nugget, fitML5$param$sill, fitML5$param$scale, fitML5$param$smooth) lines(h0, gamma_5 + 3, lwd = 2, col = "red") dev.off() rango_práctico= -log(0.5)*fitML3$param$scale rango_práctico #Plano a predecir. x_val <- seq(min(datos_base2_coords$x), max(datos_base2_coords$x), length.out = 50) y_val <- seq(min(datos_base2_coords$y), max(datos_base2_coords$y), length.out = 50) grid <- expand.grid(x = x_val, y = y_val) #Datos geo_data2 <- as.geodata(datos_base2, coords.col = c("x", "y"), data.col = "PLOmni") #KRIGGING MEDIA CONSTANTE_5. # Extraemos parámetros desde fitML5. sill_5 <- fitML5$param$sill range_5 <- fitML5$param$scale nugget_5 <- fitML5$param$nugget kappa_5 <- fitML5$param$smooth modelo_5 <- list(cov.pars = c(sill_5, range_5), nugget = nugget_5, kappa = kappa_5, cov.model = "matern") krigord_5 <- krige.conv(geodata = geo_data2, locations = grid, krige = krige.control(type.krige = "ok", cov.model = modelo_5$cov.model, cov.pars = modelo_5$cov.pars, nugget = modelo_5$nugget, kappa = modelo_5$kappa)) #Graficamos predicciones del krigging. pred_5 <- cbind(grid, pred = krigord_5$predict) x_pred5 <- sort(unique(pred_5$x)) y_pred5 <- sort(unique(pred_5$y)) z_matrixpred5 <- matrix(pred_5$pred, nrow = length(x_pred5), ncol = length(y_pred5), byrow = FALSE) pdf("krigging ordinario5.pdf") par(mar = c(5, 5, 4, 7.2)) image(x_pred5, y_pred5, z_matrixpred5, col = viridis::viridis(100), xlab = "x", ylab = "y") image.plot(z = z_matrixpred5, zlim = c(min(z_matrixpred5), max(z_matrixpred5)), col = viridis::viridis(100), legend.lab = "Path loss", horizontal = FALSE, legend.only = TRUE, legend.line = 2) dev.off() #Graficamos varianza del krigging. var_5 <- cbind(grid, var = krigord_5$krige.var) x_var5 <- sort(unique(var_5$x)) y_var5 <- sort(unique(var_5$y)) z_matrixvar5 <- matrix(krigord_5$krige.var, nrow = length(x_var5), ncol = length(y_var5), byrow = FALSE) pdf("krigging ordinariovar5.pdf") par(mar = c(5, 5, 4, 7.2)) image(x_var5, y_var5, z_matrixvar5, col = viridis::viridis(100), xlab = "x", ylab = "y") image.plot(z = z_matrixvar5, zlim = c(min(z_matrixvar5), max(z_matrixvar5)), col = viridis::viridis(100), legend.lab = "Path loss", horizontal = FALSE, legend.only = TRUE, legend.line = 2) dev.off() #PUNTOS MEDIDOS VS VARIANZA image(x_var5, y_var5, z_matrixvar5, col = viridis::viridis(100), xlab = "x", ylab = "y") points(datos_base2$x, datos_base2$y, pch = 16, cex = 0.4, col = "red") #KRIGGING MEDIA NO CONSTANTE 3. and now i need to do krigging universal with FITML3 R generate-function