Poisson process

3.2.

3.1. Exponential Interarrival

Simulasi lintasan (trajectory) proses Poisson dengan parameter laju λ=2 dan visualisasi

#specifying parameters
lambda<- 2
njumps<- 20

#defining states
N<- 0:njumps

#setting time as vector
time<- c()

#setting initial value for time
time[1]<- 0

#specifying seed
set.seed(333422)

#simulating trajectory
for (i in 2:(njumps+1))
time[i]<- time[i-1]+round((-1/lambda)*log(runif(1)),2)

Plotting trajectory

# type="n" draws empty frame with no graph
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")

Grafik tersebut menunjukkan hasil simulasi proses Poisson dengan parameter laju λ=2. Grafik berbentuk tangga karena pada proses Poisson nilai proses akan tetap konstan selama tidak ada kejadian, kemudian naik satu unit setiap kali terjadi suatu kejadian baru. Terlihat bahwa proses dimulai dari state 0 pada waktu t=0, lalu meningkat secara bertahap hingga mencapai state 20 pada sekitar waktu t≈12. Jarak antar lompatan tidak sama karena waktu kedatangan antar kejadian bersifat acak dan mengikuti distribusi eksponensial. Pada beberapa interval, seperti sekitar waktu 6 hingga 9, tidak banyak terjadi kejadian sehingga grafik tampak datar lebih panjang. Sebaliknya, di sekitar waktu 9 hingga 12 terjadi banyak lompatan dalam waktu singkat yang menunjukkan kejadian muncul lebih rapat.

3.2. Uniform Order Statistics

#specifying parameters
t<- 10
lambda<- 2

#specifying seed
set.seed(32114)

#generating N(t)
njumps<- rpois(1,lambda*t)

#defining states
N<- 0:njumps

#generating N(t) standard uniforms
u<- c()
u[1]<- 0
for(i in 2:(njumps+1))
u[i]<- runif(1)

#computing event times
time<- t*sort(u)

#plotting trajectory
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")

Grafik adalah realisasi dari suatu proses Poisson di mana, ketika diketahui bahwa terdapat sejumlah n kejadian hingga waktu t, maka waktu terjadinya kejadian-kejadian tersebut tersebar seperti order statistics dari distribusi Uniform(0, t). Hal ini berarti meskipun secara keseluruhan proses mengikuti dinamika acak (random arrivals), jika jumlah total kejadian dikunci, maka posisi waktu kejadian tersebut akan tampak tersebar acak namun merata di sepanjang interval waktu pengamatan. Pada grafik, titik-titik kejadian tersebar di sepanjang sumbu waktu tanpa pola deterministik tertentu, tetapi tetap menunjukkan akumulasi yang meningkat secara bertahap. Variasi jarak antar kejadian (interarrival times) yang terlihat tidak seragam konsisten dengan sifat acak dari distribusi uniform setelah diurutkan, di mana beberapa kejadian bisa tampak berdekatan sementara yang lain berjauhan.

3.3.

3.1. Occurence of Earthquake

library(readr)

## Warning: package 'readr' was built under R version 4.4.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.4.3

## 
## 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

#storing data
cols <- c(
  "DATE", "TIME", "ET", "GT", "MAG", "M",
  "LAT", "LON", "DEPTH", "Q", "EVID", "NPH", "NGRM"
)

eq.data <- read_table(
  "C:/Users/KOKO MEME/Documents/data.txt",
  col_names = cols,
  skip = 1
)

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   DATE = col_date(format = ""),
##   TIME = col_time(format = ""),
##   ET = col_character(),
##   GT = col_character(),
##   MAG = col_double(),
##   M = col_character(),
##   LAT = col_double(),
##   LON = col_double(),
##   DEPTH = col_double(),
##   Q = col_character(),
##   EVID = col_double(),
##   NPH = col_double(),
##   NGRM = col_double()
## )

head(eq.data)

## # A tibble: 6 × 13
##   DATE       TIME     ET    GT      MAG M       LAT   LON DEPTH Q         EVID
##   <date>     <time>   <chr> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <chr>    <dbl>
## 1 2012-01-02 21:11:29 eq    l      3    l      31.7 -116.   8.6 C     11049077
## 2 2012-01-03 14:18:56 eq    l      4.14 l      33.2 -119.  18   B     11049285
## 3 2012-01-06 02:47:51 eq    l      3.08 l      34.0 -116.   4.4 A     11050165
## 4 2012-01-06 06:50:10 eq    l      3.25 l      36.0 -118.   5   A     11050245
## 5 2012-01-06 16:58:54 eq    l      3.47 l      33.7 -116.   4.5 A     11050541
## 6 2012-01-08 01:59:54 eq    l      3.51 l      33.7 -116.   6.4 A     11051197
## # ℹ 2 more variables: NPH <dbl>, NGRM <dbl>

#creating date-time variable
datetime<- as.POSIXct(paste(as.Date(eq.data$DATE),
eq.data$TIME))

#computing lag
datetime.lag<- c(0,head(datetime, -1))

#computing interarrival times (in hours)
int.time<- (as.numeric(datetime)-as.numeric(datetime.lag))/
3600

#removing first value
int.time<- int.time[-1]

#removing immediate aftershocks (within 3 hours)
int<- int.time[int.time>3]

#plotting histogram
hist(int, main="", col="dark magenta", xlab="Interarrival
Time")

Histogram menggambarkan distribusi interarrival times (jarak waktu antar gempa) setelah data dibersihkan dari kemungkinan aftershocks. Secara visual, distribusi tersebut tampak miring ke kanan (right-skewed), di mana sebagian besar waktu antar kejadian bernilai kecil (dekat nol), dan frekuensinya menurun secara bertahap seiring bertambahnya waktu. Pola ini merupakan ciri khas dari distribusi eksponensial. Hal ini menunjukkan bahwa gempa bumi cenderung terjadi dengan interval waktu yang pendek lebih sering, sementara interval yang panjang relatif jarang terjadi. Ini konsisten dengan sifat memoryless dari distribusi eksponensial, di mana peluang terjadinya gempa berikutnya tidak bergantung pada kapan gempa sebelumnya terjadi. Dengan kata lain, meskipun beberapa interval panjang tetap muncul (terlihat dari ekor panjang di kanan histogram), mayoritas kejadian terjadi dalam waktu yang relatif dekat satu sama lain.

Uji Chi-Square

#binning interarrival times
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")))))))

#computing observed frequencies
obs<- table(binned.int)

#estimating mean for exponential distribution
mean.est<- mean(int)

#computing expected frequencies
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  42

#rounding
round(exp,1)

## [1] 318.3 189.9 113.3  67.6  40.3  24.0  35.6

#computing chi-squared statistic
chi.sq<- sum((obs-exp)^2/exp)
print(chi.sq)

## [1] 8.534049

#computing p-value
print(p.value<- 1-pchisq(chi.sq, df=5))

## [1] 0.129156

Berdasarkan hasil uji chi-squared yang diperoleh, nilai p-value lebih besar dari tingkat signifikansi umum (α = 0.05), maka tidak terdapat cukup bukti untuk menolak hipotesis nol. Dengan demikian, dapat disimpulkan bahwa distribusi interarrival time tidak berbeda secara signifikan dari distribusi eksponensial.

3.2. Scoring Goals

a

lambda_A <- 1/10
lambda_B <- 1/12

lambda_total <- lambda_A + lambda_B
expected_wait <- 1 / lambda_total
expected_wait

## [1] 5.454545

b

p_A_first <- lambda_A / (lambda_A + lambda_B)
p_B_first <- lambda_B / (lambda_A + lambda_B)

p_A_first

## [1] 0.5454545

p_B_first

## [1] 0.4545455

c

# Probability Ties
sum<- 0
for(n in 0:15){
  sum<- sum+(30^n)/(factorial(n))^2
}
sum*exp(-11)

## [1] 0.1165575

# Probability A wins
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] 7.077008

# Probability B wins
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] 4.324333

NHPP Earthquake

eq.data <- read.csv("C:/Users/KOKO MEME/Downloads/earthquakedata2010-2020 (1).csv")
tail(eq.data)

##            DATE        TIME  MAG      LAT       LON DEPTH
## 4386 2020-12-18 16:41:35.75 3.24 31.51083 -115.7487  13.8
## 4387 2020-12-22 00:30:00.24 4.05 32.33717 -115.3908  15.2
## 4388 2020-12-23 22:22:19.69 3.34 32.33583 -115.3795  15.2
## 4389 2020-12-25 20:25:59.07 3.64 33.74733 -115.9212   3.7
## 4390 2020-12-27 04:05:37.24 3.50 35.48867 -118.2883   2.7
## 4391 2020-12-28 14:34:53.60 3.37 32.43850 -115.1648  27.8

head(eq.data)

##         DATE        TIME  MAG      LAT       LON DEPTH
## 1 2010-01-01 02:33:42.82 3.24 32.45383 -115.2028   6.0
## 2 2010-01-04 00:02:17.99 3.81 32.19150 -116.6550   5.6
## 3 2010-01-04 20:07:59.93 3.39 35.71167 -117.4822   8.6
## 4 2010-01-04 20:11:25.79 3.12 35.71150 -117.4850   9.2
## 5 2010-01-10 02:55:06.75 3.04 36.40133 -117.8687   3.9
## 6 2010-01-10 06:35:53.56 4.08 32.40333 -115.1755   6.0

print(nrow(eq.data))

## [1] 4391

library(lubridate)

## Warning: package 'lubridate' was built under R version 4.4.3

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# Make Datetime
eq.data$datetime <- ymd_hms(paste(eq.data$DATE, eq.data$TIME), tz = "UTC")

# Computing lag
eq.data$datetime.lag <- c(NA, head(eq.data$datetime, -1))
class(eq.data$datetime.lag) <- class(eq.data$datetime)

# Removing first row
eq.data <- eq.data[-1, ]

print(sum(is.na(eq.data$datetime)))

## [1] 0

#computing interarrival times (in hours)
eq.data$elapsed.time<- (as.numeric(eq.data$datetime)
-as.numeric(eq.data$datetime.lag))/3600
#removing immediate aftershocks (within 1 hour)
eq.data<- eq.data[eq.data$elapsed.time>1,]

eq.data$year.month <- format(eq.data$datetime, "%Y-%m")

freq.month <- as.data.frame(table(eq.data$year.month))
colnames(freq.month) <- c("year.month", "freq")

year.month.unique <- freq.month$year.month
neq.month <- freq.month$freq
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  12

#removing outlier (July, 2019, 107th entry)
year.month.unique<- year.month.unique[-107]
neq.month<- neq.month[-107]
#computing number of days per month
library(lubridate)
day1.month <- ymd(paste(year.month.unique,"01", sep="-") )
library(Hmisc)

## Warning: package 'Hmisc' was built under R version 4.4.3

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

ndays.month<- monthDays(as.Date(day1.month, "%Y-%m-%d"))
#estimating intensity rate of earthquakes per day
lambda<- neq.month/ndays.month
#computing cumulative number of days until median day of each
month

## function (x, label = FALSE, abbr = TRUE, locale = Sys.getlocale("LC_TIME")) 
## {
##     UseMethod("month")
## }
## <bytecode: 0x0000027a277aa7e8>
## <environment: namespace:lubridate>

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]
}
#plotting lambda against median time
plot(median.time, lambda)

#regressing lambda on time
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  
##        2.16785        -3.38749         2.11080        -0.55429         0.05406  
## 
## Degrees of Freedom: 130 Total (i.e. Null);  126 Residual
## Null Deviance:       58.57 
## Residual Deviance: 37.3  AIC: 219.2

#adding fitted line
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 }
plot(median.time, lambda)
lines(median.time, lambda.fn(median.time), lwd=2,col="blue")

Grafik ini menampilkan estimasi intensity rate λ harian gempa bumi di Southern California periode 2010–2020, di mana setiap titik merepresentasikan λ bulanan yang dihitung sebagai rasio jumlah gempa per bulan terhadap jumlah hari dalam bulan tersebut, dan kurva biru adalah hasil fit polynomial degree 4. Pola yang paling mencolok adalah lonjakan λ yang sangat tinggi di awal periode (t = 0–300 hari, λ mencapai 5.2) yang mencerminkan efek aftershock dari gempa besar di awal periode pengamatan, kemudian turun drastis dan stabil di nilai rendah sekitar 0.2–0.4 pada t = 1000–3000 hari. Kurva polynomial degree 4 berhasil menangkap pola umum ini dengan bentuk U yang asimetris meskipun masih terdapat beberapa titik yang jauh dari kurva terutama di awal dan satu outlier sekitar t = 3400. Secara keseluruhan, grafik ini mengkonfirmasi bahwa rate gempa tidak konstan sepanjang waktu, sehingga pemodelan NHPP lebih tepat dibandingkan Homogeneous Poisson Process.

goodness-of-fit test

Hipotesis:
H₀: Data gempa bumi di Southern California tahun 2010–2020 mengikuti NHPP dengan intensity function λ̂(t) = 1.11674 − 1.75263(t/1000) + 1.41060(t/1000)² − 0.49950(t/1000)³ + 0.06504(t/1000)⁴
H₁: Data gempa bumi tidak mengikuti model NHPP dengan intensity function tersebut

Tolak H₀ apabila p-value > 0.05

#binning times
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")))))))))
#computing observed frequencies
obs<- table((time.binned))

#computing expected frequencies
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

## 
##   1   2   3   4   5   6   7   8   9 
## 837 197 139 141 108 112 113 119 434

round(exp,1)

## [1] 333.5 192.9 137.7 123.2 120.7 117.3 116.2 136.7 921.7

#computing chi-squared statistic
print(chi.sq<- sum((obs-exp)^2/exp))

## [1] 1024.758

#computing p-value
print(p.value<- 1-pchisq(chi.sq, df=3))

## [1] 0

Berdasarkan hasil uji goodness-of-fit, diperoleh nilai statistik χ² = 1024.758 dengan df = 3 dan p-value = 0. Karena p-value < 0.05, maka H₀ ditolak, yang berarti data gempa bumi di Southern California periode 2010–2020 tidak mengikuti model NHPP dengan intensity function polynomial degree 4 yang diestimasi.

Histogram first vs last bins

int1<- eq.data$elapsed.time[as.Date(eq.data$DATE)
<"2011/10/07"]
hist(int1, main="", xlab="", ylab="", col="light blue")

int9<- eq.data$elapsed.time[as.Date(eq.data$DATE)
>="2019/06/07"]
hist(int9, main="", xlab="", ylab="", col="light blue")

Kedua grafik ini mencerminkan histogram interarrival times dari data gempa bumi Southern California, yang digunakan untuk memverifikasi secara visual bahwa waktu antar-kedatangan dalam setiap bin mengikuti distribusi eksponensial.
Histogram pertama (bin ke-1, t = 0–400 hari) menunjukkan distribusi yang sangat right-skewed dengan puncak tinggi di interval kecil (0–25 jam) mencapai frekuensi ~550, lalu menurun tajam dan mendekati nol setelah 100 jam. Bentuk ini konsisten dengan distribusi eksponensial yang memiliki rate tinggi.
Histogram kedua (bin ke-9, t = 3200–3647 hari) juga menunjukkan bentuk right-skewed yang serupa dengan distribusi eksponensial, namun dengan puncak yang lebih rendah (~345) dan ekor yang lebih panjang hingga 300 jam.
Secara keseluruhan, kedua histogram mendukung asumsi bahwa interarrival times dalam masing-masing bin berdistribusi eksponensial₀.

Exercise 4.2

# Define functions
lambda <- function(t) -0.000025 * t^3 + 0.002 * t^2 + 0.12 * t

Lambda <- function(t) integrate(lambda, 0, t)$value

t <- seq(0, 120, length.out = 1000)

a

# Plot intensity function
plot(t, lambda(t), type = "l", col = "blue", lwd = 2,
     main = "Intensity Function λ(t)",
     xlab = "Time (days)", ylab = "λ(t)")
abline(h = 0, col = "gray", lty = 2)

# Peak of intensity
t_peak <- optimize(lambda, interval = c(0, 120), maximum = TRUE)
cat("Peak intensity at t =", round(t_peak$maximum, 4), 
    "with λ(t) =", round(t_peak$objective, 4), "\n")

## Peak intensity at t = 74.7407 with λ(t) = 9.7034

b

# Plot integrated intensity function
Lambda_vec <- sapply(t, Lambda)

plot(t, Lambda_vec, type = "l", col = "red", lwd = 2,
     main = "Integrated Intensity Function Λ(t)",
     xlab = "Time (days)", ylab = "Λ(t)")

# Average number of wildfires per season (full 120 days)
avg_total <- Lambda(120)
cat("Average number of wildfires per season:", round(avg_total, 4), "\n")

## Average number of wildfires per season: 720

c

# Middle 50% of season = t in [30, 90]
avg_middle <- Lambda(90) - Lambda(30)
cat("Average number of wildfires in middle 50% (t=30 to t=90):", 
    round(avg_middle, 4), "\n")

## Average number of wildfires in middle 50% (t=30 to t=90): 495