# PROSES POISSON & PROSES PEMBAHARUAN
# DATA JUMLAH PELANGGAN PERJAM

# PERSIAPAN
rm(list = ls())

# Panggil package
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
# Import data xl
df <- read_xlsx("D:/sem 4/stokastik/data_chat_masuk_per_jam_7hari.xlsx")

# Lihat struktur data
head(df)
## # A tibble: 6 × 4
##   tanggal    jam_mulai jam_selesai jumlah_chat
##   <chr>      <chr>     <chr>             <dbl>
## 1 2026-01-06 00:00     01:00                 0
## 2 2026-01-06 01:00     02:00                 1
## 3 2026-01-06 02:00     03:00                 2
## 4 2026-01-06 03:00     04:00                 3
## 5 2026-01-06 04:00     05:00                 0
## 6 2026-01-06 05:00     06:00                 0
str(df)
## tibble [168 × 4] (S3: tbl_df/tbl/data.frame)
##  $ tanggal    : chr [1:168] "2026-01-06" "2026-01-06" "2026-01-06" "2026-01-06" ...
##  $ jam_mulai  : chr [1:168] "00:00" "01:00" "02:00" "03:00" ...
##  $ jam_selesai: chr [1:168] "01:00" "02:00" "03:00" "04:00" ...
##  $ jumlah_chat: num [1:168] 0 1 2 3 0 0 4 8 9 7 ...
#pastikan nama kolom sesuai
#1) deskripsi data (minggu 2-4)
n_data <- nrow(df)
min_x <- min(df$jumlah_chat)
max_x <- max(df$jumlah_chat)
mean_x <- mean(df$jumlah_chat)
median_x <- median(df$jumlah_chat)
var_x <- var(df$jumlah_chat)
sd_x <- sd(df$jumlah_chat)

cat("=== DESKRIPSI DATA ===\n")
## === DESKRIPSI DATA ===
cat("Jumlah data (n) =", n_data, "\n")
## Jumlah data (n) = 168
cat("Min =", min_x, "\n")
## Min = 0
cat("Max =", max_x, "\n")
## Max = 23
cat("Mean =", mean_x, "\n")
## Mean = 7.410714
cat("Median =", median_x, "\n")
## Median = 7
cat("Varians =", var_x, "\n")
## Varians = 27.3812
cat("SD =", sd_x, "\n\n")
## SD = 5.232705
#Ringkasan Cepat
summary(df$jumlah_chat)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   3.000   7.000   7.411  11.000  23.000
#histogram (diskrit->binwidht=1)
#ggplot(df,aes(x=jumlah_pelanggan))+
# geom_histogram(aes(y=after_stat(density)),

hist(df$jumlah_chat,main="Histogram Jumlah Chat Masuk", xlab="Jumlah Chat Masuk")

boxplot(df$jumlah_chat,main="Boxplot Jumlah Chat Masuk", ylab="Jumlah Chat Masuk")

# 2) Analisis Probabilitas & Variabel Acak
# Variabel acak: X = jumlah chat masuk per jam (diskrit)
# Estimasi peluang empiris (PMF empiris)
pmf_emp <- df %>%
  count(jumlah_chat) %>%
  mutate(p_emp = n / sum(n)) %>%
  arrange(jumlah_chat)
cat("\n=== PMF EMPIRIS (beberapa baris) ===\n")
## 
## === PMF EMPIRIS (beberapa baris) ===
print(head(pmf_emp, 10))
## # A tibble: 10 × 3
##    jumlah_chat     n  p_emp
##          <dbl> <int>  <dbl>
##  1           0     9 0.0536
##  2           1    10 0.0595
##  3           2    15 0.0893
##  4           3    14 0.0833
##  5           4    10 0.0595
##  6           5    11 0.0655
##  7           6    12 0.0714
##  8           7    14 0.0833
##  9           8    11 0.0655
## 10           9    11 0.0655
# 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 (λ^) = 7.410714
# 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 = 7.410714
cat("Varians =", var_x, "\n")
## Varians = 27.3812
cat("Selisih |Mean - Varians| =", abs(mean_x - var_x), "\n")
## Selisih |Mean - Varians| = 19.97049
# Indeks dispersi: var/mean
dispersion_index <- var_x / mean_x
cat("Dispersion index (Varians/Mean) =", dispersion_index, "\n")
## Dispersion index (Varians/Mean) = 3.694813
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_chat)) +
  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_chat", y = "Probabilitas/Kepadatan")

# 5) Proses Pembaharuan (Renewal) / Waktu Antar Kejadian
cat("\n=== PROSES PEMBAHARUAN (RENEWAL) ===\n")
## 
## === PROSES PEMBAHARUAN (RENEWAL) ===
cat("Data saat ini adalah agregat per jam (jumlah chat masuk per interval waktu),\n")
## Data saat ini adalah agregat per jam (jumlah chat masuk per interval waktu),
cat("bukan waktu kedatangan tiap chat secara individual.\n")
## bukan waktu kedatangan tiap chat secara individual.
cat("Oleh karena itu, waktu antar kedatangan (interarrival time) empiris\n")
## Oleh karena itu, waktu antar kedatangan (interarrival time) empiris
cat("tidak dapat dihitung langsung dari data ini.\n")
## tidak dapat dihitung langsung dari data ini.
cat("Namun secara teori, jika proses Poisson berlaku,\n")
## Namun secara teori, jika proses Poisson berlaku,
cat("maka waktu antar kedatangan chat mengikuti distribusi Eksponensial\n")
## maka waktu antar kedatangan chat mengikuti distribusi Eksponensial
cat("dengan parameter laju (rate) λ = rata-rata jumlah chat per jam.\n")
## dengan parameter laju (rate) λ = rata-rata jumlah chat per jam.
# (Opsional) Simulasi interarrival time teoritis dari
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.1349398 jam
cat("Mean simulasi interarrival =", mean(interarrival_sim), "jam\n")
## Mean simulasi interarrival = 0.136529 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 Chat (Distribusi Eksponensial - Teoritis)",
       x = "Waktu antar chat (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_chat),
            var_per_jam = var(jumlah_chat),
            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: 24 × 4
##    jam_mulai mean_per_jam var_per_jam     n
##    <chr>            <dbl>       <dbl> <int>
##  1 00:00             1.43       1.95      7
##  2 01:00             1.43       1.62      7
##  3 02:00             2          0.667     7
##  4 03:00             1.86       1.14      7
##  5 04:00             1.57       1.62      7
##  6 05:00             1.43       1.29      7
##  7 06:00             7.29       9.57      7
##  8 07:00             7.71       5.24      7
##  9 08:00             7.43       7.62      7
## 10 09:00             6.14       3.48      7
## # ℹ 14 more rows
ggplot(df_hour, aes(x = jam_mulai, y = mean_per_jam)) +
  geom_col() +
  labs(title = "Rata-rata Jumlah Chat Masuk per Jam Mulai",
       x = "Jam Mulai", y = "Rata-rata Jumlah Chat Masuk")