setwd("/cloud/project/PROYECTO/DATOS")
datos<-read.csv("petro.csv",
header= TRUE, sep= ",", dec=".")
setwd("/cloud/project/PROYECTO/DATOS/")
datos2<-read.csv("DISCRETAS.csv",
header= TRUE, sep= ",", dec=".", check.names = TRUE)
library(knitr)
library(kableExtra)
## Estadística Inferencial: Variables Nominales y Ordinales
# Tipo de Pozo
well_type <- c(
"Oil Well", "Natural Gas Well", "Oil and Gas Well", "Private Gas Well", "Historical Oil Well",
"Brine Well", "Brine Well Oil Well", "Brine Well Natural Gas Well", "Solution Mining Well", "Source Well",
"Injection Well", "Historical Injection Well", "Natural Gas Storage Well", "Cavern Storage Well",
"Compressed Air Energy Storage Well", "Disposal Well",
"Dry Hole", "Gas Show", "Oil Show", "Oil and Gas Show", "Stratigraphic Test",
"Observation", "Test", "Unknown", "Abandoned", "Shut-in", "Other")
TablaTipo <- data.frame(
well_type = well_type,
ni = sample(1:100, length(well_type)),
hi = round(runif(length(well_type), 0.01, 1), 2))
datos$Fluido <- NA
for (i in 1:nrow(datos)) {
tipo <- as.character(datos$WELL_TYPE[i])
if (tipo %in% c("Oil Well", "Historical Oil Well")) {
datos$Fluido[i] <- "Petróleo"
} else if (tipo %in% c("Natural Gas Well", "Private Gas Well")) {
datos$Fluido[i] <- "Gas Natural"
} else if (tipo %in% c("Oil and Gas Well", "Oil and Gas Show")) {
datos$Fluido[i] <- "Petróleo y Gas"
} else if (tipo %in% c("Brine Well", "Brine Well Oil Well", "Brine Well Natural Gas Well", "Solution Mining Well", "Source Well")) {
datos$Fluido[i] <- "Salmuera/Solución"
} else if (tipo %in% c("Compressed Air Energy Storage Well")) {
datos$Fluido[i] <- "Aire Comprimido"
} else if (tipo %in% c("Injection Well", "Disposal Well", "Historical Injection Well")) {
datos$Fluido[i] <- "Agua/Residuos"
} else if (tipo %in% c("Dry Hole", "Gas Show", "Oil Show", "Stratigraphic Test")) {
datos$Fluido[i] <- "Exploración (Show)"
} else {
datos$Fluido[i] <- "Otro/No definido"
}
}
TDFfluido <- table(datos$Fluido)
Tabla <- data.frame(Tipo_de_Fluido = names(TDFfluido),
Frecuencia_Absoluta = as.vector(TDFfluido))
Tabla$Frecuencia_Relativa <- round(
Tabla$Frecuencia_Absoluta / sum(Tabla$Frecuencia_Absoluta) * 100,
2)
fila_total <- data.frame(
Tipo_de_Fluido = "Total",
Frecuencia_Absoluta = sum(Tabla$Frecuencia_Absoluta),
Frecuencia_Relativa = round(sum(Tabla$Frecuencia_Relativa), 2))
Tabla <- rbind(Tabla, fila_total)
tipo_pozo <- datos$WELL_TYPE
Tabla <- as.data.frame(table(datos$WELL_TYPE))
colnames(Tabla) <- c("Tipo_de_Fluido", "Frecuencia_Absoluta")
print(Tabla)
## Tipo_de_Fluido Frecuencia_Absoluta
## 1 Brine Well 34
## 2 Brine Well Natural Gas Well 4
## 3 Brine Well Oil Well 5
## 4 Cavern Storage Well 124
## 5 Compressed Air Energy Storage Well 2
## 6 Disposal Well 129
## 7 Dry Hole 4546
## 8 Gas Show 2728
## 9 Gas Well Oil Show 56
## 10 Historical Injection Well 49
## 11 Historical Oil Well 1210
## 12 Injection Well 235
## 13 Licensed 2
## 14 Location 936
## 15 Natural Gas Storage Well 333
## 16 Natural Gas Well 9355
## 17 Observation Well 176
## 18 Oil and Gas Show 795
## 19 Oil and Gas Well 253
## 20 Oil Show 977
## 21 Oil Well 2102
## 22 Oil Well Gas Show 53
## 23 Other 19
## 24 Private Gas Well 1057
## 25 Solution Mining Well 204
## 26 Source Well 10
## 27 Stratigraphic Test 1272
TDFfluido <- table(datos$Fluido)
Tabla <- data.frame(Tipo_de_Fluido = names(TDFfluido),
Frecuencia_Absoluta = as.vector(TDFfluido))
Tabla$Frecuencia_Relativa <- round(
Tabla$Frecuencia_Absoluta / sum(Tabla$Frecuencia_Absoluta) * 100,
2)
fila_total <- data.frame(
Tipo_de_Fluido = "Total",
Frecuencia_Absoluta = sum(Tabla$Frecuencia_Absoluta),
Frecuencia_Relativa = round(sum(Tabla$Frecuencia_Relativa), 2))
Tabla <- rbind(Tabla, fila_total)
if (knitr::is_latex_output()) {
kable(Tabla, format = "latex", caption = "Tabla Nº1: Distribución de Tipos de Fluido", booktabs = TRUE)
} else {
kable(Tabla, format = "html", caption = "Tabla Nº1: Distribución de Tipos de Fluido") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE, position = "center") %>%
column_spec(1, bold = TRUE)
}
Tabla Nº1: Distribución de Tipos de Fluido
|
Tipo_de_Fluido
|
Frecuencia_Absoluta
|
Frecuencia_Relativa
|
|
Agua/Residuos
|
413
|
1.55
|
|
Aire Comprimido
|
2
|
0.01
|
|
Exploración (Show)
|
9523
|
35.71
|
|
Gas Natural
|
9355
|
35.08
|
|
Otro/No definido
|
2760
|
10.35
|
|
Petróleo
|
3312
|
12.42
|
|
Petróleo y Gas
|
1048
|
3.93
|
|
Salmuera/Solución
|
253
|
0.95
|
|
Total
|
26666
|
100.00
|
# Probabilidades
Tabla_sin_total <- Tabla[Tabla$Tipo_de_Fluido != "Total", ]
P_TipoFluido <- Tabla_sin_total$Frecuencia_Absoluta / sum(Tabla_sin_total$Frecuencia_Absoluta)
nombres_fluido <- Tabla_sin_total$Tipo_de_Fluido
barplot(P_TipoFluido,
main = "Gráfica: \nDistribución de Probabilidad del Tipo de Fluido",
xlab = "Tipo de Fluido", ylab = "Probabilidad",
col = "darkgreen",
names.arg = nombres_fluido,
las = 3,
cex.names = 0.8)
# Modo del Pozo
Modo <- datos$WELL_MODE
Modo <- na.omit(Modo)
TDFModo <- table(Modo)
TablaModo_EI <- as.data.frame(TDFModo)
fe_Modo <- TablaModo_EI$Freq / sum(TablaModo_EI$Freq)
fe_Modo <- fe_Modo * 100
TablaModo_EI <- data.frame(TablaModo_EI, fe_Modo)
colnames(TablaModo_EI)[colnames(TablaModo_EI) == "fe_Modo"] <- "Frecuencia observada"
TablaModo_final <- TablaModo_EI[, c("Modo", "Frecuencia observada")]
nombres_modo <- TablaModo_final$Modo
P_Modo <- TablaModo_final$`Frecuencia observada`
barplot(P_Modo,
main = "Gráfica: \nDistribución de Probabilidad del Modo de Pozo",
xlab = "Modo de Pozo", ylab = "Probabilidad (%)",
col = "steelblue",
names.arg = nombres_modo,
las = 3,
cex.names = 0.8)
# Concesión
concesiones <- trimws(na.omit(datos$CONCESSION))
es_romano <- function(x) {
grepl("^M{0,4}(CM)?(CD)?(D)?(XC)?(XL)?(L)?(IX)?(IV)?(V)?I{0,3}$", x)
}
tipo_concesion <- ifelse(es_romano(concesiones),
"Numérica (Romana)",
"Alfabética")
tabla_concesion <- table(tipo_concesion)
df_tabla <- as.data.frame(tabla_concesion)
print(colnames(df_tabla))
## [1] "tipo_concesion" "Freq"
print(df_tabla)
## tipo_concesion Freq
## 1 Alfabética 11154
## 2 Numérica (Romana) 15512
colnames(df_tabla) <- c("Tipo_de_Concesion", "Frecuencia_Absoluta")
df_tabla$Frecuencia_Relativa <- round(df_tabla$Frecuencia_Absoluta / sum(df_tabla$Frecuencia_Absoluta) * 100, 2)
print(df_tabla)
## Tipo_de_Concesion Frecuencia_Absoluta Frecuencia_Relativa
## 1 Alfabética 11154 41.83
## 2 Numérica (Romana) 15512 58.17
# ====================
# GRÁFICA FRECUENCIA ABSOLUTA
# ====================
par(mar = c(5, 5, 4, 2))
barplot(df_tabla$Frecuencia_Absoluta,
names.arg = df_tabla$Tipo_de_Concesion,
col = "skyblue",
main = "Frecuencia Absoluta por Tipo de Concesión",
ylab = "Frecuencia",
cex.names = 0.9,
cex.axis = 0.9,
cex.main = 1.2,
ylim = c(0, max(df_tabla$Frecuencia_Absoluta) * 1.2),
border = "black")
# ====================
# GRÁFICA FRECUENCIA RELATIVA
# ====================
par(mar = c(5, 5, 4, 2))
barplot(df_tabla$Frecuencia_Relativa,
names.arg = df_tabla$Tipo_de_Concesion,
col = "lightgreen",
main = "Frecuencia Relativa (%) por Tipo de Concesión",
ylab = "Porcentaje",
cex.names = 0.9,
cex.axis = 0.9,
cex.main = 1.2,
ylim = c(0, 100),
border = "black")
# Tabla de frecuencias absolutas
tabla_concesion <- table(tipo_concesion)
frecuencia <- prop.table(tabla_concesion)
labels <- names(frecuencia)
expected_counts <- rep(1 / length(frecuencia), length(frecuencia))
# Gráfico Observado vs Modelo Uniforme
pos <- barplot(rbind(frecuencia, expected_counts),
beside = TRUE,
col = c("skyblue4", "skyblue"),
names.arg = labels,
ylim = c(0, max(c(frecuencia, expected_counts)) + 0.05),
main = "Comparación Observado vs Modelo Uniforme por Tipo de Concesión",
xlab = "Tipo de Concesión",
ylab = "Frecuencia Relativa (Probabilidad)",
cex.names = 0.9,
las = 1)
legend("topright",
legend = c("Modelo Observado", "Modelo Uniforme"),
fill = c("skyblue4", "skyblue"),
bty = "n",
cex = 0.8)
# Frecuencias absolutas observadas y esperadas
Fo2 <- as.numeric(tabla_concesion)
n_total <- sum(Fo2)
Fe2 <- rep(n_total / length(Fo2), length(Fo2))
# Estadístico Chi-cuadrado
x2_2 <- sum(((Fo2 - Fe2)^2) / Fe2)
k <- length(Fo2)
gl_2 <- k - 1
vc_2 <- qchisq(0.95, df = gl_2)
# Resultados
cat("Chi-cuadrado calculado (x2_2):", x2_2, "\n")
## Chi-cuadrado calculado (x2_2): 712.224
cat("Valor crítico (vc_2):", vc_2, "\n")
## Valor crítico (vc_2): 3.841459
if (x2_2 < vc_2) {
cat("Resultado: No se rechaza H0 (la distribución uniforme se ajusta bien)\n")
} else {
cat("Resultado: Se rechaza H0 (la distribución uniforme no se ajusta bien)\n")
}
## Resultado: Se rechaza H0 (la distribución uniforme no se ajusta bien)
## ESTADÍSTICA INFERENCIAL CONTINUAS
# Variable: Elevación del terreno
elevacion_terreno<-datos$GROUND_ELEVATION
elevacion_terreno<-na.omit(elevacion_terreno)
# Histograma
Histog_elevacion<-hist(elevacion_terreno, freq = FALSE,
main = "Gráfica No.1: Histograma de Elevación del terreno con Curva LogNormal",
xlab = "Elevación del terreno (m))",
ylab = "Densidad de Probabilidad",
col = "darkblue")
# LogNormal
Log_elevacion<-log(elevacion_terreno)
ulog_elevacion<-mean(Log_elevacion)
ulog_elevacion
## [1] 5.275776
sigmalog_elevacion<-sd(Log_elevacion)
sigmalog_elevacion
## [1] 0.1861014
x<-seq(min(elevacion_terreno),max(elevacion_terreno), 0.01)
curve((dlnorm(x, ulog_elevacion,sigmalog_elevacion)),col="red", add=TRUE)
# Área
Altura <- Histog_elevacion$density
Area <- sum(Altura * diff(Histog_elevacion$breaks))
Area
## [1] 1
# Frecuencias observadas y esperadas
FO <- Histog_elevacion$counts
FO
## [1] 50 242 68 14608 8792 544 235 39 27 18 5 0
## [13] 1 1
FOr <- (FO / length(elevacion_terreno)) * 100
sum(FOr)
## [1] 100
Plog<-c()
liminf_elevacion <- Histog_elevacion$breaks[-length(Histog_elevacion$breaks)]
liminf_elevacion
## [1] 0 50 100 150 200 250 300 350 400 450 500 550 600 650
limsup_elevacion <- Histog_elevacion$breaks[-1]
limsup_elevacion
## [1] 50 100 150 200 250 300 350 400 450 500 550 600 650 700
Plog <- plnorm(limsup_elevacion, ulog_elevacion, sigmalog_elevacion) -
plnorm(liminf_elevacion, ulog_elevacion, sigmalog_elevacion)
FE <- Plog*length(elevacion_terreno)
FE
## [1] 2.876750e-09 3.867291e+00 1.895612e+03 1.160279e+04 8.827552e+03
## [6] 2.035949e+03 2.425729e+02 2.018128e+01 1.387887e+00 8.693960e-02
## [11] 5.254642e-03 3.173263e-04 1.956326e-05 1.247666e-06
cor(FO,FE)
## [1] 0.9791728
X2<- sum(((FO - FE)^2)/FE)
X2
## [1] 869037252565
Vc<- qchisq(0.95, df = length(FO) - 1 - 2)
Vc
## [1] 19.67514
cat("¿Acepta H₀ (lognormal)?", X2 < Vc, "\n")
## ¿Acepta H₀ (lognormal)? FALSE
# PROFUNDIDAD ALCANZADA
prof_alcanzada<-datos$TOTAL_DEPTH
prof_alcanzada<-na.omit(prof_alcanzada)
prof_alcanzada<-as.numeric(prof_alcanzada)
Histog_profundidad<-hist(prof_alcanzada, freq = FALSE,
main = "Gráfica No.1: Histograma de Profundidad Alcanzada con Curva Exponencial",
xlab = "Elevación del terreno (m)",
ylab = "Densidad de Probabilidad",
col = "darkblue")
lambda<-1/mean(prof_alcanzada)
x<-seq(min(prof_alcanzada),max(prof_alcanzada),0.01)
curve(dexp(x,rate=lambda),col="red",lwd=2,add=TRUE)
# SUBCONJUNTOS
caja_profundidad<-boxplot(prof_alcanzada, horizontal=T,
main="Diagrama de caja de la profundidad alcanzada",
xlab = "Profundidad Alcanzada")
prof_outliers<-subset(prof_alcanzada,prof_alcanzada>=min(caja_profundidad$out))
profcomunes<-subset(prof_alcanzada,prof_alcanzada<min(caja_profundidad$out))
# Profundidad alcanzada: valores comunes-modelo Lognormal
log_profcomunes <- log(profcomunes)
ulog_profcomunes <- mean(log_profcomunes)
ulog_profcomunes
## [1] 5.580225
sigmalog_profcomunes <- sd(log_profcomunes)
sigmalog_profcomunes
## [1] 0.6883495
Histog_profcomunes<-hist(profcomunes,
main = "Gráfica No.2: Histograma de Profundidad Alcanzada con Curva Lognormal",
xlab = "Profundidad Alcanzada: valores comunes (m)",
col = "darkblue",
ylab = "Densidad de Probabilidad",
freq = FALSE)
x_log <- seq(min(log_profcomunes, na.rm = TRUE), max(log_profcomunes, na.rm = TRUE), length.out = 100)
curve(dlnorm(x, meanlog = ulog_profcomunes, sdlog = sigmalog_profcomunes),
col = "red",
add = TRUE,
lwd = 2)
# Área
Altura <- Histog_profcomunes$density
Area <- sum(Altura * diff(Histog_profcomunes$breaks))
Area
## [1] 1
# Frecuencias observadas y esperadas de profundidad alcanzada: valores comunes
FO_profcomunes <- Histog_profcomunes$counts
Plog_profcomunes<-c()
liminf_profcomunes <- Histog_profcomunes$breaks[-length(Histog_profcomunes$breaks)]
limsup_profcomunes <- Histog_profcomunes$breaks[-1]
Plog_profcomunes <- plnorm(limsup_profcomunes, ulog_profcomunes, sigmalog_profcomunes) -
plnorm(liminf_profcomunes, ulog_profcomunes, sigmalog_profcomunes)
FE_profcomunes<- Plog_profcomunes*length(profcomunes)
cor(FO_profcomunes,FE_profcomunes)
## [1] 0.8863325
X2<- sum(((FO_profcomunes - FE_profcomunes)^2)/FE_profcomunes)
X2
## [1] 3539.131
Vc<- qchisq(0.95, df = length(FO_profcomunes) - 1 - 2)
Vc
## [1] 26.29623
cat("¿Acepta H₀ (lognormal)?", X2 < Vc, "\n")
## ¿Acepta H₀ (lognormal)? FALSE
# Profundidad alcanzada: valores outliers-modelo Lognormal
log_profoutliers <- log(prof_outliers)
ulog_profoutliers <- mean(log_profoutliers)
ulog_profoutliers
## [1] 7.063837
sigmalog_profoutliers <- sd(log_profoutliers)
sigmalog_profoutliers
## [1] 0.2106759
Histog_profoutliers<-hist(prof_outliers)
y_values_for_curve_outliers <- dlnorm(x, meanlog = ulog_profoutliers, sdlog = sigmalog_profoutliers)
max_y_hist_outliers <- max(Histog_profoutliers$density, na.rm = TRUE)
max_y_curve_outliers <- max(y_values_for_curve_outliers, na.rm = TRUE)
ylim_max_outliers <- max(max_y_hist_outliers, max_y_curve_outliers) * 1.05
par(mar = c(7, 6, 5, 4) + 0.1)
Histog_profoutliers<-hist(prof_outliers,
main = "Gráfica No.2:
Histograma de Profundidad Alcanzada con Curva Lognormal",
xlab = "Profundidad Alcanzada: valores comunes (m)",
col = "darkblue",
ylab = "Densidad de Probabilidad",
freq = FALSE,
ylim = c(0, ylim_max_outliers))
x_outliers <- seq(min(log_profoutliers, na.rm = TRUE), max(log_profoutliers, na.rm = TRUE), length.out = 500)
curve(dlnorm(x, meanlog = ulog_profoutliers, sdlog = sigmalog_profoutliers),
col = "red",
add = TRUE,
lwd = 2)
# Área
Altura <- Histog_profoutliers$density
Area <- sum(Altura * diff(Histog_profoutliers$breaks))
Area
## [1] 1
# Frecuencias observadas y esperadas
FO_profoutliers<- Histog_profoutliers$counts
FO_profoutliers <- (Histog_profoutliers$counts/length(prof_outliers))*100
Plog_profoutliers<-c()
liminf_profoutliers <- Histog_profoutliers$breaks[-length(Histog_profoutliers$breaks)]
limsup_profoutliers <- Histog_profoutliers$breaks[-1]
Plog_profoutliers <- plnorm(limsup_profoutliers, ulog_profoutliers, sigmalog_profoutliers) -
plnorm(liminf_profoutliers, ulog_profoutliers, sigmalog_profoutliers)
FE_profoutliers<- Plog_profoutliers*length(prof_outliers)
FE_profoutliers[FE_profoutliers == 0] <- 1e-6
cor(FO_profoutliers,FE_profoutliers)
## [1] 0.9854689
X2_profoutliers<- sum(((FO_profoutliers - FE_profoutliers)^2)/FE_profoutliers)
X2_profoutliers
## [1] 478828.9
Vc_profoutliers<- qchisq(0.95, df = length(FO_profoutliers) - 1 - 2)
Vc_profoutliers
## [1] 12.59159
cat("¿Acepta H₀ (lognormal)?", X2_profoutliers < Vc_profoutliers, "\n")
## ¿Acepta H₀ (lognormal)? FALSE