PROJECT 1 PROSES STOKASTIK EKSPEKTASI BERSYARAT DISKRIT DAN KONTINU

Simulasi Proses Stokastik Curah Hujan

Jumlah hari: 1000

Kondisi cuaca:
0 = Cerah/Berawan
1 = Hujan Sedang
2 = Hujan Lebat/Ekstrem

1. LOAD PACKAGE

# PROJECT 1 PROSES STOKASTIK EKSPEKTASI BERSYARAT DISKRIT DAN KONTINU
# Simulasi Proses Stokastik Curah Hujan
# Jumlah hari        : 1000
# Kondisi cuaca     : 0 = Cerah/Berawan, 1 = Hujan Sedang, 2 = Hujan Lebat/Ekstrem

# 1. LOAD PACKAGE
library(ggplot2)
library(igraph)

2. MENENTUKAN PARAMETER MODEL

# 2. MENENTUKAN PARAMETER MODEL
# 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")
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

# 3. SIMULASI KONDISI CUACA
set.seed(123)
# Membuat vektor untuk menyimpan kondisi cuaca
S = numeric(n)
# Menentukan kondisi awal
S[1] = S0
# Simulasi kondisi cuaca berdasarkan matriks transisi
for (t in 1:(n - 1)) {
  S[t + 1] <- sample(
    x = 0:2,
    size = 1,
    prob = P[S[t] + 1, ]
  )
}

4. SIMULASI CURAH HUJAN

# 4. SIMULASI CURAH HUJAN
# Membuat vektor untuk menyimpan curah hujan
R = numeric(n)
# Menghasilkan curah hujan berdasarkan kondisi cuaca
for (t in 1:n) {
  if (S[t] == 0) {
    # Kondisi 0: Uniform(0,5)
    R[t] = runif(1, min = 0, max = 5)
  } else if (S[t] == 1) {
    # Kondisi 1: Exponential(lambda = 1/20)
    R[t] = rexp(1, rate = 1 / 20)
  } else {
    # Kondisi 2: Gamma(k = 2, theta = 30)
    R[t] = rgamma(1, shape = 2, scale = 30)
  }
}

5. MEMBUAT DATA FRAME HASIL SIMULASI

# 5. MEMBUAT DATA FRAME HASIL SIMULASI
data_simulasi <- data.frame(Hari = 1:n, Kondisi = S, Curah_Hujan = R)
# Memberikan 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")
)
head(data_simulasi)

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

6. HASIL SIMULASI KONDISI CUACA

# 6. HASIL SIMULASI KONDISI CUACA
jumlah_kondisi = table(data_simulasi$Kondisi)
jumlah_kondisi

## 
##   0   1   2 
## 340 433 227

proporsi_kondisi = prop.table(jumlah_kondisi)
proporsi_kondisi

## 
##     0     1     2 
## 0.340 0.433 0.227

persentase_kondisi = proporsi_kondisi * 100
persentase_kondisi

## 
##    0    1    2 
## 34.0 43.3 22.7

7. RATA-RATA CURAH HUJAN HASIL SIMULASI

# 7. RATA-RATA CURAH HUJAN HASIL SIMULASI
mean_per_kondisi = tapply(data_simulasi$Curah_Hujan, data_simulasi$Kondisi, mean)
mean_per_kondisi

##         0         1         2 
##  2.618053 20.586376 58.026557

mean_simulasi = mean(data_simulasi$Curah_Hujan)
mean_simulasi

## [1] 22.97607

8. PELUANG CURAH HUJAN LEBIH DARI 50 MM HASIL SIMULASI

# 8. PELUANG CURAH HUJAN LEBIH DARI 50 MM HASIL SIMULASI
jumlah_lebih_50 = sum(data_simulasi$Curah_Hujan > 50)
jumlah_lebih_50

## [1] 148

proporsi_lebih_50 = mean(data_simulasi$Curah_Hujan > 50)
proporsi_lebih_50

## [1] 0.148

persentase_lebih_50 = proporsi_lebih_50 * 100
persentase_lebih_50

## [1] 14.8

9. NILAI TEORITIS EKSPEKTASI CURAH HUJAN

# 9. NILAI TEORITIS EKSPEKTASI CURAH HUJAN
# Ekspektasi curah hujan bersyarat
# S0 : Uniform(0,5)       = E(R|S0) = 2.5
# S1 : Exponential(1/20)  = E(R|S1) = 20
# S2 : Gamma(2,30)        = E(R|S2) = 60
teori_mean_kondisi <- c(2.5, 20, 60)
names(teori_mean_kondisi) <- c("S0", "S1", "S2")
teori_mean_kondisi

##   S0   S1   S2 
##  2.5 20.0 60.0

10. DISTRIBUSI STASIONER

# 10. DISTRIBUSI STASIONER
# Distribusi stasioner: pi = (58/159, 68/159, 33/159)
pi = c(58 / 159, 68 / 159, 33 / 159)
names(pi) = c("S0", "S1", "S2")
pi

##        S0        S1        S2 
## 0.3647799 0.4276730 0.2075472

11. NILAI TEORITIS RATA-RATA CURAH HUJAN JANGKA PANJANG

# 11. NILAI TEORITIS RATA-RATA CURAH HUJAN JANGKA PANJANG
teori_mean_total = sum(teori_mean_kondisi * pi)
teori_mean_total

## [1] 21.91824

12. NILAI TEORITIS P(R > 50 | S = i)

# 12. NILAI TEORITIS P(R > 50 | S = i)
# Kondisi S0: Uniform(0,5), sehingga P(R > 50 | S0) = 0
tail_S0 = 0

# Kondisi S1: Exponential(lambda = 1/20), P(R > 50 | S1) = exp(-50/20)
tail_S1 = exp(-50 / 20)

# Kondisi S2: Gamma(k=2, theta=30), P(R > 50 | S2) = exp(-50/30) * (1 + 50/30)
tail_S2 = exp(-50 / 30) * (1 + 50 / 30)

# Menggabungkan peluang bersyarat
teori_tail_kondisi = c(tail_S0, tail_S1, tail_S2)
names(teori_tail_kondisi) = c("S0", "S1", "S2")
teori_tail_kondisi

##        S0        S1        S2 
## 0.0000000 0.0820850 0.5036683

13. NILAI TEORITIS P(R > 50) JANGKA PANJANG

# 13. NILAI TEORITIS P(R > 50) JANGKA PANJANG
teori_tail_total = sum(teori_tail_kondisi * pi)
teori_tail_total

## [1] 0.1396405

teori_persentase_lebih_50 <- teori_tail_total * 100
teori_persentase_lebih_50

## [1] 13.96405

14. PERHITUNGAN ERROR RELATIF

# 14. PERHITUNGAN ERROR RELATIF
# Error relatif rata-rata curah hujan
error_mean = abs(mean_simulasi - teori_mean_total) / teori_mean_total * 100
error_mean

## [1] 4.826245

# Error relatif peluang R > 50
error_tail = abs(proporsi_lebih_50 - teori_tail_total) / teori_tail_total * 100
error_tail

## [1] 5.986475

15. TABEL PERBANDINGAN TEORITIS DAN SIMULASI

# 15. TABEL PERBANDINGAN TEORITIS DAN SIMULASI
data_perbandingan <- data.frame(
  Ukuran = c("Rata-rata Curah Hujan", "P(R > 50 mm)"),
  Teoritis = c(teori_mean_total, teori_persentase_lebih_50),
  Simulasi = c(mean_simulasi, persentase_lebih_50),
  Error_Relatif = c(error_mean, error_tail)
)
data_perbandingan

##                  Ukuran Teoritis Simulasi Error_Relatif
## 1 Rata-rata Curah Hujan 21.91824 22.97607      4.826245
## 2          P(R > 50 mm) 13.96405 14.80000      5.986475

16. PERBANDINGAN DISTRIBUSI KONDISI

# 16. PERBANDINGAN DISTRIBUSI KONDISI
data_kondisi = data.frame(
  Kondisi = c("S0", "S1", "S2"),
  Teoritis = pi * 100,
  Simulasi = as.numeric(persentase_kondisi)
)
data_kondisi

##    Kondisi Teoritis Simulasi
## S0      S0 36.47799     34.0
## S1      S1 42.76730     43.3
## S2      S2 20.75472     22.7

17. VISUALISASI 1 DIAGRAM TRANSISI KONDISI CUACA

# 17. VISUALISASI 1 DIAGRAM TRANSISI KONDISI CUACA
edges = expand.grid(from = c("S0", "S1", "S2"), to = c("S0", "S1", "S2"))
edges$prob = as.vector(P)
edges$label = sprintf("%.2f", edges$prob)
nama_kondisi = c("S0", "S1", "S2")
g = graph_from_data_frame(edges, directed = TRUE, vertices = data.frame(name = nama_kondisi, label = nama_kondisi))
layout_matrix = matrix(c(0, 1, 1, 0, 2, 1), ncol = 2, byrow = TRUE)
plot(g, layout = layout_matrix, vertex.size = 35, vertex.label.cex = 1.2, edge.arrow.size = 0.4,
     edge.label = edges$label, edge.label.cex = 0.8, main = "Diagram Transisi Kondisi Cuaca")

18. VISUALISASI 2 DISTRIBUSI TEORITIS DAN HASIL SIMULASI

# 18. VISUALISASI 2 DISTRIBUSI TEORITIS DAN HASIL SIMULASI
# Kondisi S0: Uniform(0,5)
x0 = seq(0, 5, length.out = 300)
x1 = seq(0, quantile(data_simulasi$Curah_Hujan[data_simulasi$Kondisi == 1], 0.99), length.out = 300)
x2 = seq(0, quantile(data_simulasi$Curah_Hujan[data_simulasi$Kondisi == 2], 0.99), length.out = 300)

teori_s0 = data.frame(Curah_Hujan = x0, Density = dunif(x0, min = 0, max = 5), Kondisi_Label = "S0 - Cerah/Berawan")
teori_s1 = data.frame(Curah_Hujan = x1, Density = dexp(x1, rate = 1 / 20), Kondisi_Label = "S1 - Hujan Sedang")
teori_s2 = data.frame(Curah_Hujan = x2, Density = dgamma(x2, shape = 2, scale = 30), Kondisi_Label = "S2 - Hujan Lebat/Ekstrem")
data_teoritis <- rbind(teori_s0, teori_s1, teori_s2)

ggplot(data_simulasi, aes(x = Curah_Hujan)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, alpha = 0.6) +
  geom_line(data = data_teoritis, aes(x = Curah_Hujan, y = Density), linewidth = 1) +
  facet_wrap(~Kondisi_Label, scales = "free") +
  labs(title = "Distribusi Curah Hujan Teoritis dan Hasil Simulasi", x = "Curah Hujan (mm)", y = "Kepadatan") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), strip.text = element_text(face = "bold"))

19. VISUALISASI 3 PERBANDINGAN NILAI TEORITIS DAN SIMULASI

# 19. VISUALISASI 3 PERBANDINGAN NILAI TEORITIS DAN SIMULASI
data_perbandingan_grafik <- data.frame(
  Ukuran = rep(c("Rata-rata Curah Hujan", "Proporsi R > 50 mm"), each = 2),
  Metode = rep(c("Teoritis", "Simulasi"), 2),
  Nilai = c(teori_mean_total, mean_simulasi, teori_persentase_lebih_50, persentase_lebih_50)
)

ggplot(data_perbandingan_grafik, aes(x = Metode, y = Nilai)) +
  geom_col(width = 0.6, alpha = 0.8) +
  geom_text(aes(label = round(Nilai, 2)), vjust = -0.4, size = 4) +
  facet_wrap(~Ukuran, scales = "free_y") +
  labs(title = "Perbandingan Hasil Teoritis dan Simulasi", x = NULL, y = "Nilai") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"), strip.text = element_text(face = "bold"))

20. RINGKASAN HASIL

# 20. RINGKASAN HASIL
cat("\nJumlah hari simulasi :", n)

## 
## Jumlah hari simulasi : 1000

cat("\n\nDistribusi kondisi hasil simulasi:\n")

## 
## 
## Distribusi kondisi hasil simulasi:

print(jumlah_kondisi)

## 
##   0   1   2 
## 340 433 227

cat("\nProporsi kondisi hasil simulasi:\n")

## 
## Proporsi kondisi hasil simulasi:

print(proporsi_kondisi)

## 
##     0     1     2 
## 0.340 0.433 0.227

cat("\nPersentase kondisi hasil simulasi:\n")

## 
## Persentase kondisi hasil simulasi:

print(persentase_kondisi)

## 
##    0    1    2 
## 34.0 43.3 22.7

cat("\nRata-rata curah hujan berdasarkan kondisi:\n")

## 
## Rata-rata curah hujan berdasarkan kondisi:

print(mean_per_kondisi)

##         0         1         2 
##  2.618053 20.586376 58.026557

cat("\nRata-rata curah hujan simulasi :", mean_simulasi, "mm")

## 
## Rata-rata curah hujan simulasi : 22.97607 mm

cat("\nRata-rata curah hujan teoritis :", teori_mean_total, "mm")

## 
## Rata-rata curah hujan teoritis : 21.91824 mm

cat("\nError relatif rata-rata :", error_mean, "%")

## 
## Error relatif rata-rata : 4.826245 %

cat("\n\nJumlah R > 50 mm :", jumlah_lebih_50)

## 
## 
## Jumlah R > 50 mm : 148

cat("\nProporsi R > 50 mm simulasi :", proporsi_lebih_50)

## 
## Proporsi R > 50 mm simulasi : 0.148

cat("\nPersentase R > 50 mm simulasi :", persentase_lebih_50, "%")

## 
## Persentase R > 50 mm simulasi : 14.8 %

cat("\nP(R > 50 mm) teoritis :", teori_tail_total)

## 
## P(R > 50 mm) teoritis : 0.1396405

cat("\nPersentase teoritis :", teori_persentase_lebih_50, "%")

## 
## Persentase teoritis : 13.96405 %

cat("\nError relatif peluang :", error_tail, "%")

## 
## Error relatif peluang : 5.986475 %

Simulasi Komputasi Project 1 Proses Stokastik

Restika (H051241055)

2026-09-27

PROJECT 1 PROSES STOKASTIK EKSPEKTASI BERSYARAT DISKRIT DAN KONTINU

Simulasi Proses Stokastik Curah Hujan

1. LOAD PACKAGE

2. MENENTUKAN PARAMETER MODEL

3. SIMULASI KONDISI CUACA

4. SIMULASI CURAH HUJAN

5. MEMBUAT DATA FRAME HASIL SIMULASI

6. HASIL SIMULASI KONDISI CUACA

7. RATA-RATA CURAH HUJAN HASIL SIMULASI

8. PELUANG CURAH HUJAN LEBIH DARI 50 MM HASIL SIMULASI

9. NILAI TEORITIS EKSPEKTASI CURAH HUJAN

10. DISTRIBUSI STASIONER

11. NILAI TEORITIS RATA-RATA CURAH HUJAN JANGKA PANJANG

12. NILAI TEORITIS P(R > 50 | S = i)

13. NILAI TEORITIS P(R > 50) JANGKA PANJANG

14. PERHITUNGAN ERROR RELATIF

15. TABEL PERBANDINGAN TEORITIS DAN SIMULASI

16. PERBANDINGAN DISTRIBUSI KONDISI

17. VISUALISASI 1 DIAGRAM TRANSISI KONDISI CUACA

18. VISUALISASI 2 DISTRIBUSI TEORITIS DAN HASIL SIMULASI

19. VISUALISASI 3 PERBANDINGAN NILAI TEORITIS DAN SIMULASI

20. RINGKASAN HASIL