Poisson Process & NHPP
Najwa Aulia
2026-05-25
Poisson Process
Poisson Process adalah proses stokastik yang menghitung banyaknya kejadian (events) yang terjadi dalam suatu interval waktu, di mana kejadian-kejadian tersebut bersifat acak dan independen satu sama lain. Syarat utamanya: laju rata-rata kejadian (\(\lambda\)) bersifat konstan.
3.2 Simulasi di R
Terdapat dua metode untuk mensimulasikan trajektori Poisson Process:
- Metode 1 — Exponential Interarrivals: Waktu antar-kedatangan dibangkitkan satu per satu dari distribusi eksponensial. Kita menentukan batas jumlah lompatan.
- Metode 2 — Uniform Order Statistics: Berdasarkan teorema bahwa jika \(N(t) = n\), waktu-waktu kejadian terdistribusi seperti order statistics dari \(\text{Uniform}(0, t)\). Kita menentukan batas waktu, bukan jumlah lompatan.
Metode 1: Exponential Interarrivals
Ide dasarnya: setiap waktu antar-kedatangan dibangkitkan dari \(\text{Exp}(\lambda)\) menggunakan transformasi invers, yaitu \(X = -\frac{1}{\lambda} \ln(U)\) dengan \(U \sim \text{Uniform}(0,1)\).
lambda <- 2
njumps <- 20N <- 0:njumps
time <- c()
time[1] <- 0set.seed(333422)for (i in 2:(njumps+1))
time[i] <- time[i-1] + round((-1/lambda) * log(runif(1)), 2)# Setiap tangga horizontal merepresentasikan satu selang waktu tanpa kejadian;
# titik solid = saat kejadian terjadi, lingkaran kosong = batas kanan interval
plot(time, N, type = "n", xlab = "Time", ylab = "State",
panel.first = grid())
segments(time[-length(time)], N[-length(time)],
time[-1] - 0.07, N[-length(time)], lwd = 2, col = "blue")
points(time, N, pch = 20, col = "blue")
points(time[-1], N[-length(time)], pch = 1, col = "blue")
Metode 2: Uniform Order Statistics
Ide dasarnya: jika sudah diketahui ada \(n\) kejadian dalam \([0, t]\), maka posisi waktu kejadian-kejadian itu setara dengan order statistics dari \(n\) bilangan acak \(\text{Uniform}(0, t)\). Langkah pertama adalah membangkitkan \(N(t) \sim \text{Poisson}(\lambda t)\), lalu menyebar waktu-waktu kejadiannya secara acak.
t <- 10
lambda <- 2set.seed(32114)njumps <- rpois(1, lambda * t)N <- 0:njumpsu <- c()
u[1] <- 0for (i in 2:(njumps+1))
u[i] <- runif(1)time <- t * sort(u)# Bandingkan dengan Metode 1: jumlah lompatan di sini bersifat acak (dari rpois),
# sedangkan horizon waktu yang tetap adalah t = 10
plot(time, N, type = "n", xlab = "Time", ylab = "State",
panel.first = grid())
segments(time[-length(time)], N[-length(time)],
time[-1] - 0.07, N[-length(time)], lwd = 2, col = "blue")
points(time, N, pch = 20, col = "blue")
points(time[-1], N[-length(time)], pch = 1, col = "blue")
3.3 Aplikasi Poisson Process
Aplikasi 1: Seismologi (Gempa Bumi)
Tujuan: Membuktikan bahwa data gempa bumi di Southern California (2012–2018) mengikuti Poisson Process menggunakan uji chi-squared goodness-of-fit pada distribusi waktu antar-kedatangan.
Hipotesis:
- \(H_0\): Waktu antar-kedatangan gempa berdistribusi eksponensial (sesuai Poisson Process)
- \(H_1\): Tidak berdistribusi eksponensial
eq.data <- read.table(file = "C:/Users/hi/Downloads/data.txt",
header = TRUE, sep = "", fill = TRUE, comment.char = "")datetime <- as.POSIXct(paste(eq.data[[1]], eq.data[[2]]),
format = "%Y/%m/%d %H:%M:%S", tz = "UTC")datetime <- sort(datetime)
datetime.lag <- c(NA, head(datetime, -1))int.time <- as.numeric(difftime(datetime, datetime.lag, units = "hours"))int.time <- int.time[-1]# Hanya menyimpan interarrival > 3 jam untuk menghilangkan aftershock langsung
int <- int.time[int.time > 3]# Histogram idealnya mendekati distribusi eksponensial (ekor panjang ke kanan)
hist(int,
main = "Interarrival Time Gempa Bumi",
col = "darkmagenta",
xlab = "Interarrival Time (jam)",
breaks = 30)
# Binning interarrival times ke dalam 7 kelas interval
binned.int <- as.factor(ifelse(int < 40, "1",
ifelse(int >= 40 & int < 80, "2",
ifelse(int >= 80 & int < 120, "3",
ifelse(int >= 120 & int < 160, "4",
ifelse(int >= 160 & int < 200, "5",
ifelse(int >= 200 & int < 240, "6", "7")))))))
# Frekuensi observasi per bin
obs <- table(binned.int)
# Estimasi mean untuk distribusi eksponensial
mean.est <- mean(int)
# Frekuensi ekspektasi berdasarkan CDF eksponensial: P(a <= X < b) = e^(-a/mu) - e^(-b/mu)
exp <- c(1:7)
exp[1] <- length(int) * (1 - exp(-40/mean.est))
exp[2] <- length(int) * (exp(-40/mean.est) - exp(-80/mean.est))
exp[3] <- length(int) * (exp(-80/mean.est) - exp(-120/mean.est))
exp[4] <- length(int) * (exp(-120/mean.est) - exp(-160/mean.est))
exp[5] <- length(int) * (exp(-160/mean.est) - exp(-200/mean.est))
exp[6] <- length(int) * (exp(-200/mean.est) - exp(-240/mean.est))
exp[7] <- length(int) * exp(-240/mean.est)
obs## binned.int
## 1 2 3 4 5 6 7
## 339 179 117 49 38 25 42round(exp, 1)## [1] 318.3 189.9 113.3 67.6 40.3 24.0 35.6# Statistik chi-squared: semakin kecil, semakin mirip distribusi observasi dengan ekspektasi
print(chi.sq <- sum((obs - exp)^2 / exp))## [1] 8.534049# df = jumlah bin - 1 - jumlah parameter yang diestimasi = 7 - 1 - 1 = 5
print(p.value <- 1 - pchisq(chi.sq, df = 5))## [1] 0.129156Kesimpulan: Karena \(p\text{-value} > 0.05\), kita gagal tolak \(H_0\). Artinya, tidak ada bukti cukup untuk menyatakan bahwa data menyimpang dari distribusi eksponensial — gempa bumi di Southern California (2012–2018) konsisten dengan Poisson Process.
Aplikasi 2: Sports Analytics (Ice Hockey)
Tujuan: Menghitung probabilitas hasil akhir pertandingan ice hockey 60 menit — seri, Tim A menang, atau Tim B menang — menggunakan deret tak hingga (dipotong di \(n = 15\) untuk konvergensi).
Asumsi: Gol Tim A ~ \(\text{Poisson}(\lambda_A)\) dan Tim B ~ \(\text{Poisson}(\lambda_B)\), keduanya independen. Probabilitas kedua tim mencetak jumlah gol yang sama besar (\(n\)) pada akhir pertandingan dengan total durasi waktu \(T\) dirumuskan secara numerik melalui deret konvergensi:
\[P(\text{Game Ties}) = e^{-(\lambda_A + \lambda_B)T} \sum_{n=0}^{\infty} \frac{(\lambda_A T)^n (\lambda_B T)^n}{(n!)^2}\]
# P(seri): P(N_A = N_B = n) untuk semua n >= 0
# Konvergen setelah ~15 suku karena suku ke-n melibatkan (lambda_A * lambda_B)^n / (n!)^2
sum <- 0
for (n in 0:15)
sum <- sum + 30^n / (factorial(n))^2
sum * exp(-11)## [1] 0.1165575Probabilitas Tim A memenangkan pertandingan dengan selisih gol sebesar \(k\) (\(k \ge 1\)), di mana jumlah gol tim yang kalah (Tim B) adalah \(n\), dihitung menggunakan perulangan bersarang (nested loop):
\[P(\text{Tim A Menang}) = e^{-(\lambda_A + \lambda_B)T} \sum_{n=0}^{\infty} \left[ \frac{(\lambda_A T)^n}{n!} \times \sum_{k=1}^{\infty} \frac{(\lambda_B T)^{n+k}}{(n+k)!} \right]\]
# P(Tim A menang): P(N_A = n+k, N_B = n) untuk n >= 0 dan k >= 1
# Dua loop: n untuk skor dasar Tim B, k untuk keunggulan Tim A
sum.n <- 0
for (n in 0:15) {
sum.k <- 0
for (k in 1:15)
sum.k <- sum.k + 6^k / factorial(n + k)
sum.n <- sum.n + 30^n / factorial(n) * sum.k
}
sum.n * exp(-11)## [1] 0.5589743# P(Tim B menang): P(N_B = n+k, N_A = n) untuk n >= 0 dan k >= 1
# Simetris dengan perhitungan Tim A, tapi menggunakan lambda_B^k
sum.n <- 0
for (n in 0:15) {
sum.k <- 0
for (k in 1:15)
sum.k <- sum.k + 5^k / factorial(n + k)
sum.n <- sum.n + 30^n / factorial(n) * sum.k
}
sum.n * exp(-11)## [1] 0.3244495Catatan: Ketiga probabilitas (seri, A menang, B menang) seharusnya menjumlah ke ~1 sebagai sanity check.
Nonhomogeneous Poisson Process (NHPP)
Pada NHPP, laju \(\lambda(t)\) tidak konstan melainkan bergantung pada waktu. Ini cocok untuk memodelkan fenomena yang intensitasnya berubah-ubah, seperti telepon masuk ke klinik, gempa bumi, atau kebakaran hutan. Didefinisikan melalui:
- Intensity function \(\lambda(t)\): laju sesaat kejadian pada waktu \(t\)
- Integrated rate function \(\Lambda(t) = \int_0^t \lambda(s)\,ds\): ekspektasi kumulatif kejadian hingga waktu \(t\)
4.1 Definisi & Visualisasi
Incoming Calls to Doctor’s Office
Contoh ini mengilustrasikan NHPP dengan laju yang bervariasi sepanjang jam kerja klinik. Grafik pertama menunjukkan laju sesaat, grafik kedua menunjukkan akumulasi kumulatifnya.
# Intensity function: laju panggilan masuk (panggilan/jam) berubah tiap interval waktu
t <- c(0, 1.5, 3, 4, 8)
lambda <- c(10, 10, 5, 8, 4)
plot(t, lambda, type = "n", col = "blue", xlim = c(0, 8),
ylim = c(0, 12), xlab = "Time", ylab = "Intensity rate")
segments(t[-5] + 0.07, lambda[-1], t[-1], lambda[-1], lwd = 2, col = "blue")
points(t, lambda, cex = 1.2, pch = 19, col = "blue")
points(t[-5], lambda[-1], cex = 1.2, pch = 1, col = "blue")
# Integrated rate function Lambda(t): luas di bawah kurva lambda(t)
# Mewakili ekspektasi total panggilan hingga waktu t
t <- c(0, 1.5, 3, 4, 8)
Lambda <- c(0, 15, 22.5, 30.5, 46.5)
plot(t, Lambda, type = "l", lwd = 2, col = "blue", xlim = c(0, 8),
ylim = c(0, 50), xlab = "time", ylab = "integrated rate function")
points(t, Lambda, cex = 1.2, pch = 16, col = "blue")
4.2 Simulasi NHPP
Metode 1: Exponential Interarrivals (NHPP)
Metode ini membangkitkan waktu antar-kedatangan langsung dari NHPP menggunakan inverse transform terhadap \(\Lambda(t)\). Berbeda dengan Poisson homogen, di sini formulanya melibatkan invers dari \(\Lambda(t)\).
# Spesifikasi parameter: 20 lompatan, tanpa batas waktu
njumps <- 20
N <- 0:njumps
time <- c()
set.seed(76855)# Membangkitkan uniform random untuk transformasi invers
u <- c()
for (i in 1:njumps)
u[i] <- runif(1)# Transformasi invers: waktu dihitung dari solusi Lambda(t_{i}) - Lambda(t_{i-1}) = -log(U)
# Untuk lambda(t) = (t+10)/10, inverse Lambda memberikan bentuk sqrt(...)
time[1] <- 0
time[2] <- sqrt(100 - 20 * log(u[1])) - 10
for (i in 3:(njumps + 1)) {
time[i] <- sqrt((time[i-1] + 10)^2 - 20 * log(u[i-1])) - 10
}plot(time, N, type = "n", xlab = "Time", ylab = "State",
panel.first = grid())
segments(time[-length(time)], N[-length(time)],
time[-1] - 0.07, N[-length(time)], lwd = 2, col = "blue")
points(time, N, pch = 20, col = "blue")
points(time[-1], N[-length(time)], pch = 1, col = "blue")
Metode 2: Uniform Order Statistics (NHPP)
Analog dengan Metode 2 di Poisson homogen, hanya saja transformasi dari \(U\) ke waktu kejadian menggunakan invers dari \(\Lambda(t)\) bukan sekadar skala linear.
# t = 10, Lambda(t) = t + 0.05*t^2 (integrated rate)
t <- 10
Lambda <- t + 0.05 * t^2
set.seed(997755)njumps <- rpois(1, Lambda)
N <- 0:njumpsu <- c()
u[1] <- 0
for (i in 2:(njumps + 1))
u[i] <- runif(1)# Invers dari Lambda(t) = t + 0.05t^2 = 10*sqrt(1 + 3*u) - 10
time <- 10 * sqrt(1 + 3 * sort(u)) - 10plot(time, N, type = "n", xlab = "Time", ylab = "State",
panel.first = grid())
segments(time[-length(time)], N[-length(time)],
time[-1] - 0.07, N[-length(time)], lwd = 2, col = "blue")
points(time, N, pch = 20, col = "blue")
points(time[-1], N[-length(time)], pch = 1, col = "blue")
Thinning
Thinning adalah teknik alternatif untuk mensimulasikan NHPP dengan \(\lambda(t)\) yang kompleks. Caranya: bangkitkan dahulu Poisson homogen dengan laju \(\lambda^* \geq \lambda(t)\) untuk semua \(t\) (majorizing rate), lalu terima setiap kejadian dengan probabilitas \(\lambda(t_i)/\lambda^*\).
# lambda(t) = 20 + 20*sin(pi*t): laju berosilasi antara 0 dan 40
# lambda* = 40: batas atas konstan (majorizing rate)
lambda <- function(t) { 20 + 20 * sin(pi * t) }
lambda.star <- function(t) 40
Lambda.star <- function(t) 40 * t
set.seed(2866514)# Langkah 1: bangkitkan Poisson homogen dengan lambda* = 40 di [0,10]
njumps <- rpois(1, Lambda.star(10))
u <- c()
u[1] <- 0
for (i in 2:(njumps + 1))
u[i] <- runif(1)
time.star <- 10 * sort(u)# Langkah 2: thinning — terima kejadian ke-i jika U <= lambda(t_i)/lambda*(t_i)
accepted <- c()
time <- c()
accepted[1] <- 1
time[1] <- 0
for (i in 2:(njumps + 1)) {
if (runif(1) <= lambda(time.star[i]) / lambda.star(time.star[i]))
accepted[i] <- 1 else accepted[i] <- 0
}
time <- time.star[-which(accepted == 0)]
N <- 0:(length(time) - 1)plot(time, N, type = "n", xlab = "Time", ylab = "State",
panel.first = grid())
segments(time[-length(time)], N[-length(time)],
time[-1] - 0.07, N[-length(time)], lwd = 2, col = "blue")
points(time, N, ylim = c(0, 120), pch = 20, col = "blue")
points(time[-1], N[-length(time)], pch = 1, col = "blue")
4.3 Aplikasi NHPP
Aplikasi: Seismologi (NHPP)
Tujuan: Memodelkan laju gempa bumi di Southern California (2010–2020) sebagai NHPP (bukan homogen), karena laju gempa berdenyut seiring waktu. Intensitas \(\lambda(t)\) diestimasi menggunakan regresi polinomial orde 4 pada rata-rata gempa per bulan.
eq.data <- read.csv(file = "C:/Users/hi/Downloads/earthquakedata2010-2020.csv",
header = TRUE, sep = ",")
eq.data$datetime <- as.POSIXct(paste(as.Date(eq.data$DATE), eq.data$TIME))
eq.data$datetime.lag <- c(0, head(eq.data$datetime, -1))
eq.data <- eq.data[-1, ]
eq.data$elapsed.time <- (as.numeric(eq.data$datetime) - as.numeric(eq.data$datetime.lag)) / 3600
eq.data <- eq.data[eq.data$elapsed.time > 1, ]
eq.data$year.month <- format(as.Date(eq.data$DATE), "%Y-%m")
freq.month <- data.frame(table(eq.data$year.month))
year.month.unique <- freq.month[, 1]
neq.month <- freq.month[, 2]
neq.month## [1] 20 20 30 158 110 79 52 42 38 43 29 35 29 19 19 26 18 22
## [19] 24 11 12 16 19 6 10 15 10 16 19 12 20 27 14 11 8 10
## [37] 9 8 8 10 20 14 7 14 14 12 7 12 8 8 17 9 7 11
## [55] 21 9 10 14 5 9 5 10 7 7 14 9 4 2 2 9 5 16
## [73] 10 13 7 8 4 13 6 10 17 11 3 10 6 5 10 11 9 10
## [91] 10 8 8 8 6 9 9 5 6 11 9 7 7 11 11 10 8 7
## [109] 14 11 3 9 5 18 122 35 20 11 25 21 15 12 12 16 19 31
## [127] 10 15 13 23 13 12eq.data <- read.table(file = "C:/Users/hi/Downloads/SearchResults.txt",
header = TRUE, sep = "", fill = TRUE, skip = 10,
comment.char = "")
names(eq.data)[1] <- "DATE"
names(eq.data)[2] <- "TIME"
datetime_string <- paste(eq.data$DATE, eq.data$TIME)
eq.data$datetime <- as.POSIXct(datetime_string, format = "%Y/%m/%d %H:%M:%S")
head(eq.data$datetime)## [1] "2010-02-19 19:54:28 WIB" "2010-02-19 19:56:23 WIB"
## [3] "2010-02-26 15:02:31 WIB" "2010-02-27 07:10:48 WIB"
## [5] "2010-02-27 13:56:38 WIB" "2010-02-27 15:01:35 WIB"eq.data$datetime.lag <- c(0, head(eq.data$datetime, -1))
eq.data <- eq.data[-1, ]
eq.data$elapsed.time <- (as.numeric(eq.data$datetime) - as.numeric(eq.data$datetime.lag)) / 3600
eq.data <- eq.data[eq.data$elapsed.time > 1, ]
eq.data$year.month <- format(as.Date(eq.data$DATE), "%Y-%m")
freq.month <- data.frame(table(eq.data$year.month))
year.month.unique <- freq.month[, 1]
neq.month <- freq.month[, 2]
neq.month## [1] 5 20 102 45 43 15 15 6 11 9 19 8 6 9 8 10 14 9
## [19] 5 7 9 9 2 6 9 8 13 9 4 8 25 5 7 7 6 4
## [37] 5 7 6 13 9 4 11 9 6 5 7 8 5 17 7 6 10 17
## [55] 7 8 8 3 9 5 2 5 5 10 7 4 1 2 5 3 11 7
## [73] 10 5 6 2 12 6 8 15 7 3 7 4 3 8 9 9 10 6
## [91] 7 6 7 4 8 6 4 5 11 8 6 4 6 6 8 5 7 14
## [109] 8 3 5 4 16 121 35 19 11 24 20 10 10 10 14 17 29 10
## [127] 11# Menghapus outlier Juli 2019 (entry ke-107: gempa besar Ridgecrest M7.1)
year.month.unique <- year.month.unique[-107]
neq.month <- neq.month[-107]
library(lubridate)
day1.month <- ymd(paste(year.month.unique, "01", sep = "-"))
library(Hmisc)
ndays.month <- monthDays(as.Date(day1.month, "%Y-%m-%d"))
# Estimasi laju: rata-rata gempa per hari di setiap bulan
lambda <- neq.month / ndays.month
# Waktu median setiap bulan sebagai representasi titik tengah
median.time <- c()
ndays.total <- c()
median.time[1] <- ndays.month[1] / 2
ndays.total[1] <- ndays.month[1]
for (i in 2:length(ndays.month)) {
median.time[i] <- ndays.total[i-1] + ndays.month[i] / 2
ndays.total[i] <- ndays.total[i-1] + ndays.month[i]
}
plot(median.time, lambda)
# Regresi polinomial orde 4 — variabel waktu diskalakan ke ribuan hari (/1000)
# agar koefisien numerik lebih stabil
median.time.re <- median.time / 1000
median.time.sq <- median.time.re^2
median.time.cu <- median.time.re^3
median.time.qd <- median.time.re^4
glm(lambda ~ median.time.re + median.time.sq + median.time.cu + median.time.qd)##
## Call: glm(formula = lambda ~ median.time.re + median.time.sq + median.time.cu +
## median.time.qd)
##
## Coefficients:
## (Intercept) median.time.re median.time.sq median.time.cu median.time.qd
## 1.10409 -2.04328 1.59241 -0.51806 0.06171
##
## Degrees of Freedom: 125 Total (i.e. Null); 121 Residual
## Null Deviance: 28.61
## Residual Deviance: 23.51 AIC: 158# Fungsi fitted lambda: hasil GLM dikodekan langsung sebagai fungsi numerik
plot(median.time, lambda, xlab = "Median Time", ylab = "Lambda")
lambda.fn <- function(t) {
1.11674 - 1.75263*(t/1000) + 1.41060*(t/1000)^2 - 0.49950*(t/1000)^3 + 0.06504*(t/1000)^4
}
lines(median.time, lambda.fn(median.time), lwd = 2, col = "blue")
# Uji goodness-of-fit: data dibagi menjadi 9 bin waktu yang sama-rata ekspektasinya
# Bin ditentukan berdasarkan tanggal kalender agar interval memiliki ekspektasi serupa
time.binned <- as.factor(ifelse(as.Date(eq.data$DATE) < "2011/10/07", "1",
ifelse(as.Date(eq.data$DATE) >= "2011/10/07" & as.Date(eq.data$DATE) < "2012/11/10", "2",
ifelse(as.Date(eq.data$DATE) >= "2012/11/10" & as.Date(eq.data$DATE) < "2013/12/15", "3",
ifelse(as.Date(eq.data$DATE) >= "2013/12/15" & as.Date(eq.data$DATE) < "2015/01/19", "4",
ifelse(as.Date(eq.data$DATE) >= "2015/01/19" & as.Date(eq.data$DATE) < "2016/02/23", "5",
ifelse(as.Date(eq.data$DATE) >= "2016/02/23" & as.Date(eq.data$DATE) < "2017/03/29", "6",
ifelse(as.Date(eq.data$DATE) >= "2017/03/29" & as.Date(eq.data$DATE) < "2018/05/03", "7",
ifelse(as.Date(eq.data$DATE) >= "2018/05/03" & as.Date(eq.data$DATE) < "2019/06/07", "8",
"9")))))))))
obs <- table(time.binned)
exp <- c()
exp[1] <- integrate(lambda.fn, 0, 400)$value
exp[2] <- integrate(lambda.fn, 400, 800)$value
exp[3] <- integrate(lambda.fn, 800, 1200)$value
exp[4] <- integrate(lambda.fn, 1200, 1600)$value
exp[5] <- integrate(lambda.fn, 1600, 2000)$value
exp[6] <- integrate(lambda.fn, 2000, 2400)$value
exp[7] <- integrate(lambda.fn, 2400, 2800)$value
exp[8] <- integrate(lambda.fn, 2800, 3200)$value
exp[9] <- sum(obs) - sum(exp)
obs## time.binned
## 1 2 3 4 5 6 7 8 9
## 367 115 91 113 72 89 94 91 348round(exp, 1)## [1] 333.5 192.9 137.7 123.2 120.7 117.3 116.2 136.7 101.7print(chi.sq <- sum((obs - exp)^2 / exp))## [1] 693.8352# df = jumlah bin - 1 - jumlah parameter GLM yang diestimasi = 9 - 1 - 5 = 3
print(p.value <- 1 - pchisq(chi.sq, df = 3))## [1] 0# Histogram interarrival time di periode awal (sebelum Okt 2011)
int1 <- eq.data$elapsed.time[as.Date(eq.data$DATE) < "2011/10/07"]
hist(int1, main = "", xlab = "", ylab = "", col = "light blue")
# Histogram interarrival time di periode akhir (setelah Jun 2019)
int9 <- eq.data$elapsed.time[as.Date(eq.data$DATE) >= "2019/06/07"]
hist(int9, main = "", xlab = "", ylab = "", col = "light blue")
Exercise 4.2: Wildfire NHPP
Konteks: Diberikan fungsi intensitas kebakaran hutan (wildfires) dalam satu musim selama 120 hari: \[\lambda(t) = -0.000025t^3 + 0.002t^2 + 0.12t\]
a. Intensity Function
Pertanyaan: Plot fungsi intensitas \(\lambda(t)\). Kapan puncak laju kebakaran terjadi?
t <- seq(0, 120, by = 0.1)
lambda_fn <- function(t) { -0.000025*t^3 + 0.002*t^2 + 0.12*t }
# Garis putus-putus hijau menandai puncak intensitas di sekitar t = 74.35 hari
plot(t, lambda_fn(t), type = "l", col = "red", lwd = 2,
xlab = "Days (t)", ylab = "Intensity Rate (lambda)",
main = "Intensity Function lambda(t)")
abline(v = 74.35, col = "darkgreen", lty = 2)
Insight: Intensitas meningkat dari \(t = 0\), mencapai puncak sekitar hari ke-74.35, kemudian menurun — mencerminkan pola musim kebakaran yang intensif di pertengahan dan mereda di akhir musim.
b. Integrated Intensity
Pertanyaan: Plot \(\Lambda(t)\). Berapa rata-rata total kebakaran per musim?
\[\Lambda(t) = \int_0^t \lambda(s)\,ds = -0.00000625t^4 + \frac{0.002}{3}t^3 + 0.06t^2\]
Lambda_fn <- function(t) { -0.00000625*t^4 + (0.002/3)*t^3 + 0.06*t^2 }
plot(t, Lambda_fn(t), type = "l", col = "blue", lwd = 2,
xlab = "Days (t)", ylab = "Expected Cumulative Wildfires (Lambda)",
main = "Integrated Intensity Function Lambda(t)")
avg_total <- Lambda_fn(120)
print(paste("Rata-rata total kebakaran per musim:", avg_total))## [1] "Rata-rata total kebakaran per musim: 720"Insight: \(\Lambda(120)\) memberikan ekspektasi total kejadian selama satu musim penuh. Kurva yang melandai di ujung mencerminkan penurunan intensitas setelah puncak musim.
c. Middle 50% of Season
Pertanyaan: Berapa rata-rata kebakaran yang terjadi di 50% tengah musim (hari ke-30 sampai hari ke-90)?
# Middle 50%: hari ke-30 (kuartil 1) hingga hari ke-90 (kuartil 3)
avg_middle_50 <- Lambda_fn(90) - Lambda_fn(30)
print(paste("Rata-rata kebakaran di middle 50% musim:", round(avg_middle_50, 4)))## [1] "Rata-rata kebakaran di middle 50% musim: 495"Insight: Dengan membandingkan \(\Lambda(90) - \Lambda(30)\) terhadap \(\Lambda(120)\), kita bisa melihat proporsi kebakaran yang terkonsentrasi di tengah musim — konsisten dengan puncak intensitas di sekitar hari ke-74.