Final Project Proses Stokastik CTMC

data$Time.Arrival <- parse_date_time(data$Time.Arrival, orders = "H:M:S")
data$Service.Duration <- hms(data$Service.Duration)

interarrival <- as.numeric(diff(data$Time.Arrival), units = "secs")
k1 <- ceiling(log2(length(interarrival)) + 1)
binwidth1 <- (max(interarrival) - min(interarrival)) / k1
p1 <- ggplot(data.frame(interarrival), aes(x = interarrival)) +
  geom_histogram(binwidth = binwidth1, fill = "blue", color = "black") +
  labs(title = "Distribusi Waktu Antar Kedatangan (Inter-Arrival Times)",
       x = "Waktu Antar Kedatangan (detik)", y = "Frekuensi") +
  theme_minimal()
p1

ks.test(interarrival, "pexp", rate=1/mean(interarrival))

## Warning in ks.test.default(interarrival, "pexp", rate = 1/mean(interarrival)):
## ties should not be present for the one-sample Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  interarrival
## D = 0.080798, p-value = 0.5378
## alternative hypothesis: two-sided

service_time <- as.numeric(data$Service.Duration, units = "secs")
service_time

##   [1]   7 134  54  54  63  42 273  60 198  54 443  45 120  43  61 115 284 139
##  [19]  72 235  55 272 197 234  36 201 110  83  73  34  73 117  11 113 109  54
##  [37]  41  29  10 179 155  11 102  30   8  22 519  59 156  76  89  23 218 238
##  [55]  73  16  96  10  61  49 128  42  36 138 310 181 113  88  55  86  27  86
##  [73]  34  62 152  79  95 230  25 118 155  15  27 144  77  85 214 354 115  97
##  [91] 181  47  44 125 107 244  39  82 105  56

k <- ceiling(log2(length(service_time)) + 1)
binwidth <- (max(service_time) - min(service_time)) / k
p2 <- ggplot(data.frame(service_time), aes(x = service_time)) +
  geom_histogram(binwidth = binwidth, fill = "forestgreen", color = "black") +
  labs(title = "Distribusi Waktu Pelayanan (Service Times)",
       x = "Waktu Pelayanan (detik)", y = "Frekuensi") +
  theme_minimal()
p2

ks.test(service_time, "pexp", rate=1/mean(service_time))

## Warning in ks.test.default(service_time, "pexp", rate = 1/mean(service_time)):
## ties should not be present for the one-sample Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  service_time
## D = 0.11721, p-value = 0.1281
## alternative hypothesis: two-sided

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

lambda <- 1 / mean(interarrival)  # kedatangan per detik
mu <- 1 / mean(service_time)      # pelayanan per detik
lambda <- lambda*3600
mu <- mu*3600
rho <- lambda/mu
cat("λ =", round(lambda, 4), "kedatangan/jam\n")

## λ = 40.2666 kedatangan/jam

cat("μ =", round(mu, 4), "pelayanan/jam\n")

## μ = 32.9188 pelayanan/jam

cat("ρ =", round(rho, 4), "\n")

## ρ = 1.2232

Diagram Rate

library(DiagrammeR)

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

m <- 15

# Buat list node
node_str <- paste0("  ", 0:m, " [label = '", 0:m, "']", collapse = "\n")

# Buat list edge
edge_str <- ""
for (i in 0:(m-1)) {
  edge_str <- paste0(edge_str,
                     sprintf("  %d -> %d [label = 'λ = %.2f']\n", i, i+1, lambda),
                     sprintf("  %d -> %d [label = 'μ = %.2f']\n", i+1, i, mu))
}

# Gabungkan dan tampilkan graf
grViz(paste0("
digraph CTMC {
  graph [layout = dot, rankdir = LR]
  node [shape = circle, style = filled, color = lightgray]

", node_str, "\n\n", edge_str, "
}
"))

m_full <- 15
n_full <- m_full + 1

# Buat matrix rate penuh
R_full <- matrix(0, nrow = n_full, ncol = n_full)
dimnames(R_full) <- list(paste0("[", 0:m_full, ",]"), paste0("[,", 0:m_full, "]"))

# Isi transisi antar state
for (i in 1:m_full) {
  R_full[i, i+1] <- lambda      # arrival
  R_full[i+1, i] <- mu          # service
}
R_full

##          [,0]     [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
## [0,]   0.0000 40.26664  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [1,]  32.9188  0.00000 40.26664  0.00000  0.00000  0.00000  0.00000  0.00000
## [2,]   0.0000 32.91880  0.00000 40.26664  0.00000  0.00000  0.00000  0.00000
## [3,]   0.0000  0.00000 32.91880  0.00000 40.26664  0.00000  0.00000  0.00000
## [4,]   0.0000  0.00000  0.00000 32.91880  0.00000 40.26664  0.00000  0.00000
## [5,]   0.0000  0.00000  0.00000  0.00000 32.91880  0.00000 40.26664  0.00000
## [6,]   0.0000  0.00000  0.00000  0.00000  0.00000 32.91880  0.00000 40.26664
## [7,]   0.0000  0.00000  0.00000  0.00000  0.00000  0.00000 32.91880  0.00000
## [8,]   0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 32.91880
## [9,]   0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [10,]  0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [11,]  0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [12,]  0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [13,]  0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [14,]  0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [15,]  0.0000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
##           [,8]     [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]
## [0,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [1,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [2,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [3,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [4,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [5,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [6,]   0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [7,]  40.26664  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [8,]   0.00000 40.26664  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
## [9,]  32.91880  0.00000 40.26664  0.00000  0.00000  0.00000  0.00000  0.00000
## [10,]  0.00000 32.91880  0.00000 40.26664  0.00000  0.00000  0.00000  0.00000
## [11,]  0.00000  0.00000 32.91880  0.00000 40.26664  0.00000  0.00000  0.00000
## [12,]  0.00000  0.00000  0.00000 32.91880  0.00000 40.26664  0.00000  0.00000
## [13,]  0.00000  0.00000  0.00000  0.00000 32.91880  0.00000 40.26664  0.00000
## [14,]  0.00000  0.00000  0.00000  0.00000  0.00000 32.91880  0.00000 40.26664
## [15,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 32.91880  0.00000

m_generator <- 15
n_generator <- m_generator+1
Q_full <- matrix(0, nrow = n_generator, ncol = n_generator)
dimnames(Q_full) <- list(paste0("[", 0:m_generator, ",]"), paste0("[,", 0:m_generator, "]"))
for (i in 1:m_generator) {
  Q_full[i, i+1] <- lambda      # arrival
  Q_full[i+1, i] <- mu          # service
} 
for (i in 1:n_generator) {
  Q_full[i, i] <- -sum(Q_full[i, ])
}
Q_full

##            [,0]      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [0,]  -40.26664  40.26664   0.00000   0.00000   0.00000   0.00000   0.00000
## [1,]   32.91880 -73.18544  40.26664   0.00000   0.00000   0.00000   0.00000
## [2,]    0.00000  32.91880 -73.18544  40.26664   0.00000   0.00000   0.00000
## [3,]    0.00000   0.00000  32.91880 -73.18544  40.26664   0.00000   0.00000
## [4,]    0.00000   0.00000   0.00000  32.91880 -73.18544  40.26664   0.00000
## [5,]    0.00000   0.00000   0.00000   0.00000  32.91880 -73.18544  40.26664
## [6,]    0.00000   0.00000   0.00000   0.00000   0.00000  32.91880 -73.18544
## [7,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  32.91880
## [8,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [9,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [10,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [11,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [12,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [13,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [14,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [15,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
##            [,7]      [,8]      [,9]     [,10]     [,11]     [,12]     [,13]
## [0,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [1,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [2,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [3,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [4,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [5,]    0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [6,]   40.26664   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
## [7,]  -73.18544  40.26664   0.00000   0.00000   0.00000   0.00000   0.00000
## [8,]   32.91880 -73.18544  40.26664   0.00000   0.00000   0.00000   0.00000
## [9,]    0.00000  32.91880 -73.18544  40.26664   0.00000   0.00000   0.00000
## [10,]   0.00000   0.00000  32.91880 -73.18544  40.26664   0.00000   0.00000
## [11,]   0.00000   0.00000   0.00000  32.91880 -73.18544  40.26664   0.00000
## [12,]   0.00000   0.00000   0.00000   0.00000  32.91880 -73.18544  40.26664
## [13,]   0.00000   0.00000   0.00000   0.00000   0.00000  32.91880 -73.18544
## [14,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  32.91880
## [15,]   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
##           [,14]     [,15]
## [0,]    0.00000   0.00000
## [1,]    0.00000   0.00000
## [2,]    0.00000   0.00000
## [3,]    0.00000   0.00000
## [4,]    0.00000   0.00000
## [5,]    0.00000   0.00000
## [6,]    0.00000   0.00000
## [7,]    0.00000   0.00000
## [8,]    0.00000   0.00000
## [9,]    0.00000   0.00000
## [10,]   0.00000   0.00000
## [11,]   0.00000   0.00000
## [12,]   0.00000   0.00000
## [13,]  40.26664   0.00000
## [14,] -73.18544  40.26664
## [15,]  32.91880 -32.91880

Uniformisasi

library(expm)      # untuk matrix exponentiation

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

## Loading required package: Matrix

## 
## Attaching package: 'expm'

## The following object is masked from 'package:Matrix':
## 
##     expm

library(ggplot2)   # untuk visualisasi

v <- max(abs(diag(Q_full)))

# Transition probability matrix P
P <- diag(n_generator) + Q_full / v
P

##            [,0]      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [0,]  0.4497999 0.5502001 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [1,]  0.4497999 0.0000000 0.5502001 0.0000000 0.0000000 0.0000000 0.0000000
## [2,]  0.0000000 0.4497999 0.0000000 0.5502001 0.0000000 0.0000000 0.0000000
## [3,]  0.0000000 0.0000000 0.4497999 0.0000000 0.5502001 0.0000000 0.0000000
## [4,]  0.0000000 0.0000000 0.0000000 0.4497999 0.0000000 0.5502001 0.0000000
## [5,]  0.0000000 0.0000000 0.0000000 0.0000000 0.4497999 0.0000000 0.5502001
## [6,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.4497999 0.0000000
## [7,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.4497999
## [8,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [9,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [10,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [11,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [12,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [13,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [14,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [15,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
##            [,7]      [,8]      [,9]     [,10]     [,11]     [,12]     [,13]
## [0,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [1,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [3,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [4,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [5,]  0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [6,]  0.5502001 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [7,]  0.0000000 0.5502001 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [8,]  0.4497999 0.0000000 0.5502001 0.0000000 0.0000000 0.0000000 0.0000000
## [9,]  0.0000000 0.4497999 0.0000000 0.5502001 0.0000000 0.0000000 0.0000000
## [10,] 0.0000000 0.0000000 0.4497999 0.0000000 0.5502001 0.0000000 0.0000000
## [11,] 0.0000000 0.0000000 0.0000000 0.4497999 0.0000000 0.5502001 0.0000000
## [12,] 0.0000000 0.0000000 0.0000000 0.0000000 0.4497999 0.0000000 0.5502001
## [13,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.4497999 0.0000000
## [14,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.4497999
## [15,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
##           [,14]     [,15]
## [0,]  0.0000000 0.0000000
## [1,]  0.0000000 0.0000000
## [2,]  0.0000000 0.0000000
## [3,]  0.0000000 0.0000000
## [4,]  0.0000000 0.0000000
## [5,]  0.0000000 0.0000000
## [6,]  0.0000000 0.0000000
## [7,]  0.0000000 0.0000000
## [8,]  0.0000000 0.0000000
## [9,]  0.0000000 0.0000000
## [10,] 0.0000000 0.0000000
## [11,] 0.0000000 0.0000000
## [12,] 0.0000000 0.0000000
## [13,] 0.5502001 0.0000000
## [14,] 0.0000000 0.5502001
## [15,] 0.4497999 0.5502001

# Function untuk menghitung P(t) dengan uniformisasi
uniformized_prob <- function(t, P, v, initial_state, tol = 1e-6, max_k = 1000) {
  n <- nrow(P)
  probs <- numeric(n)
  probs[initial_state + 1] <- 1

  result <- numeric(n)
  cumulative_weight <- 0
  k <- 0

  while (k <= max_k && cumulative_weight < 1 - tol) {
    weight <- dpois(k, lambda = v * t)  # lebih stabil daripada manual
    result <- result + weight * (probs %*% (P %^% k))
    cumulative_weight <- cumulative_weight + weight
    k <- k + 1
  }

  return(as.vector(result))
}
t <- 1  # jam
initial_state <- 0
pt1 <- uniformized_prob(t = t, P = P, v = v, initial_state = initial_state)
df_pt1 <- data.frame(State = 0:m, Probability = pt1)
print(df_pt1)

##    State Probability
## 1      0  0.02805811
## 2      1  0.03326795
## 3      2  0.03841686
## 4      3  0.04335609
## 5      4  0.04795351
## 6      5  0.05211045
## 7      6  0.05578060
## 8      7  0.05899061
## 9      8  0.06186258
## 10     9  0.06463873
## 11    10  0.06770877
## 12    11  0.07164096
## 13    12  0.07721769
## 14    13  0.08547668
## 15    14  0.09775851
## 16    15  0.11576096

sum(df_pt1$Probability)

## [1] 0.9999991

P_t <- matrix(0, nrow = n_generator, ncol = n_generator)
for (i in 0:m) {
  P_t[i + 1, ] <- uniformized_prob(t = 1, P = P, v = v, initial_state = i)
}
P_t

##              [,1]        [,2]        [,3]       [,4]       [,5]       [,6]
##  [1,] 0.028058113 0.033267954 0.038416863 0.04335609 0.04795351 0.05211045
##  [2,] 0.027197234 0.032267452 0.037305877 0.04217534 0.04675448 0.05095393
##  [3,] 0.025675517 0.030498318 0.035340086 0.04008413 0.04462825 0.04889986
##  [4,] 0.023688965 0.028187458 0.032769596 0.03734539 0.04183802 0.04619744
##  [5,] 0.021419778 0.025545781 0.029826834 0.03420344 0.03862824 0.04307806
##  [6,] 0.019029092 0.022760001 0.026717975 0.03087559 0.03521720 0.03974932
##  [7,] 0.016652333 0.019987385 0.023617378 0.02754675 0.03179210 0.03639109
##  [8,] 0.014397047 0.017353256 0.020664922 0.02436664 0.02850642 0.03315313
##  [9,] 0.012342900 0.014950947 0.017965795 0.02144948 0.02547932 0.03015427
## [10,] 0.010543399 0.012843663 0.015592338 0.01887543 0.02279662 0.02748264
## [11,] 0.009028828 0.011067748 0.013587309 0.01669366 0.02051319 0.02519721
## [12,] 0.007809919 0.009636774 0.011968099 0.01492621 0.01865619 0.02333003
## [13,] 0.006881776 0.008545974 0.010731352 0.01357251 0.01722904 0.02188928
## [14,] 0.006227734 0.007776617 0.009857612 0.01261396 0.01621563 0.02086286
## [15,] 0.005822850 0.007300025 0.009315687 0.01201842 0.01558469 0.02022226
## [16,] 0.005636919 0.007081072 0.009066528 0.01174432 0.01529392 0.01992660
##             [,7]       [,8]       [,9]      [,10]      [,11]      [,12]
##  [1,] 0.05578060 0.05899061 0.06186258 0.06463873 0.06770877 0.07164096
##  [2,] 0.05473470 0.05812850 0.06126017 0.06437241 0.06785338 0.07226786
##  [3,] 0.05287338 0.05659011 0.06018033 0.06388821 0.06809956 0.07337319
##  [4,] 0.05041669 0.05455080 0.05873859 0.06322735 0.06840077 0.07481012
##  [5,] 0.04756876 0.05217312 0.05704180 0.06242769 0.06871337 0.07644207
##  [6,] 0.04451397 0.04960521 0.05518904 0.06152674 0.06900162 0.07814900
##  [7,] 0.04141415 0.04697955 0.05327174 0.06056336 0.06924056 0.07983158
##  [8,] 0.03840675 0.04441161 0.05137317 0.05957794 0.06941707 0.08141285
##  [9,] 0.03560362 0.04199862 0.04956707 0.05861126 0.06952904 0.08283815
## [10,] 0.03309069 0.03981831 0.04791590 0.05770230 0.06958327 0.08407322
## [11,] 0.03092823 0.03792820 0.04646902 0.05688575 0.06959246 0.08510144
## [12,] 0.02915197 0.03636533 0.04526123 0.05618948 0.06957217 0.08592029
## [13,] 0.02777497 0.03514687 0.04431200 0.05563264 0.06953781 0.08653788
## [14,] 0.02679024 0.03427153 0.04362571 0.05522454 0.06950246 0.08696958
## [15,] 0.02617394 0.03372186 0.04319275 0.05496457 0.06947537 0.08723518
## [16,] 0.02588901 0.03346722 0.04299162 0.05484310 0.06946144 0.08735668
##            [,13]      [,14]      [,15]     [,16]
##  [1,] 0.07721769 0.08547668 0.09775851 0.1157610
##  [2,] 0.07839286 0.08725834 0.10019406 0.1188825
##  [3,] 0.08047630 0.09042461 0.10452765 0.1244396
##  [4,] 0.08320939 0.09459431 0.11024585 0.1317784
##  [5,] 0.08635188 0.09941366 0.11687233 0.1402923
##  [6,] 0.08968939 0.10456462 0.12397724 0.1494331
##  [7,] 0.09303829 0.10977062 0.13118372 0.1587185
##  [8,] 0.09624839 0.11479995 0.13817233 0.1677376
##  [9,] 0.09920345 0.11946737 0.14468351 0.1761543
## [10,] 0.10182015 0.12363388 0.15051842 0.1837089
## [11,] 0.10404560 0.12720501 0.15553791 0.1902175
## [12,] 0.10585408 0.13012780 0.15965985 0.1955697
## [13,] 0.10724302 0.13238658 0.16285459 0.1997228
## [14,] 0.10822874 0.13399784 0.16513886 0.2026952
## [15,] 0.10884220 0.13500440 0.16656833 0.2045566
## [16,] 0.10912483 0.13546924 0.16722917 0.2054174

ggplot(df_pt1, aes(x = factor(State), y = Probability)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Probabilitas Transisi dari State 0 ke State 0–15 pada t = 1 Jam",
       x = "State Tujuan", y = "Probabilitas") +
  theme_minimal()

occupancy_time_matrix <- function(t, P, v, tol = 1e-6, max_k = 1000) {
  n <- nrow(P)
  O <- matrix(0, n, n)  # hasil occupancy time
  cumulative_weight <- 0
  k <- 0
  Pm_sum <- diag(n)     # sum_{m=0}^{k-1} P^m

  while (k <= max_k && 1 - cumulative_weight > tol) {
    weight <- dpois(k, lambda = v * t)
    if (k > 0) {
      Pm_sum <- Pm_sum + (P %^% k)  # update akumulasi P^m untuk m=0..k
    }
    O <- O + weight * Pm_sum
    cumulative_weight <- cumulative_weight + weight
    k <- k + 1
  }

  O_matrix <- (1 / v) * O
  dimnames(O_matrix) <- list(paste0("from_", 0:(n - 1)), paste0("to_", 0:(n - 1)))
  message(sprintf("Selesai pada k = %d dengan bobot kumulatif = %.10f", k, cumulative_weight))
  return(O_matrix)
}

O <- occupancy_time_matrix(t = 1, P = P, v = v)

## Selesai pada k = 118 dengan bobot kumulatif = 0.9999991014

print(round(O, 5))

##            to_0    to_1    to_2    to_3    to_4    to_5    to_6    to_7    to_8
## from_0  0.10581 0.09988 0.09363 0.08714 0.08049 0.07379 0.06716 0.06073 0.05464
## from_1  0.08166 0.10070 0.09457 0.08820 0.08166 0.07507 0.06853 0.06218 0.05616
## from_2  0.06258 0.07731 0.09625 0.09010 0.08377 0.07736 0.07100 0.06480 0.05893
## from_3  0.04761 0.05895 0.07366 0.09263 0.08658 0.08044 0.07431 0.06834 0.06267
## from_4  0.03595 0.04462 0.05598 0.07078 0.08991 0.08409 0.07826 0.07257 0.06718
## from_5  0.02694 0.03353 0.04227 0.05376 0.06874 0.08813 0.08267 0.07732 0.07229
## from_6  0.02005 0.02502 0.03171 0.04060 0.05231 0.06758 0.08736 0.08243 0.07782
## from_7  0.01482 0.01856 0.02366 0.03052 0.03965 0.05168 0.06739 0.08777 0.08366
## from_8  0.01090 0.01371 0.01759 0.02289 0.03001 0.03950 0.05201 0.06840 0.08971
## from_9  0.00800 0.01011 0.01307 0.01717 0.02276 0.03028 0.04032 0.05360 0.07105
## from_10 0.00589 0.00748 0.00976 0.01297 0.01739 0.02343 0.03158 0.04249 0.05698
## from_11 0.00439 0.00561 0.00740 0.00994 0.01352 0.01845 0.02520 0.03435 0.04664
## from_12 0.00336 0.00432 0.00576 0.00785 0.01081 0.01497 0.02072 0.02860 0.03931
## from_13 0.00269 0.00348 0.00469 0.00647 0.00903 0.01266 0.01774 0.02477 0.03441
## from_14 0.00230 0.00299 0.00406 0.00566 0.00798 0.01130 0.01598 0.02249 0.03150
## from_15 0.00212 0.00277 0.00378 0.00530 0.00751 0.01069 0.01519 0.02147 0.03019
##            to_9   to_10   to_11   to_12   to_13   to_14   to_15
## from_0  0.04905 0.04419 0.04029 0.03772 0.03693 0.03859 0.04364
## from_1  0.05066 0.04587 0.04208 0.03965 0.03908 0.04105 0.04656
## from_2  0.05356 0.04893 0.04534 0.04319 0.04302 0.04559 0.05193
## from_3  0.05752 0.05313 0.04984 0.04810 0.04851 0.05192 0.05946
## from_4  0.06232 0.05826 0.05538 0.05419 0.05537 0.05986 0.06892
## from_5  0.06779 0.06416 0.06182 0.06133 0.06347 0.06927 0.08017
## from_6  0.07379 0.07070 0.06902 0.06940 0.07269 0.08007 0.09310
## from_7  0.08020 0.07776 0.07690 0.07831 0.08297 0.09217 0.10763
## from_8  0.08691 0.08526 0.08536 0.08799 0.09423 0.10550 0.12369
## from_9  0.09385 0.09312 0.09434 0.09838 0.10641 0.12000 0.14120
## from_10 0.07613 0.10127 0.10376 0.10939 0.11944 0.13560 0.16009
## from_11 0.06305 0.08483 0.11358 0.12098 0.13325 0.15221 0.18025
## from_12 0.05375 0.07311 0.09891 0.13309 0.14777 0.16976 0.20159
## from_13 0.04753 0.06526 0.08906 0.12081 0.16293 0.18814 0.22399
## from_14 0.04382 0.06057 0.08317 0.11345 0.15381 0.20726 0.24732
## from_15 0.04215 0.05846 0.08051 0.11014 0.14970 0.20219 0.27147

First Passages Times

first_passage_time <- function(Q, target_state) {
  n <- nrow(Q)
  idx <- setdiff(1:n, target_state + 1)  # +1 karena R 1-based
  Q_sub <- Q[idx, idx]
  b <- rep(-1, length(idx))
  f_sub <- solve(Q_sub, b)

  f <- rep(NA, n)
  f[idx] <- f_sub
  f[target_state + 1] <- 0
  names(f) <- paste0("f_", 0:(n - 1), "_", target_state)
  return(f)
}

Pertama kali mencapai state 5 dari state 0

f_0_to_5 <- first_passage_time(Q_full, target_state = 5)
f_0_to_5

##      f_0_5      f_1_5      f_2_5      f_3_5      f_4_5      f_5_5      f_6_5 
## 0.29340772 0.26857327 0.22343613 0.16170115 0.08639708 0.00000000 0.88450204 
##      f_7_5      f_8_5      f_9_5     f_10_5     f_11_5     f_12_5     f_13_5 
## 1.58276612 2.12877685 2.55031736 2.87010090 3.10669605 3.27528297 3.38827177 
##     f_14_5     f_15_5 
## 3.45580798 3.48618576

distribution limit

library(expm)  # untuk matrix exponential

steady_state_distribution <- function(Q) {
  n <- nrow(Q)
  A <- t(Q)           # transpos Q untuk sistem linier
  A[n, ] <- 1         # ganti baris terakhir dengan syarat sum π = 1
  b <- c(rep(0, n - 1), 1)
  pi <- solve(A, b)   # menyelesaikan π Q = 0 dengan constraint
  return(pi)
}

pi <- steady_state_distribution(Q_full)
round(pi, 5)

##    [0,]    [1,]    [2,]    [3,]    [4,]    [5,]    [6,]    [7,]    [8,]    [9,] 
## 0.00925 0.01132 0.01385 0.01694 0.02072 0.02534 0.03100 0.03792 0.04638 0.05673 
##   [10,]   [11,]   [12,]   [13,]   [14,]   [15,] 
## 0.06940 0.08489 0.10384 0.12701 0.15537 0.19004

sum(pi)

## [1] 1

Expected Cost

c <- c(0,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30)
c_mat <- matrix(c, nrow = 1)
Mt <- c%*%O
Mt

##          to_0     to_1    to_2     to_3     to_4     to_5     to_6     to_7
## [1,] 9.878231 12.27507 14.5268 16.64502 18.64912 20.56882 22.44757 24.34735
##          to_8     to_9    to_10    to_11    to_12 to_13    to_14   to_15
## [1,] 26.35501 28.59102 31.22113 34.47212 38.65248 44.18 51.61766 61.7207

Expected Cost Infinite Time

expected_cost <- sum(pi * c)

expected_cost

## [1] 29.72238

Final Project Proses Stokastik CTMC

Rochmad Budi Sampurna

2025-07-05