# PERSIAPAN
rm(list = ls())

# Panggil package
library(dplyr)
library(ggplot2)

# Import data CSV
df <- read.csv2("C:/Users/user/OneDrive/Documents/URGENT/csv penting/data_kedatangan_pelanggan_per_jam.csv")
df
##       tanggal jam_mulai jam_selesai jumlah_pelanggan
## 1  06/01/2026     10.00       11.00                9
## 2  06/01/2026     11.00       12.00               33
## 3  06/01/2026     12.00       13.00               26
## 4  06/01/2026     13.00       14.00               17
## 5  07/01/2026     10.00       11.00               10
## 6  07/01/2026     11.00       12.00               19
## 7  07/01/2026     12.00       13.00               21
## 8  07/01/2026     13.00       14.00               14
## 9  08/01/2026     10.00       11.00               12
## 10 08/01/2026     11.00       12.00               24
## 11 08/01/2026     12.00       13.00               25
## 12 08/01/2026     13.00       14.00               16
## 13 09/01/2026     10.00       11.00               13
## 14 09/01/2026     11.00       12.00               14
## 15 09/01/2026     12.00       13.00               35
## 16 09/01/2026     13.00       14.00               16
## 17 10/01/2026     10.00       11.00               14
## 18 10/01/2026     11.00       12.00               23
## 19 10/01/2026     12.00       13.00               23
## 20 10/01/2026     13.00       14.00               22
# Lihat struktur data
head(df)
##      tanggal jam_mulai jam_selesai jumlah_pelanggan
## 1 06/01/2026     10.00       11.00                9
## 2 06/01/2026     11.00       12.00               33
## 3 06/01/2026     12.00       13.00               26
## 4 06/01/2026     13.00       14.00               17
## 5 07/01/2026     10.00       11.00               10
## 6 07/01/2026     11.00       12.00               19
str(df)
## 'data.frame':    20 obs. of  4 variables:
##  $ tanggal         : chr  "06/01/2026" "06/01/2026" "06/01/2026" "06/01/2026" ...
##  $ jam_mulai       : chr  "10.00" "11.00" "12.00" "13.00" ...
##  $ jam_selesai     : chr  "11.00" "12.00" "13.00" "14.00" ...
##  $ jumlah_pelanggan: int  9 33 26 17 10 19 21 14 12 24 ...
#pastikan nama kolom sesuai
#1) deskripsi data (minggu 2-4)
n_data <- nrow(df)
min_x <- min(df$jumlah_pelanggan)
max_x <- max(df$jumlah_pelanggan)
mean_x <- mean(df$jumlah_pelanggan)
median_x<- median(df$jumlah_pelanggan)
var_x <- var(df$jumlah_pelanggan)
sd_x <- sd(df$jumlah_pelanggan)

cat("=== DESKRIPSI DATA ===\n")
## === DESKRIPSI DATA ===
cat("Jumlah data (n) =", n_data, "\n")
## Jumlah data (n) = 20
cat("Min =", min_x, "\n")
## Min = 9
cat("Max =", max_x, "\n")
## Max = 35
cat("Mean =", mean_x, "\n")
## Mean = 19.3
cat("Median =", median_x, "\n")
## Median = 18
cat("Varians =", var_x, "\n")
## Varians = 50.95789
cat("SD =", sd_x, "\n\n")
## SD = 7.13848
#Ringkasan Cepat
summary(df$jumlah_pelanggan)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.00   14.00   18.00   19.30   23.25   35.00
#histogram (diskrit->binwidht=1)
#ggplot(df,aes(x=jumlah_pelanggan))+
# geom_histogram(aes(y=after_stat(density)),

hist(df$jumlah_pelanggan,main="Histogram Jumlah Pelanggan", xlab="jumlah pelanggan")

boxplot(df$jumlah_pelanggan,main="Boxplot jumlah pelanggan", ylab="jumlah pelanggan")

# 3) Estimasi Parameter Distribusi Poisson
# Untuk Poisson, lambda = E[X] -> estimasi pakai rata-rata sampel
lambda_hat <- mean_x
cat("\n=== ESTIMASI PARAMETER ===\n")
## 
## === ESTIMASI PARAMETER ===
cat("Estimasi lambda (λ^) =", lambda_hat, "\n")
## Estimasi lambda (λ^) = 19.3
# 4) Uji Kesesuaian Asumsi Proses Poisson
# 4a) Cek ciri Poisson: mean ~ varians (diagnostik awal)
cat("\n=== CEK CIRI POISSON (diagnostik awal) ===\n")
## 
## === CEK CIRI POISSON (diagnostik awal) ===
cat("Mean =", mean_x, "\n")
## Mean = 19.3
cat("Varians =", var_x, "\n")
## Varians = 50.95789
cat("Selisih |Mean - Varians| =", abs(mean_x - var_x), "\n")
## Selisih |Mean - Varians| = 31.65789
# Indeks dispersi: var/mean
dispersion_index <- var_x / mean_x
cat("Dispersion index (Varians/Mean) =", dispersion_index, "\n")
## Dispersion index (Varians/Mean) = 2.640305
if (dispersion_index > 1.2) {
  cat("Indikasi: OVERDISPERSION (data lebih bervariasi dari
Poisson)\n")
} else if (dispersion_index < 0.8) {
  cat("Indikasi: UNDERDISPERSION\n")
} else {
  cat("Indikasi: Mendekati Poisson (sekitar 1)\n")
}
## Indikasi: OVERDISPERSION (data lebih bervariasi dari
## Poisson)
# 4b) Visual perbandingan data vs Poisson teoretis (OPSIONAL)
df_pois <- data.frame(
  k = 0:max(max_x),
  prob = dpois(0:max(max_x), lambda = lambda_hat)
)
ggplot(df, aes(x = jumlah_pelanggan)) +
  geom_histogram(aes(y = after_stat(density)),
                 binwidth = 1, boundary = -0.5) +
  geom_point(data = df_pois, aes(x = k, y = prob)) +
  geom_line(data = df_pois, aes(x = k, y = prob)) +
  labs(title = "Perbandingan Data dengan Distribusi Poisson Teoretis",
       x = "Jumlah pelanggan", y = "Probabilitas/Kepadatan")

# 5) Proses Pembaharuan (Renewal) / Waktu Antar Kejadian
cat("\n=== PROSES PEMBAHARUAN (RENEWAL) ===\n")
## 
## === PROSES PEMBAHARUAN (RENEWAL) ===
cat("Data saat ini agregat per jam (jumlah pelanggan per interval),\n")
## Data saat ini agregat per jam (jumlah pelanggan per interval),
cat("bukan waktu kedatangan tiap pelanggan.\n")
## bukan waktu kedatangan tiap pelanggan.
cat("Sehingga interarrival time empiris tidak bisa dihitung langsung.\n")
## Sehingga interarrival time empiris tidak bisa dihitung langsung.
cat("Namun secara teori: jika proses Poisson berlaku,\n")
## Namun secara teori: jika proses Poisson berlaku,
cat("waktu antar kedatangan ~ Eksponensial(rate = λ).\n")
## waktu antar kedatangan ~ Eksponensial(rate = λ).
# (Opsional) Simulasi interarrival time teoritis dari Exp(rate=lambda_hat)set.seed(123)
m <- 2000
interarrival_sim <- rexp(m, rate = lambda_hat)  # satuan: "jam" jika λ per jam
cat("\nContoh ringkasan interarrival (simulasi teoritis):\n")
## 
## Contoh ringkasan interarrival (simulasi teoritis):
cat("Mean teoritis interarrival ~ 1/λ =", 1/lambda_hat, "jam\n")
## Mean teoritis interarrival ~ 1/λ = 0.05181347 jam
cat("Mean simulasi interarrival =", mean(interarrival_sim), "jam\n")
## Mean simulasi interarrival = 0.05064549 jam
# Plot distribusi interarrival simulasi (teoritis)
ggplot(data.frame(t = interarrival_sim), aes(x = t)) +  
geom_histogram(aes(y = after_stat(density)), bins = 30) +  
labs(title = "Simulasi Waktu Antar Kedatangan (Eksponensial) - Teoritis", 
x = "Interarrival time (jam)", y = "Kepadatan")

# 6) (Nilai plus) Cek heterogenitas: rata-rata per jam_mulai / per tanggal
# Kalau ada jam tertentu yang lebih ramai, Poisson homogen mungkin gagal.

df_hour <- df %>%  
group_by(jam_mulai) %>%  
summarise(mean_per_jam = mean(jumlah_pelanggan),
var_per_jam = var(jumlah_pelanggan),
n = n(), .groups = "drop") %>% 
arrange(jam_mulai)
cat("\n=== RINGKASAN PER JAM MULAI (cek heterogenitas) ===\n")
## 
## === RINGKASAN PER JAM MULAI (cek heterogenitas) ===
print(df_hour)
## # A tibble: 4 × 4
##   jam_mulai mean_per_jam var_per_jam     n
##   <chr>            <dbl>       <dbl> <int>
## 1 10.00             11.6         4.3     5
## 2 11.00             22.6        49.3     5
## 3 12.00             26          29       5
## 4 13.00             17           9       5
ggplot(df_hour, aes(x = jam_mulai, y = mean_per_jam)) +  
geom_col() +  labs(title = "Rata-rata Jumlah Pelanggan per Jam Mulai", x = "Jam mulai", y = "Mean jumlah pelanggan")