# Sintaks Project mandiri PROSES STOKASTIK
# BAGIAN D - SIMULASI KOMPUTASI
# 1. LIBRARY
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(igraph)
## Warning: package 'igraph' was built under R version 4.4.3
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
# 2. PARAMETER SIMULASI
# Jumlah hari simulasi
n <- 1000
# Kondisi awal
S0 <- 0
# Matriks probabilitas transisi
P <- matrix(c(
0.65, 0.30, 0.05,
0.25, 0.55, 0.20,
0.10, 0.40, 0.50
), nrow = 3, byrow = TRUE)
rownames(P) <- c("S0", "S1", "S2")
colnames(P) <- c("S0", "S1", "S2")
# Menampilkan matriks transisi
P
## S0 S1 S2
## S0 0.65 0.30 0.05
## S1 0.25 0.55 0.20
## S2 0.10 0.40 0.50
# 3. SIMULASI KONDISI CUACA
# Seed agar hasil simulasi dapat direproduksi
set.seed(123)
# Menyimpan kondisi cuaca
S <- numeric(n)
# Kondisi awal
S[1] <- S0
# Membangkitkan kondisi hari berikutnya
for (t in 1:(n - 1)) {
S[t + 1] <- sample(
x = 0:2,
size = 1,
prob = P[S[t] + 1, ]
)
}
# Melihat 20 kondisi pertama
head(S, 20)
## [1] 0 0 1 1 2 0 0 0 1 0 0 2 2 1 0 0 1 1 1 1
# 4. SIMULASI CURAH HUJAN
R <- numeric(n)
for (t in 1:n) {
if (S[t] == 0) {
# S0 : Uniform(0,5)
R[t] <- runif(
1,
min = 0,
max = 5
)
} else if (S[t] == 1) {
# S1 : Exponential(lambda = 1/20)
R[t] <- rexp(
1,
rate = 1 / 20
)
} else {
# S2 : Gamma(k = 2, theta = 30)
R[t] <- rgamma(
1,
shape = 2,
scale = 30
)
}
}
# 5. MEMBUAT DATA FRAME
data_simulasi <- data.frame(
Hari = 1:n,
Kondisi = S,
Curah_Hujan = R
)
# Label kondisi
data_simulasi$Kondisi_Label <- factor(
data_simulasi$Kondisi,
levels = c(0, 1, 2),
labels = c(
"S0 - Cerah/Berawan",
"S1 - Hujan Sedang",
"S2 - Hujan Lebat/Ekstrem"
)
)
# Melihat 10 data pertama
head(data_simulasi, 10)
## Hari Kondisi Curah_Hujan Kondisi_Label
## 1 1 0 0.5441204 S0 - Cerah/Berawan
## 2 2 0 1.3681137 S0 - Cerah/Berawan
## 3 3 1 3.7546773 S1 - Hujan Sedang
## 4 4 1 33.3554574 S1 - Hujan Sedang
## 5 5 2 91.9143284 S2 - Hujan Lebat/Ekstrem
## 6 6 0 2.3894341 S0 - Cerah/Berawan
## 7 7 0 3.8684606 S0 - Cerah/Berawan
## 8 8 0 1.4770004 S0 - Cerah/Berawan
## 9 9 1 42.5898266 S1 - Hujan Sedang
## 10 10 0 2.2026656 S0 - Cerah/Berawan
# 6. JUMLAH DAN PROPORSI KONDISI CUACA
jumlah_kondisi <- table(
data_simulasi$Kondisi
)
proporsi_kondisi <- prop.table(
jumlah_kondisi
)
persentase_kondisi <- proporsi_kondisi * 100
jumlah_kondisi
##
## 0 1 2
## 340 433 227
proporsi_kondisi
##
## 0 1 2
## 0.340 0.433 0.227
persentase_kondisi
##
## 0 1 2
## 34.0 43.3 22.7
# Membuat tabel hasil
tabel_kondisi <- data.frame(
Kondisi = c(
"Cerah/Berawan",
"Hujan Sedang",
"Hujan Lebat/Ekstrem"
),
Jumlah_Hari = as.numeric(
jumlah_kondisi
),
Proporsi = as.numeric(
proporsi_kondisi
),
Persentase = as.numeric(
persentase_kondisi
)
)
tabel_kondisi
## Kondisi Jumlah_Hari Proporsi Persentase
## 1 Cerah/Berawan 340 0.340 34.0
## 2 Hujan Sedang 433 0.433 43.3
## 3 Hujan Lebat/Ekstrem 227 0.227 22.7
# VISUALISASI 1
# DISTRIBUSI KONDISI CUACA
# Warna kondisi cuaca
warna_state <- c(
"S0 - Cerah/Berawan" = "#F4B942",
"S1 - Hujan Sedang" = "#4C9BD3",
"S2 - Hujan Lebat/Ekstrem" = "#263F66"
)
# Urutan kondisi
tabel_kondisi$Kondisi <- factor(
tabel_kondisi$Kondisi,
levels = c(
"Cerah/Berawan",
"Hujan Sedang",
"Hujan Lebat/Ekstrem"
)
)
ggplot(
tabel_kondisi,
aes(
x = Kondisi,
y = Persentase,
fill = Kondisi
)
) +
geom_col(
width = 0.65,
color = "white",
linewidth = 0.5
) +
geom_text(
aes(
label = paste0(
round(Persentase, 1),
"%"
)
),
vjust = -0.5,
fontface = "bold",
size = 4
) +
scale_fill_manual(
values = c(
"Cerah/Berawan" = "#F4B942",
"Hujan Sedang" = "#4C9BD3",
"Hujan Lebat/Ekstrem" = "#263F66"
)
) +
labs(
title = "Distribusi Kondisi Cuaca Hasil Simulasi",
x = "Kondisi Cuaca",
y = "Persentase (%)"
) +
scale_y_continuous(
expand = expansion(
mult = c(0, 0.10)
)
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
axis.title = element_text(
face = "bold"
),
axis.text.x = element_text(
angle = 0,
hjust = 0.5
),
legend.position = "none",
panel.grid.minor = element_blank()
)
# 7. DIAGRAM TRANSISI CUACA
# 7.1 DATA TRANSISI
edges <- data.frame(
from = c(
"S0", "S0", "S0",
"S1", "S1", "S1",
"S2", "S2", "S2"
),
to = c(
"S0", "S1", "S2",
"S0", "S1", "S2",
"S0", "S1", "S2"
),
prob = c(
P[1,1], P[1,2], P[1,3],
P[2,1], P[2,2], P[2,3],
P[3,1], P[3,2], P[3,3]
)
)
# Label probabilitas
edges$label <- sprintf(
"%.2f",
edges$prob
)
# 7.2 DATA NODE
vertices <- data.frame(
name = c(
"S0",
"S1",
"S2"
),
label = c(
"S0\nCerah/Berawan",
"S1\nHujan Sedang",
"S2\nHujan Lebat/Ekstrem"
)
)
# 7.3 MEMBUAT GRAPH
g <- graph_from_data_frame(
d = edges,
directed = TRUE,
vertices = vertices
)
# 7.4 WARNA NODE
warna_node <- c(
S0 = "#F4B942",
S1 = "#4C9BD3",
S2 = "#263F66"
)
V(g)$color <- warna_node[
V(g)$name
]
V(g)$size <- 8
V(g)$frame.color <- "white"
V(g)$frame.width <- 2
V(g)$label.color <- "white"
V(g)$label.cex <- 1
# 7.5 PENGATURAN PANAH
E(g)$color <- "#64748B"
E(g)$width <- 2
E(g)$arrow.size <- 0.35
E(g)$label <- edges$label
E(g)$label.color <- "black"
E(g)$label.cex <- 2
# 7.6 MEMBUAT PANAH TIDAK BERTUMPUK
E(g)$curved <- c(
# S0 -> S0
0.6,
# S0 -> S1
0.18,
# S0 -> S2
-0.18,
# S1 -> S0
-0.18,
# S1 -> S1
0.6,
# S1 -> S2
0.18,
# S2 -> S0
0.18,
# S2 -> S1
-0.18,
# S2 -> S2
0.6
)
# 7.7 POSISI NODE
layout_cuaca <- matrix(
c(
-1.2, -0.8, # S0
1.2, -0.8, # S1
0.0, 1.0 # S2
),
ncol = 2,
byrow = TRUE
)
# 7.8 MENAMPILKAN DIAGRAM
plot(
g,
layout = layout_cuaca,
# Node
vertex.color = V(g)$color,
vertex.size = V(g)$size,
vertex.label = V(g)$label,
vertex.label.color = V(g)$label.color,
vertex.label.cex = V(g)$label.cex,
vertex.frame.color = V(g)$frame.color,
vertex.frame.width = V(g)$frame.width,
# Panah
edge.color = E(g)$color,
edge.width = E(g)$width,
edge.arrow.size = E(g)$arrow.size,
edge.curved = E(g)$curved,
# Probabilitas
edge.label = E(g)$label,
edge.label.color = E(g)$label.color,
edge.label.cex = E(g)$label.cex,
# Tampilan
rescale = TRUE,
asp = 0.90,
margin = 0.88
)
# 8. RATA-RATA CURAH HUJAN SETIAP KONDISI
mean_per_kondisi <- tapply(
data_simulasi$Curah_Hujan,
data_simulasi$Kondisi,
mean
)
mean_per_kondisi
## 0 1 2
## 2.618053 20.586376 58.026557
# Tabel rata-rata
tabel_mean <- data.frame(
Kondisi = c(
"Cerah/Berawan",
"Hujan Sedang",
"Hujan Lebat/Ekstrem"
),
Rata_rata_Simulasi = as.numeric(
mean_per_kondisi
),
Rata_rata_Teoritis = c(
2.5,
20,
60
)
)
tabel_mean
## Kondisi Rata_rata_Simulasi Rata_rata_Teoritis
## 1 Cerah/Berawan 2.618053 2.5
## 2 Hujan Sedang 20.586376 20.0
## 3 Hujan Lebat/Ekstrem 58.026557 60.0
# VISUALISASI 2
# PERBANDINGAN RATA-RATA CURAH HUJAN
# Data untuk grafik
data_mean_plot <- data.frame(
Kondisi = rep(
c(
"S0",
"S1",
"S2"
),
each = 2
),
Label = rep(
c(
"Cerah/Berawan",
"Hujan Sedang",
"Hujan Lebat/Ekstrem"
),
each = 2
),
Metode = rep(
c(
"Teoritis",
"Simulasi"
),
times = 3
),
Nilai = c(
2.5,
as.numeric(mean_per_kondisi)[1],
20,
as.numeric(mean_per_kondisi)[2],
60,
as.numeric(mean_per_kondisi)[3]
)
)
ggplot(
data_mean_plot,
aes(
x = Kondisi,
y = Nilai,
fill = Metode
)
) +
geom_col(
position = position_dodge(
width = 0.75
),
width = 0.62
) +
geom_text(
aes(
label = round(
Nilai,
2
)
),
position = position_dodge(
width = 0.75
),
vjust = -0.4,
fontface = "bold",
size = 3.5
) +
scale_fill_manual(
values = c(
"Teoritis" = "#263F66",
"Simulasi" = "#F4B942"
)
) +
scale_x_discrete(
labels = c(
"S0" = "S0\nCerah/Berawan",
"S1" = "S1\nHujan Sedang",
"S2" = "S2\nHujan Lebat"
)
) +
labs(
title = "Perbandingan Rata-Rata Curah Hujan",
x = "Kondisi Cuaca",
y = "Rata-rata Curah Hujan (mm)",
fill = "Metode"
) +
scale_y_continuous(
expand = expansion(
mult = c(0, 0.12)
)
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
axis.title = element_text(
face = "bold"
),
legend.position = "bottom",
panel.grid.minor = element_blank()
)
# 9. DATA DISTRIBUSI TEORITIS
# S0 - Uniform(0,5)
x0 <- seq(
0,
5,
length.out = 300
)
teori_s0 <- data.frame(
Curah_Hujan = x0,
Density = dunif(
x0,
min = 0,
max = 5
),
Kondisi_Label = "S0 - Cerah/Berawan"
)
# S1 - Exponential(lambda = 1/20)
x1 <- seq(
0,
quantile(
data_simulasi$Curah_Hujan[
data_simulasi$Kondisi == 1
],
0.99
),
length.out = 300
)
teori_s1 <- data.frame(
Curah_Hujan = x1,
Density = dexp(
x1,
rate = 1 / 20
),
Kondisi_Label = "S1 - Hujan Sedang"
)
# S2 - Gamma(k = 2, theta = 30)
x2 <- seq(
0,
quantile(
data_simulasi$Curah_Hujan[
data_simulasi$Kondisi == 2
],
0.99
),
length.out = 300
)
teori_s2 <- data.frame(
Curah_Hujan = x2,
Density = dgamma(
x2,
shape = 2,
scale = 30
),
Kondisi_Label = "S2 - Hujan Lebat/Ekstrem"
)
# Menggabungkan data teoritis
data_teoritis <- rbind(
teori_s0,
teori_s1,
teori_s2
)
# VISUALISASI 3
# DISTRIBUSI CURAH HUJAN TEORITIS DAN SIMULASI
ggplot(
data_simulasi,
aes(
x = Curah_Hujan
)
) +
# Histogram simulasi
geom_histogram(
aes(
y = after_stat(density),
fill = Kondisi_Label
),
bins = 30,
alpha = 0.45,
color = "white",
linewidth = 0.2
) +
# Kurva teoritis
geom_line(
data = data_teoritis,
aes(
x = Curah_Hujan,
y = Density,
color = Kondisi_Label
),
linewidth = 1.2
) +
# Rata-rata teoritis
geom_vline(
data = data.frame(
Kondisi_Label = c(
"S0 - Cerah/Berawan",
"S1 - Hujan Sedang",
"S2 - Hujan Lebat/Ekstrem"
),
Mean = c(
2.5,
20,
60
)
),
aes(
xintercept = Mean
),
linetype = "dashed",
color = "#374151",
linewidth = 0.7
) +
# Facet
facet_wrap(
~Kondisi_Label,
scales = "free"
) +
# Warna histogram
scale_fill_manual(
values = warna_state
) +
# Warna kurva
scale_color_manual(
values = warna_state
) +
labs(
title = "Distribusi Curah Hujan Teoritis dan Hasil Simulasi",
subtitle = "Histogram menunjukkan hasil simulasi dan garis menunjukkan fungsi densitas teoritis",
x = "Curah Hujan (mm)",
y = "Kepadatan"
) +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
plot.subtitle = element_text(
hjust = 0.5,
size = 10
),
strip.text = element_text(
face = "bold",
size = 10
),
legend.position = "none",
panel.grid.minor = element_blank()
)
# 10. RATA-RATA CURAH HUJAN KESELURUHAN
mean_simulasi <- mean(
data_simulasi$Curah_Hujan
)
# Nilai teoritis
mean_teoritis <-
(2.5 * 0.3647799) +
(20 * 0.4276729) +
(60 * 0.2075472)
mean_simulasi
## [1] 22.97607
mean_teoritis
## [1] 21.91824
# Galat relatif
error_mean <- abs(
mean_simulasi - mean_teoritis
) / mean_teoritis * 100
error_mean
## [1] 4.826241
# 11. PROPORSI CURAH HUJAN > 50 MM
jumlah_lebih_50 <- sum(
data_simulasi$Curah_Hujan > 50
)
proporsi_lebih_50 <- mean(
data_simulasi$Curah_Hujan > 50
)
persentase_lebih_50 <-
proporsi_lebih_50 * 100
jumlah_lebih_50
## [1] 148
proporsi_lebih_50
## [1] 0.148
persentase_lebih_50
## [1] 14.8
# Nilai teoritis
teori_prob_50 <-
(0 * 0.3647799) +
(0.082084999 * 0.4276729) +
(0.503668 * 0.2075472)
teori_persentase_50 <-
teori_prob_50 * 100
teori_prob_50
## [1] 0.1396404
teori_persentase_50
## [1] 13.96404
# Galat relatif
error_50 <- abs(
proporsi_lebih_50 - teori_prob_50
) / teori_prob_50 * 100
error_50
## [1] 5.98651
# 12. TABEL PERBANDINGAN
tabel_perbandingan <- data.frame(
Ukuran = c(
"Rata-rata Curah Hujan (mm)",
"P(R > 50 mm)"
),
Teoritis = c(
mean_teoritis,
teori_prob_50 * 100
),
Simulasi = c(
mean_simulasi,
persentase_lebih_50
),
Galat_Relatif = c(
error_mean,
error_50
)
)
tabel_perbandingan
## Ukuran Teoritis Simulasi Galat_Relatif
## 1 Rata-rata Curah Hujan (mm) 21.91824 22.97607 4.826241
## 2 P(R > 50 mm) 13.96404 14.80000 5.986510
# 13. VISUALISASI PERBANDINGAN HASIL TEORITIS DAN SIMULASI
data_perbandingan_plot <- data.frame(
Ukuran = rep(
c(
"Rata-rata Curah Hujan",
"Proporsi R > 50 mm"
),
each = 2
),
Metode = rep(
c(
"Teoritis",
"Simulasi"
),
2
),
Nilai = c(
mean_teoritis,
mean_simulasi,
teori_prob_50 * 100,
persentase_lebih_50
)
)
ggplot(
data_perbandingan_plot,
aes(
x = Metode,
y = Nilai,
fill = Metode
)
) +
geom_col(
width = 0.60
) +
geom_text(
aes(
label = round(
Nilai,
2
)
),
vjust = -0.4,
fontface = "bold",
size = 4
) +
facet_wrap(
~Ukuran,
scales = "free_y"
) +
scale_fill_manual(
values = c(
"Teoritis" = "#263F66",
"Simulasi" = "#F4B942"
)
) +
labs(
title = "Perbandingan Hasil Teoritis dan Simulasi",
x = NULL,
y = "Nilai",
fill = "Metode"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(
hjust = 0.5,
face = "bold",
size = 14
),
strip.text = element_text(
face = "bold",
size = 11
),
legend.position = "bottom",
panel.grid.minor = element_blank(),
axis.title = element_text(
face = "bold"
)
)
# SELESAI