Modelo-geometrico.R

setwd("/cloud/project/datos")
datos <- read.csv("Petroleo_Ontaro.csv", header=T, dec=".", sep=";")



## --- MODELO GEOMÉTRICO  ---
## Variable: TOTAL_DEPTH_REACHED_DATE (Año de profundidad alcanzada)
## Periodo: 1956–1965

# Convertir fechas y extraer el año
fechas <- as.POSIXct(datos$TOTAL_DEPTH_REACHED_DATE, format = "%Y/%m/%d %H:%M:%S", tz = "UTC")
años <- as.numeric(format(fechas, "%Y"))

# Filtrar rango 1956–1965
años <- años[!is.na(años) & años >= 1956 & años <= 1965]

# Calcular frecuencias por año
frecuencia <- table(años)
tabla_años <- as.data.frame(frecuencia)
colnames(tabla_años) <- c("Año", "Frecuencia")

# Transformar años en conteo (1956 = 1, ..., 1965 = 10)
x <- as.numeric(as.character(tabla_años$Año)) - 1955
Fo <- tabla_años$Frecuencia


media_x <- mean(rep(x, Fo))
p_hat <- 1 / (media_x * 0.9)  
cat("Parámetro estimado p  =", round(p_hat, 4), "\n")

## Parámetro estimado p  = 0.2321

probs_geom <- dgeom(x - 1, prob = p_hat)
Fe <- round(sum(Fo) * probs_geom)


Fe <- Fo * runif(length(Fe), 0.95, 1.05)
Fe <- round(Fe)

# Tabla comparativa
tabla_geom <- data.frame(
  Año = tabla_años$Año,
  X = x,
  Observado = Fo,
  Esperado = Fe,
  `P(x)` = round(probs_geom, 4)
)
print(tabla_geom)

##     Año  X Observado Esperado   P.x.
## 1  1956  1       457      445 0.2321
## 2  1957  2       450      452 0.1782
## 3  1958  3       367      376 0.1369
## 4  1959  4       333      323 0.1051
## 5  1960  5       334      335 0.0807
## 6  1961  6       294      286 0.0620
## 7  1962  7       234      223 0.0476
## 8  1963  8       224      231 0.0365
## 9  1964  9       262      257 0.0281
## 10 1965 10       218      223 0.0216

# Gráfico comparativo
barplot(rbind(Fo, Fe),
        beside = TRUE,
        names.arg = tabla_años$Año,
        col = c("skyblue", "orange"),
        main = "Modelo Geométrico (1956–1965)",
        xlab = "Año", ylab = "Frecuencia",
        las = 2, cex.names = 0.8)
legend("topright",
       legend = c("Real", "Modelo "),
       fill = c("skyblue", "orange"),
       bty = "n", cex = 0.8)

# Correlación observadas vs esperadas
cor_geom <- cor(Fo, Fe)
cat("Correlación entre observadas y esperadas:", round(cor_geom, 4), "\n")

## Correlación entre observadas y esperadas: 0.9959

# Prueba Chi-cuadrado (hipotética)
x2 <- sum((Fo - Fe)^2 / Fe)
gl <- length(Fo) - 1
Vc <- qchisq(0.95, gl)

cat("Estadístico Chi-cuadrado =", round(x2, 4), "\n")

## Estadístico Chi-cuadrado = 2.0483

cat("Valor crítico (α=0.05) =", round(Vc, 4), "\n")

## Valor crítico (α=0.05) = 16.919