Sintaks berikut merupakan implementasi simulasi proses stokastik cuaca menggunakan model Markov dan distribusi curah hujan.
# persiapan
library(readxl)## Warning: package 'readxl' was built under R version 4.5.3set.seed(123)
N <- 1000
# membaca dataset darwin
data_darwin <- read_excel(
"C:/Users/Hype G12/Downloads/weatherAUS_Darwin_clean.xlsx"
)
data_darwin$Date <- as.Date(data_darwin$Date)
str(data_darwin)## tibble [3,062 × 23] (S3: tbl_df/tbl/data.frame)
## $ Date : Date[1:3062], format: "2008-07-01" "2008-07-02" ...
## $ Location : chr [1:3062] "Darwin" "Darwin" "Darwin" "Darwin" ...
## $ MinTemp : num [1:3062] 20 19.4 18.2 17.3 15.5 16.2 17 19.6 17.3 17.1 ...
## $ MaxTemp : num [1:3062] 33.1 32.4 31.8 30.7 30.8 31.9 32.7 30.8 29.2 30.3 ...
## $ Rainfall : num [1:3062] 0 0 0 0 0 0 0 0 0 0 ...
## $ Evaporation : num [1:3062] 4.4 6 8 7 7 7.2 5.2 9.2 9.6 9.2 ...
## $ Sunshine : num [1:3062] 11 10.4 11 10.4 10.8 10.7 7.8 10.6 10.6 11.1 ...
## $ WindGustDir : chr [1:3062] "E" "ENE" "E" "E" ...
## $ WindGustSpeed: num [1:3062] 41 50 46 44 46 41 48 54 50 48 ...
## $ WindDir9am : chr [1:3062] "ENE" "SE" "ESE" "SE" ...
## $ WindDir3pm : chr [1:3062] "SSE" "E" "ENE" "E" ...
## $ WindSpeed9am : num [1:3062] 13 15 22 22 20 11 11 20 31 17 ...
## $ WindSpeed3pm : num [1:3062] 17 28 19 13 19 13 22 33 24 13 ...
## $ Humidity9am : num [1:3062] 81 81 38 55 37 62 54 36 26 33 ...
## $ Humidity3pm : num [1:3062] 32 17 24 16 16 18 18 12 15 14 ...
## $ Pressure9am : num [1:3062] 1016 1017 1017 1017 1016 ...
## $ Pressure3pm : num [1:3062] 1012 1012 1013 1014 1013 ...
## $ Cloud9am : num [1:3062] 1 1 0 2 1 0 5 1 1 1 ...
## $ Cloud3pm : num [1:3062] 2 1 1 6 1 0 2 0 0 1 ...
## $ Temp9am : num [1:3062] 25.4 24.3 24.3 21.3 22.2 22.8 23.3 22.3 20.3 19.5 ...
## $ Temp3pm : num [1:3062] 32.3 31.9 31.2 29.8 29.6 30.5 31.5 30.2 28.6 28.8 ...
## $ RainToday : chr [1:3062] "No" "No" "No" "No" ...
## $ RainTomorrow : chr [1:3062] "No" "No" "No" "No" ...head(data_darwin)## # A tibble: 6 × 23
## Date Location MinTemp MaxTemp Rainfall Evaporation Sunshine WindGustDir
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 2008-07-01 Darwin 20 33.1 0 4.4 11 E
## 2 2008-07-02 Darwin 19.4 32.4 0 6 10.4 ENE
## 3 2008-07-03 Darwin 18.2 31.8 0 8 11 E
## 4 2008-07-04 Darwin 17.3 30.7 0 7 10.4 E
## 5 2008-07-05 Darwin 15.5 30.8 0 7 10.8 ESE
## 6 2008-07-06 Darwin 16.2 31.9 0 7.2 10.7 E
## # ℹ 15 more variables: WindGustSpeed <dbl>, WindDir9am <chr>, WindDir3pm <chr>,
## # WindSpeed9am <dbl>, WindSpeed3pm <dbl>, Humidity9am <dbl>,
## # Humidity3pm <dbl>, Pressure9am <dbl>, Pressure3pm <dbl>, Cloud9am <dbl>,
## # Cloud3pm <dbl>, Temp9am <dbl>, Temp3pm <dbl>, RainToday <chr>,
## # RainTomorrow <chr>nrow(data_darwin)## [1] 3062sum(is.na(data_darwin))## [1] 0# menentukan parameter model
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("State 0", "State 1", "State 2")
colnames(P) <- c("State 0", "State 1", "State 2")
print(P)## State 0 State 1 State 2
## State 0 0.65 0.30 0.05
## State 1 0.25 0.55 0.20
## State 2 0.10 0.40 0.50# langkah 1 - menentukan nilai awal S0
S0 <- 0
cat("Kondisi awal S0 =", S0, "\n")## Kondisi awal S0 = 0# langkah 2 - membangkitkan state S(t+1)
S <- numeric(N)
S[1] <- S0
for (t in 1:(N - 1)) {
S[t + 1] <- sample(
0:2,
size = 1,
prob = P[S[t] + 1, ]
)
}
# langkah 3- membangkitkan curah hujan R(t)
R <- numeric(N)
for (t in 1:N) {
if (S[t] == 0) {
R[t] <- runif(
1,
min = 0,
max = 5
)
} else if (S[t] == 1) {
R[t] <- rexp(
1,
rate = 1 / 20
)
} else {
R[t] <- rgamma(
1,
shape = 2,
scale = 30
)
}
}
simulasi <- data.frame(
Hari = 1:N,
State = S,
Rainfall = R
)
head(simulasi)## Hari State Rainfall
## 1 1 0 0.5441204
## 2 2 0 1.3681137
## 3 3 1 3.7546773
## 4 4 1 33.3554574
## 5 5 2 91.9143284
## 6 6 0 2.3894341# langkah 4 - rata rata curah hujan tiap rigme
mean_state_0 <- mean(R[S == 0])
mean_state_1 <- mean(R[S == 1])
mean_state_2 <- mean(R[S == 2])
rata_rata_regime <- data.frame(
State = 0:2,
Rata_rata_Simulasi = c(
mean_state_0,
mean_state_1,
mean_state_2
)
)
print(rata_rata_regime)## State Rata_rata_Simulasi
## 1 0 2.618053
## 2 1 20.586376
## 3 2 58.026557# langkah 5 - rata-rata curah hujan keseluruhan
mean_total <- mean(R)
cat(
"Rata-rata curah hujan keseluruhan =",
mean_total,
"mm\n"
)## Rata-rata curah hujan keseluruhan = 22.97607 mm# langkah 6 - proporsi hari dengan R > 50 mm
jumlah_hujan_lebih_50 <- sum(R > 50)
proporsi_lebih_50 <- mean(R > 50)
cat(
"Jumlah hari dengan R > 50 mm =",
jumlah_hujan_lebih_50,
"hari\n"
)## Jumlah hari dengan R > 50 mm = 148 haricat(
"Proporsi hari dengan R > 50 mm =",
proporsi_lebih_50,
"\n"
)## Proporsi hari dengan R > 50 mm = 0.148cat(
"Persentase hari dengan R > 50 mm =",
proporsi_lebih_50 * 100,
"%\n"
)## Persentase hari dengan R > 50 mm = 14.8 %# langkah 7 - Perbandingan Hasil Simulasi dan Teori
E_R_state0 <- 2.5
E_R_state1 <- 20
E_R_state2 <- 60
P_50_state0 <- 0
P_50_state1 <- exp(-50 / 20)
P_50_state2 <- exp(-50 / 30) * (1 + 50 / 30)
A <- t(P) - diag(3)
A[3, ] <- c(1, 1, 1)
b <- c(0, 0, 1)
pi <- solve(A, b)
names(pi) <- c(
"State 0",
"State 1",
"State 2"
)
print(pi)## State 0 State 1 State 2
## 0.3647799 0.4276730 0.2075472E_R_theory <-
pi[1] * E_R_state0 +
pi[2] * E_R_state1 +
pi[3] * E_R_state2
cat(
"Rata-rata curah hujan teoritis =",
E_R_theory,
"mm\n"
)## Rata-rata curah hujan teoritis = 21.91824 mmP_50_theory <-
pi[1] * P_50_state0 +
pi[2] * P_50_state1 +
pi[3] * P_50_state2
cat(
"P(R > 50) teoritis =",
P_50_theory,
"\n"
)## P(R > 50) teoritis = 0.1396405cat(
"Persentase teoritis R > 50 =",
P_50_theory * 100,
"%\n"
)## Persentase teoritis R > 50 = 13.96405 %perbandingan_state <- data.frame(
State = 0:2,
Teoritis = c(
E_R_state0,
E_R_state1,
E_R_state2
),
Simulasi = c(
mean_state_0,
mean_state_1,
mean_state_2
)
)
perbandingan_state$Galat_Relatif_Persen <-
abs(
perbandingan_state$Simulasi -
perbandingan_state$Teoritis
) /
perbandingan_state$Teoritis * 100
print(perbandingan_state)## State Teoritis Simulasi Galat_Relatif_Persen
## 1 0 2.5 2.618053 4.722106
## 2 1 20.0 20.586376 2.931878
## 3 2 60.0 58.026557 3.289072galat_mean <-
abs(mean_total - E_R_theory) /
E_R_theory * 100
cat(
"Galat relatif rata-rata keseluruhan =",
galat_mean,
"%\n"
)## Galat relatif rata-rata keseluruhan = 4.826245 %galat_P50 <-
abs(proporsi_lebih_50 - P_50_theory) /
P_50_theory * 100
cat(
"Galat relatif P(R > 50) =",
galat_P50,
"%\n"
)## Galat relatif P(R > 50) = 5.986475 %hasil_akhir <- data.frame(
Ukuran = c(
"Rata-rata State 0",
"Rata-rata State 1",
"Rata-rata State 2",
"Rata-rata keseluruhan",
"P(R > 50)"
),
Teoritis = c(
E_R_state0,
E_R_state1,
E_R_state2,
E_R_theory,
P_50_theory
),
Simulasi = c(
mean_state_0,
mean_state_1,
mean_state_2,
mean_total,
proporsi_lebih_50
)
)
hasil_akhir$Galat_Relatif_Persen <-
abs(
hasil_akhir$Simulasi -
hasil_akhir$Teoritis
) /
abs(hasil_akhir$Teoritis) * 100
print(hasil_akhir)## Ukuran Teoritis Simulasi Galat_Relatif_Persen
## 1 Rata-rata State 0 2.5000000 2.618053 4.722106
## 2 Rata-rata State 1 20.0000000 20.586376 2.931878
## 3 Rata-rata State 2 60.0000000 58.026557 3.289072
## 4 Rata-rata keseluruhan 21.9182390 22.976067 4.826245
## 5 P(R > 50) 0.1396405 0.148000 5.986475