Simulasi Rantai Markov: Kebangkrutan Penjudi

Asal komponen dan lisensi. Bagian “Simulasi rantai” menerjemahkan serta mengadaptasi baris 601–666 dari source/05-Markov-chains.Rmd karya Gordan Žitković pada komit e2b35ad91a3689454ae6455e8ffc510a90760c0d. Donornya dilepas dengan CC0 1.0 Universal; terjemahan dan adaptasi bahasa Indonesia ini dilepas dengan CC BY 4.0. Benih deterministik, pembanding peluang eksak, dan bagian yang berlabel Tambahan asli merupakan materi baru yang juga dilepas dengan CC BY 4.0.

Simulasi rantai

Matriks transisi memakai urutan keadaan \(S=(0,1,2,3)\). Baris pertama dan terakhir menyatakan keadaan menyerap. Fungsi langkah_berikutnya memilih baris yang sesuai dengan keadaan sekarang, lalu sample menarik keadaan berikutnya dari peluang pada baris itu. Potongan berikut hanya memakai R dasar dan menetapkan benih secara eksplisit.

S = c(0L, 1L, 2L, 3L)
P = matrix(c(1,   0,   0,   0,
             1/3, 0,   2/3, 0,
             0,   1/3, 0,   2/3,
             0,   0,   0,   1),
           byrow=TRUE, ncol=4)

horizon = 100L
nsim = 1000L
seed = 20260822L

langkah_berikutnya = function(s) {
  i = match(s, S)
  sample(S, prob=P[i, ], size=1)
}

lintasan_tunggal = function(keadaan_awal) {
  lintasan = integer(horizon)
  terakhir = keadaan_awal
  for (n in seq_len(horizon)) {
    lintasan[n] = langkah_berikutnya(terakhir)
    terakhir = lintasan[n]
  }
  lintasan
}

set.seed(seed)
lintasan = t(replicate(nsim, lintasan_tunggal(1L)))
berhasil = lintasan[, horizon] == 3L

distribusi = c(0, 1, 0, 0)
for (n in seq_len(horizon)) {
  distribusi = as.numeric(distribusi %*% P)
}
peluang_eksak_horizon = distribusi[4]
peluang_akhir = (1 - (1/2)^1) / (1 - (1/2)^3)
celah_ekor_eksak = (4/7) * (2/9)^50

hasil = data.frame(
  seed = seed,
  simulasi = nsim,
  horizon = horizon,
  keadaan_awal = 1L,
  batas_atas = 3L,
  berhasil = sum(berhasil),
  taksiran = sprintf("%.12f", mean(berhasil)),
  eksak_sampai_horizon = sprintf("%.12f", peluang_eksak_horizon),
  peluang_akhir = sprintf("%.12f", peluang_akhir),
  celah_ekor_eksak = sprintf("%.12e", celah_ekor_eksak),
  galat_mutlak = sprintf("%.12f", abs(mean(berhasil) - peluang_eksak_horizon)),
  stringsAsFactors = FALSE,
  check.names = FALSE
)
write.csv(hasil, file="", row.names=FALSE, quote=FALSE)

Karena keadaan \(0\) dan \(3\) menyerap, syarat lintasan[, horizon] == 3 ekuivalen dengan “telah mencapai \(3\) paling lambat pada waktu \(T\).” Kolom eksak_sampai_horizon dihitung dengan memperbarui distribusi keadaan melalui perkalian matriks sebanyak \(100\) kali. Kolom peluang_akhir adalah peluang mencapai \(3\) sebelum \(0\) tanpa batas waktu. Pada horizon genap \(100\), peluang yang belum terserap dan kelak mencapai \(3\) adalah \((4/7)(2/9)^{50}=1{,}248349703776\times10^{-33}\). Kolom celah_ekor_eksak mempertahankan bilangan positif yang sangat kecil ini dalam notasi ilmiah; eksak_sampai_horizon dan peluang_akhir tetap tampil sama pada 12 angka desimal. Dalam aritmetika titik-mengambang, iterasi matriks dan representasi \(4/7\) justru berbeda sekitar \(2{,}22\times10^{-16}\) akibat pembulatan—jauh lebih besar daripada celah eksak di atas—sehingga selisih bilangan mesin itu tidak dapat dipakai untuk mengukur celah horizon.

Tambahan asli: persamaan harmonik dan pembanding eksak