Diharapkan nomor saya akan aktif setelah menggambar kartu sampai saya mendapatkan kartu as, 2, 3, dan sebagainya

Saya mengalami beberapa masalah dalam menyelesaikan yang berikut ini.

Anda mengambil kartu dari tumpukan kartu 52 kartu standar tanpa penggantian sampai Anda mendapatkan kartu As. Anda menarik dari apa yang tersisa sampai Anda mendapatkan 2. Anda melanjutkan dengan 3. Berapa angka yang diharapkan Anda akan berada setelah seluruh dek habis?

Wajar membiarkannya

• $T_i = \text{first position of card whose value is }i$
• $U_i = \text{last position of card whose value is }i$

Jadi masalah pada dasarnya sama dengan mencari tahu probabilitas bahwa Anda akan berada di ketika dek habis, yaitu:$k$

$$Pr(T_1<\cdots<T_k \cap U_{k+1} < T_k)$$

Saya bisa melihatnya

$$Pr(T_1<\cdots<T_k) = 1/k! \\ \text{and} \\ Pr(U_{k+1} < T_k) = 1/70$$

tetapi tidak bisa melangkah lebih jauh ...

mengikuti ide @ gung, saya percaya nilai yang diharapkan adalah 5,84? dan dari interpretasi saya tentang komentar, saya berasumsi "A" adalah nilai yang hampir mustahil (kecuali empat kartu terakhir di geladak adalah kartu As). berikut adalah hasil dari iterasi 100.000 simulasi monte carlo

results
    2     3     4     5     6     7     8     9     J     K     Q     T 
 1406  7740 16309 21241 19998 15127  9393  4906   976   190   380  2334 

dan inilah kode R jika Anda ingin bermain dengannya ..

# monte carlo card-drawing functions from here
# http://streaming.stat.iastate.edu/workshops/r-intro/lectures/5-Rprogramming.pdf

# create a straightforward deck of cards
create_deck <-
    function( ){
        suit <- c( "H" , "C" , "D" , "S" )
        rank <- c( "A" , 2:9 , "T" , "J" , "Q" , "K" )
        deck <- NULL
        for ( r in rank ) deck <- c( deck , paste( r , suit ) )
        deck
    }

# construct a function to shuffle everything
shuffle <- function( deck ){ sample( deck , length( deck ) ) }

# draw one card at a time
draw_cards <-
    function( deck , start , n = 1 ){
        cards <- NULL

        for ( i in start:( start + n - 1 ) ){
            if ( i <= length( deck ) ){
                cards <- c( cards , deck[ i ] )
            }
        }

        return( cards )
    }

# create an empty vector for your results
results <- NULL

# run your simulation this many times..
for ( i in seq( 100000 ) ){
    # create a new deck
    sdeck <- shuffle( create_deck() )

    d <- sdeck[ grep('A|2' , sdeck ) ]
    e <- identical( grep( "2" , d ) , 1:4 )

    # loop through ranks in this order
    rank <- c( "A" , 2:9 , "T" , "J" , "Q" , "K" )

    # start at this position
    card.position <- 0

    # start with a blank current.draw
    current.draw <- ""

    # start with a blank current rank
    this.rank <- NULL

    # start with the first rank
    rank.position <- 1

    # keep drawing until you find the rank you wanted
    while( card.position < 52 ){

        # increase the position by one every time
        card.position <- card.position + 1

        # store the current draw for testing next time
        current.draw <- draw_cards( sdeck , card.position )

        # if you draw the current rank, move to the next.
        if ( grepl( rank[ rank.position ] , current.draw ) ) rank.position <- rank.position + 1

        # if you have gone through every rank and are still not out of cards,
        # should it still be a king?  this assumes yes.
        if ( rank.position == length( rank ) ) break        

    }

    # store the rank for this iteration.
    this.rank <- rank[ rank.position ]

    # at the end of the iteration, store the result
    results <- c( results , this.rank )

}

# print the final results
table( results )

# make A, T, J, Q, K numerics
results[ results == 'A' ] <- 1
results[ results == 'T' ] <- 10
results[ results == 'J' ] <- 11
results[ results == 'Q' ] <- 12
results[ results == 'K' ] <- 13
results <- as.numeric( results )

# and here's your expected value after 100,000 simulations.
mean( results )

Untuk simulasi, sangat penting untuk menjadi benar dan cepat. Kedua tujuan ini menyarankan penulisan kode yang menargetkan kemampuan inti dari lingkungan pemrograman serta kode yang sesingkat dan sesederhana mungkin, karena kesederhanaan memberikan kejelasan dan kejelasan mempromosikan kebenaran. Inilah usaha saya untuk mencapai keduanya di R:

#
# Simulate one play with a deck of `n` distinct cards in `k` suits.
#
sim <- function(n=13, k=4) {
  deck <- sample(rep(1:n, k)) # Shuffle the deck
  deck <- c(deck, 1:n)        # Add sentinels to terminate the loop
  k <- 0                      # Count the cards searched for
  for (j in 1:n) {
    k <- k+1                          # Count this card
    deck <- deck[-(1:match(j, deck))] # Deal cards until `j` is found
    if (length(deck) < n) break       # Stop when sentinels are reached
  }
  return(k)                   # Return the number of cards searched
}

Menerapkan ini dengan cara yang dapat direproduksi dapat dilakukan dengan replicatefungsi setelah mengatur seed number acak, seperti pada

> set.seed(17);  system.time(d <- replicate(10^5, sim(13, 4)))
   user  system elapsed 
   5.46    0.00    5.46

Itu lambat, tetapi cukup cepat untuk melakukan simulasi yang cukup panjang (dan karenanya tepat) berulang kali tanpa menunggu. Ada beberapa cara untuk menunjukkan hasilnya. Mari kita mulai dengan artinya:

> n <- length(d)
> mean(d)
[1] 5.83488

> sd(d) / sqrt(n)
[1] 0.005978956

Yang terakhir adalah kesalahan standar: kami berharap rata-rata yang disimulasikan berada dalam dua atau tiga SE dari nilai sebenarnya. Itu menempatkan harapan sebenarnya di suatu tempat antara dan$5.817$$5.853$ .

Kami juga mungkin ingin melihat tabulasi frekuensi (dan mereka kesalahan standar). Kode berikut sedikit memberi sedikit tabulasi:

u <- table(d)
u.se <- sqrt(u/n * (1-u/n)) / sqrt(n)
cards <- c("A", "2", "3", "4", "5", "6", "7", "8", "9", "T", "J", "Q", "K")
dimnames(u) <- list(sapply(dimnames(u), function(x) cards[as.integer(x)]))
print(rbind(frequency=u/n, SE=u.se), digits=2)

Berikut hasilnya:

                2       3      4      5      6      7       8       9       T       J       Q       K
frequency 0.01453 0.07795 0.1637 0.2104 0.1995 0.1509 0.09534 0.04995 0.02249 0.01009 0.00345 0.00173
SE        0.00038 0.00085 0.0012 0.0013 0.0013 0.0011 0.00093 0.00069 0.00047 0.00032 0.00019 0.00013

---

Bagaimana kita tahu simulasi itu benar? Salah satu caranya adalah mengujinya secara mendalam untuk masalah yang lebih kecil. Untuk alasan itu kode ini ditulis untuk menyerang generalisasi kecil dari masalah, mengganti kartu yang berbeda dengan dan kartu dengan . Namun, untuk pengujian penting untuk dapat memberi makan kode dek dalam urutan yang telah ditentukan. Mari kita tulis antarmuka yang sedikit berbeda dengan algoritma yang sama:$13$n$4$k

draw <- function(deck) {
  n <- length(sentinels <- sort(unique(deck)))
  deck <- c(deck, sentinels)
  k <- 0
  for (j in sentinels) {
    k <- k+1
    deck <- deck[-(1:match(j, deck))]
    if (length(deck) < n) break
  }
  return(k)
}

(Dimungkinkan untuk digunakan drawdi simmana - mana, tetapi kerja ekstra yang dilakukan di awal drawmembuatnya dua kali lebih lambat sim.)

Kita dapat menggunakan ini dengan menerapkannya pada setiap pengocokan berbeda dari dek yang diberikan. Karena tujuan di sini hanya beberapa tes satu kali, efisiensi dalam menghasilkan shuffles tidak penting. Berikut ini cara kasar yang cepat:

n <- 4 # Distinct cards
k <- 2 # Number of suits
d <- expand.grid(lapply(1:(n*k), function(i) 1:n))
e <- apply(d, 1, function(x) var(tabulate(x))==0)
g <- apply(d, 1, function(x) length(unique(x))==n)
d <- d[e & g,]

Sekarang dadalah bingkai data yang barisnya berisi semua shuffles. Terapkan drawke setiap baris dan hitung hasilnya:

d$result <- apply(as.matrix(d), 1, draw)
    (counts <- table(d$result))

Output (yang akan kita gunakan dalam tes formal sebentar lagi) adalah

   2    3    4 
 420  784 1316 

(Nilai mudah dimengerti, ngomong-ngomong: kita masih akan mengerjakan kartu jika dan hanya jika semua pasangan mendahului semua kartu As. Peluang terjadinya ini (dengan dua setelan) adalah . Dari berbeda, miliki properti ini.)2 1 / ( 2 + $420$$2$25202520/6=$1/\binom{2+2}{2} = 1/6$$2520$$2520/6 = 420$

Kami dapat menguji output dengan uji chi-squared. Untuk tujuan ini saya menerapkan kali untuk kasus ini kartu berbeda dalam :sim n = 4 k = $10,000$$n = 4$$k = 2$

>set.seed(17)
>d.sim <- replicate(10^4, sim(n, k))
>print((rbind(table(d.sim) / length(d.sim), counts / dim(d)[1])), digits=3)

         2     3     4
[1,] 0.168 0.312 0.520
[2,] 0.167 0.311 0.522

> chisq.test(table(d.sim), p=counts / dim(d)[1])

    Chi-squared test for given probabilities

data:  table(d.sim) 
X-squared = 0.2129, df = 2, p-value = 0.899

Karena sangat tinggi, kami tidak menemukan perbedaan yang signifikan antara apa yang dikatakan dan nilai-nilai yang dihitung oleh enumerasi lengkap. Mengulangi latihan ini untuk beberapa nilai dan lainnya (kecil) menghasilkan hasil yang sebanding, memberi kami banyak alasan untuk percaya ketika diterapkan pada dan .n k n = 13 k = $p$sim$n$$k$sim$n=13$$k=4$

Akhirnya, uji chi-kuadrat dua sampel akan membandingkan keluaran simke keluaran yang dilaporkan dalam jawaban lain:

>y <- c(1660,8414,16973,21495,20021,14549,8957,4546,2087,828,313,109)
>chisq.test(cbind(u, y))

data:  cbind(u, y) 
X-squared = 142.2489, df = 11, p-value < 2.2e-16

Statistik chi-squared yang sangat besar menghasilkan nilai-p yang pada dasarnya nol: tanpa keraguan, simtidak setuju dengan jawaban lainnya. Ada dua kemungkinan resolusi perselisihan: satu (atau keduanya!) Jawaban ini salah atau mereka menerapkan interpretasi yang berbeda dari pertanyaan tersebut. Sebagai contoh, saya telah menafsirkan "setelah dek habis" berarti setelah mengamati kartu terakhir dan, jika diperbolehkan, memperbarui "nomor Anda akan berada di" sebelum mengakhiri prosedur. Dapat dibayangkan bahwa langkah terakhir tidak dimaksudkan untuk diambil. Mungkin beberapa perbedaan penafsiran yang halus seperti itu akan menjelaskan ketidaksepakatan, pada titik mana kita dapat memodifikasi pertanyaan untuk memperjelas apa yang ditanyakan.

Ada jawaban yang tepat (dalam bentuk produk matriks, disajikan pada poin 4 di bawah). Algoritma yang cukup efisien untuk menghitungnya ada, berasal dari pengamatan ini:

1. Acak acak kartu dapat dihasilkan dengan mengocok kartu secara acak dan kemudian secara acak memotong kartu tersisa di dalamnya.N $N+k$$N$$k$

2. Dengan hanya mengocok kartu As, dan kemudian (menerapkan pengamatan pertama) menyelingi keduanya, lalu bertiga, dan seterusnya, masalah ini dapat dilihat sebagai rantai dari tiga belas langkah.

3. Kita perlu melacak lebih dari nilai kartu yang kita cari. Ketika melakukan ini, kita tidak perlu memperhitungkan posisi tanda relatif terhadap semua kartu, tetapi hanya posisinya relatif terhadap kartu yang nilainya sama atau lebih kecil.

   Bayangkan menempatkan tanda pada kartu as pertama, dan kemudian menandai dua yang pertama ditemukan setelahnya, dan seterusnya. (Jika pada suatu tahap dek habis tanpa menampilkan kartu yang sedang kami cari, kami akan membiarkan semua kartu tidak ditandai.) Biarkan "tempat" dari setiap tanda (bila ada) adalah jumlah kartu dengan nilai yang sama atau lebih rendah yang dibagikan ketika tanda dibuat (termasuk kartu yang ditandai itu sendiri). Tempat-tempat berisi semua informasi penting.

4. Tempat setelah tanda dibuat adalah angka acak. Untuk dek tertentu, urutan tempat-tempat ini membentuk proses stokastik. Ini sebenarnya adalah proses Markov (dengan matriks transisi variabel). Karena itu jawaban yang tepat dapat dihitung dari dua belas perkalian matriks.$i^\text{th}$

Menggunakan ide-ide ini, mesin ini memperoleh nilai (komputasi dalam floating point presisi ganda) dalam detik. Perkiraan nilai persis ini akurat untuk semua digit yang ditampilkan.1 / 9 $5.8325885529019965$$1/9$

$$\frac{1982600579265894785026945331968939023522542569}{339917784579447928182134345929899510000000000}$$

Sisa dari posting ini memberikan perincian, menyajikan implementasi kerja (dalam R), dan diakhiri dengan beberapa komentar tentang pertanyaan dan efisiensi solusi.

---

Menghasilkan serutan acak dari sebuah geladak

Sebenarnya lebih jelas secara konseptual dan tidak lebih rumit secara matematis untuk mempertimbangkan "dek" (alias multiset ) dari kartu yang ada dari denominasi terendah, dari terendah berikutnya, dan seterusnya . (Pertanyaan yang diajukan menyangkut dek yang ditentukan oleh vektor- .)$N = k_1+k_2+\cdots+k_m$$k_1$$k_2$$13$$(4,4,\ldots,4)$

A "acak acak" kartu adalah satu permutasi diambil secara seragam dan acak dari permutasi kartuShuffles ini jatuh ke dalam kelompok konfigurasi yang setara karena "ace" antara mereka sendiri tidak mengubah apa pun, "dua" antara mereka sendiri juga tidak mengubah apa pun, dan sebagainya. Oleh karena itu setiap kelompok permutasi yang terlihat identik ketika kartu-kartu tersebut diabaikan berisipermutasi. Kelompok-kelompok ini, yang jumlahnya diberikan oleh koefisien multinomial$N$$N! = N\times(N-1)\times\cdots\times 2\times 1$$N$$k_1$$k_2$$k_1!\times k_2!\times \cdots \times k_m!$

$$\binom{N}{k_1,k_2,\ldots,k_m} = \frac{N!}{k_1!k_2!\cdots k_m!},$$

disebut "kombinasi" dari dek.

Ada cara lain untuk menghitung kombinasi. Kartu pertama dapat membentuk kombinasi. Mereka meninggalkan "slot" di antara dan di sekelilingnya tempat kartu berikutnya dapat ditempatkan. Kami dapat menunjukkan ini dengan diagram di mana " " menunjuk salah satu kartu dan " " menunjuk sebuah slot yang dapat menampung antara dan kartu tambahan:$k_1$$k_1!/k_1! = 1$$k_1+1$$k_2$$*$$k_1$$\_$$0$$k_2$

$$\underbrace{\_*\_*\_\cdots\_*\_}_{k_1\text{ stars}}$$

Ketika kartu tambahan diselingi, pola bintang dan kartu baru kartu menjadi dua subset. Jumlah himpunan bagian yang berbeda adalah .$k_2$$k_1+k_2$$\binom{k_1+k_2}{k_1,k_2} = \frac{(k_1+k_2)!}{k_1!k_2!}$

Mengulangi prosedur ini dengan "bertiga," kami menemukan ada cara untuk menyelinginya di antara kartu pertama . Karenanya jumlah total cara berbeda untuk mengatur kartu dengan cara ini sama dengan$k_3$$\binom{(k_1+k_2)+k_3}{k_1+k_2,k_3}= \frac{(k_1+k_2+k_3)!}{(k_1+k_2)!k_3!}$$k_1+k_2$$k_1+k_2+k_3$

$$1\times\frac{(k_1+k_2)!}{k_1!k_2!}\times\frac{(k_1+k_2+k_3)!}{(k_1+k_2)!k_3!} = \frac{(k_1+k_2+k_3)!}{k_1!k_2!k_3!}.$$

Setelah menyelesaikan kartu terakhir dan terus melipatgandakan fraksi teleskop ini, kami menemukan bahwa jumlah kombinasi berbeda yang diperoleh sama dengan jumlah total kombinasi yang dihitung sebelumnya, . Karenanya, kami tidak mengabaikan kombinasi. Itu berarti proses berurutan mengocok kartu dengan benar menangkap probabilitas setiap kombinasi, dengan asumsi bahwa pada setiap tahap, setiap cara yang mungkin berbeda untuk menyelingi kartu-kartu baru di antara kartu yang lama diambil dengan probabilitas yang sama merata.$k_n$$\binom{N}{k_1,k_2,\ldots,k_m}$

Proses tempat

Awalnya, ada ace dan jelas yang pertama ditandai. Pada tahap selanjutnya ada , tempatnya (jika kartu yang ditandai ada) sama dengan (beberapa nilai dari hingga ), dan kami akan menyelingi kartu sekitar mereka. Kita dapat memvisualisasikan ini dengan diagram seperti$k_1$$n = k_1 + k_2 + \cdots + k_{j-1}$$p$$1$$n$$k=k_j$

$$\underbrace{\_*\_*\_\cdots\_*\_}_{p-1\text{ stars}}\odot\underbrace{\_*\_\cdots\_*\_}_{n-p\text{ stars}}$$

di mana " " menunjukkan simbol yang saat ini ditandai. Bersyarat pada nilai tempat , kami ingin menemukan probabilitas bahwa tempat berikutnya akan sama dengan (beberapa nilai dari hingga ; menurut aturan permainan, tempat berikutnya harus datang setelah , dari mana ). Jika kita dapat menemukan berapa banyak cara yang ada untuk menyelingi kartu baru di tempat kosong sehingga tempat berikutnya sama dengan , maka kita dapat membagi dengan jumlah total cara untuk menyelingi kartu-kartu ini (sama dengan , seperti yang telah kita lihat) untuk mendapatkan$\odot$$p$$q$$1$$n+k$$p$$q\ge p+1$$k$$q$$\binom{n+k}{k}$probabilitas transisi bahwa tempat berubah dari ke . (Akan ada juga kemungkinan transisi untuk tempat tersebut hilang sama sekali ketika tidak ada kartu baru yang mengikuti kartu yang ditandai, tetapi tidak perlu menghitung ini secara eksplisit.)$p$$q$

Mari kita perbarui diagram untuk mencerminkan situasi ini:

$$\underbrace{\_*\_*\_\cdots\_*\_}_{p-1\text{ stars}}\odot\underbrace{**\cdots*}_{s\text{ stars}}\ \vert\ \underbrace{\_*\_\cdots\_*\_}_{n-p-s\text{ stars}}$$

Bilah vertikal " " menunjukkan di mana kartu baru pertama muncul setelah kartu yang ditandai: karena itu tidak ada kartu baru yang muncul di antara dan (dan karenanya tidak ada slot yang ditampilkan dalam interval itu). Kita tidak tahu berapa banyak bintang dalam interval ini, jadi saya baru saja menyebutnya (yang mungkin nol) tidak diketahui akan hilang begitu kita menemukan hubungan antara itu dan .$\vert$$\odot$$\vert$$s$$s$$q$

Misalkan, kemudian, kami menyelingi kartu baru di sekitar bintang sebelum dan then-- secara independen dari yang --Kami menyelingi sisa kartu baru di sekitar bintang setelah . Ada$j$$\odot$$k-j-1$$\vert$

$$\tau_{n,k}(s,p) = \binom{(p-1)+j}{j}\binom{(n-p-s) + (k-j)-1}{k-j-1}$$

cara untuk melakukan ini. Perhatikan, meskipun - ini adalah bagian tersulit dari analisis - bahwa tempat sama dengan karena$\vert$$p+s+j+1$

• Ada kartu "lama" di atau sebelum tanda.$p$
• Ada kartu lama setelah tanda tapi sebelum .$s$$\vert$
• Ada kartu baru sebelum tanda.$j$
• Ada kartu baru yang diwakili oleh itu sendiri.$\vert$

Dengan demikian, memberi kami informasi tentang transisi dari tempat ke tempat . Ketika kami melacak informasi ini dengan cermat untuk semua nilai yang mungkin dari , dan menjumlahkan semua kemungkinan (terpisah) ini, kami memperoleh probabilitas bersyarat tempat berikut tempat ,$\tau_{n,k}(s,p)$$p$$q=p+s+j+1$$s$$q$$p$

$${\Pr}_{n,k}(q|p) = \left(\sum_j \binom{p-1+j}{j}\binom{n+k-q}{k-j-1}\right) / \binom{n+k}{k}$$

di mana jumlah dimulai pada dan berakhir pada . (Panjang variabel dari jumlah ini menunjukkan ada tidak mungkin menjadi formula tertutup untuk itu sebagai fungsi dari dan , kecuali dalam kasus khusus.)$j=\max(0, q-(n+1))$$j=\min(k-1, q-(p+1)$$n, k, q,$$p$

Algoritma

Awalnya ada probabilitas bahwa tempat itu akan menjadi dan probabilitas itu akan memiliki nilai lain yang mungkin dalam . Ini dapat diwakili oleh vektor .$1$$1$$0$$2, 3, \ldots, k_1$$p_1 = (1, 0, \ldots, 0)$

Setelah kartu berikutnya , vektor diperbarui ke dengan mengalikannya (di sebelah kiri) dengan matriks transisi . Ini diulangi sampai semua telah ditempatkan. Pada setiap tahap , jumlah entri dalam vektor probabilitas adalah kemungkinan beberapa kartu telah ditandai. Apa pun yang tersisa untuk membuat nilai sama dengan oleh karena itu adalah kesempatan bahwa tidak ada kartu yang ditandai setelah langkah$k_2$$p_1$$p_2$$(\Pr_{k_1,k_2}(q|p), 1\le p\le k_1, 1\le q\le k_2)$$k_1+k_2+\cdots+k_m$$j$$p_j$$1$$j$. Perbedaan berturut-turut dalam nilai-nilai ini karena itu memberi kita probabilitas bahwa kita tidak dapat menemukan kartu tipe untuk ditandai: yaitu distribusi probabilitas dari nilai kartu yang kita cari ketika tumpukan kartu habis di akhir permainan. .$j$

---

Penerapan

RKode berikut mengimplementasikan algoritma. Ini sejajar dengan diskusi sebelumnya. Pertama, perhitungan probabilitas transisi dilakukan oleh t.matrix(tanpa normalisasi dengan pembagian dengan , membuatnya lebih mudah untuk melacak perhitungan saat menguji kode):$\binom{n+k}{k}$

t.matrix <- function(q, p, n, k) {
  j <- max(0, q-(n+1)):min(k-1, q-(p+1))
  return (sum(choose(p-1+j,j) * choose(n+k-q, k-1-j))
}

Ini digunakan transitionuntuk memperbarui ke . Ini menghitung matriks transisi dan melakukan perkalian. Ia juga menangani perhitungan vektor awal jika argumennya adalah vektor kosong:$p_{j-1}$$p_j$$p_1$p

#
# `p` is the place distribution: p[i] is the chance the place is `i`.
#
transition <- function(p, k) {
  n <- length(p)
  if (n==0) {
    q <- c(1, rep(0, k-1))
  } else {
    #
    # Construct the transition matrix.
    #
    t.mat <- matrix(0, nrow=n, ncol=(n+k))
    #dimnames(t.mat) <- list(p=1:n, q=1:(n+k))
    for (i in 1:n) {
      t.mat[i, ] <- c(rep(0, i), sapply((i+1):(n+k), 
                                        function(q) t.matrix(q, i, n, k)))
    }
    #
    # Normalize and apply the transition matrix.
    #
    q <- as.vector(p %*% t.mat / choose(n+k, k))
  }
  names(q) <- 1:(n+k)
  return (q)
}

Kita sekarang dapat dengan mudah menghitung probabilitas non-mark pada setiap tahap untuk setiap dek:

#
# `k` is an array giving the numbers of each card in order;
# e.g., k = rep(4, 13) for a standard deck.
#
# NB: the *complements* of the p-vectors are output.
#
game <- function(k) {
  p <- numeric(0)
  q <- sapply(k, function(i) 1 - sum(p <<- transition(p, i)))
  names(q) <- names(k)
  return (q)
}

Ini untuk dek standar:

k <- rep(4, 13)
names(k) <- c("A", 2:9, "T", "J", "Q", "K")
(g <- game(k))

Outputnya adalah

         A          2          3          4          5          6          7          8          9          T          J          Q          K 
0.00000000 0.01428571 0.09232323 0.25595013 0.46786622 0.66819134 0.81821790 0.91160622 0.96146102 0.98479430 0.99452614 0.99818922 0.99944610

Menurut aturan, jika seorang raja ditandai maka kita tidak akan mencari kartu lebih lanjut: ini berarti nilai harus ditingkatkan menjadi . Setelah melakukan itu, perbedaannya memberikan distribusi "nomor Anda akan ketika dek habis":$0.9994461$$1$

> g[13] <- 1; diff(g)
          2           3           4           5           6           7           8           9           T           J           Q           K 
0.014285714 0.078037518 0.163626897 0.211916093 0.200325120 0.150026562 0.093388313 0.049854807 0.023333275 0.009731843 0.003663077 0.001810781

(Bandingkan ini dengan keluaran yang saya laporkan dalam jawaban terpisah yang menggambarkan simulasi Monte-Carlo: semuanya tampak sama, hingga jumlah variasi acak yang diharapkan.)

Nilai yang diharapkan segera:

> sum(diff(g) * 2:13)
[1] 5.832589

Semua mengatakan, ini membutuhkan hanya selusin baris kode yang dapat dieksekusi. Saya telah memeriksanya dengan perhitungan tangan untuk nilai kecil (hingga ). Dengan demikian, jika ada perbedaan antara kode dan analisis masalah sebelumnya, percayakan kode tersebut (karena analisis tersebut mungkin memiliki kesalahan ketik).$k$$3$

---

Catatan

Hubungan dengan urutan lainnya

Ketika ada satu kartu masing-masing, distribusi adalah urutan kebalikan dari seluruh angka:

> 1/diff(game(rep(1,10)))
[1]      2      3      8     30    144    840   5760  45360 403200

Nilai di tempat adalah(mulai dari tempat ). Ini adalah urutan A001048 dalam Ensiklopedia Online Urutan Bilangan Bulat. Dengan demikian, kita mungkin berharap untuk formula tertutup untuk deck dengan konstan (deck "cocok") yang akan menggeneralisasi urutan ini, yang dengan sendirinya memiliki beberapa makna mendalam. (Misalnya, ia menghitung ukuran kelas konjugasi terbesar dalam kelompok permutasi dan juga terkait dengan koefisien trinomial .) (Sayangnya, timbal balik dalam generalisasi untuk biasanya tidak bilangan bulat.)$i$$i! + (i-1)!$$i=1$$k_i$$k\gt 1$

Permainan sebagai proses stokastik

Analisis kami memperjelas bahwa koefisien awal vektor , , adalah konstan. Misalnya, mari kita lacak output saat memproses setiap kelompok kartu:$i$$p_j$$j\ge i$game

> sapply(1:13, function(i) game(rep(4,i)))

[[1]]
[1] 0

[[2]]
[1] 0.00000000 0.01428571

[[3]]
[1] 0.00000000 0.01428571 0.09232323

[[4]]
[1] 0.00000000 0.01428571 0.09232323 0.25595013

...

[[13]]
 [1] 0.00000000 0.01428571 0.09232323 0.25595013 0.46786622 0.66819134 0.81821790 0.91160622 0.96146102 0.98479430 0.99452614 0.99818922 0.99944610

Sebagai contoh, nilai kedua dari vektor final (menggambarkan hasil dengan setumpuk penuh 52 kartu) sudah muncul setelah kelompok kedua diproses (dan sama dengan ). Dengan demikian, jika Anda menginginkan informasi hanya tentang tanda naik melalui nilai kartu , Anda hanya perlu melakukan perhitungan untuk setumpuk kartu .$1/\binom{8}{4}=1/70$$j^\text{th}$$k_1+k_2+\cdots+k_j$

Karena peluang untuk tidak menandai kartu nilai semakin cepat mendekati ketika meningkat, setelah jenis kartu dalam empat setelan, kita hampir mencapai nilai pembatas untuk ekspektasi. Memang, nilai pembatas sekitar (dihitung untuk setumpuk kartu , di mana titik kesalahan pembulatan presisi ganda mencegah melangkah lebih jauh).$j$$1$$j$$13$$5.833355$$4 \times 32$

Pengaturan waktu

Melihat algoritma yang diterapkan pada vektor- , kami melihat waktunya harus proporsional dengan dan - menggunakan batas atas mentah - tidak lebih buruk daripada proporsional dengan . Dengan menghitung semua perhitungan untuk hingga dan hingga , dan menganalisis hanya mereka yang mengambil waktu yang relatif lama ( detik atau lebih lama), saya memperkirakan waktu perhitungan sekitar , mendukung penilaian batas atas ini.( k , k , ... , k ) k 2 m 3 k = 1 7 n = 10 30 1 / 2 O ( k 2 n 2,9$m$$(k,k, \ldots, k)$$k^2$$m^3$$k=1$$7$$n=10$$30$$1/2$$O(k^2 n^{2.9})$

Salah satu penggunaan asimptotik ini adalah memproyeksikan waktu perhitungan untuk masalah yang lebih besar. Sebagai contoh, melihat bahwa kasus membutuhkan waktu sekitar detik, kami akan memperkirakan bahwa kasus (sangat menarik) akan memakan waktu sekitar detik. (Sebenarnya butuh detik.)1.31 k = 1 , n = 100 1,31 ( 1 / 4 ) 2 ( 100 / 30 ) 2,9 ≈ 2,7 $k=4, n=30$$1.31$$k=1, n=100$$1.31(1/4)^2(100/30)^{2.9}\approx 2.7$$2.87$

Meretas Monte Carlo sederhana di Perl dan menemukan sekitar .$5.8329$

#!/usr/bin/perl

use strict;

my @deck = (1..13) x 4;

my $N = 100000; # Monte Carlo iterations.

my $mean = 0;

for (my $i = 1; $i <= $N; $i++) {
    my @d = @deck;
    fisher_yates_shuffle(\@d);
    my $last = 0;
        foreach my $c (@d) {
        if ($c == $last + 1) { $last = $c }
    }
    $mean += ($last + 1) / $N;
}

print $mean, "\n";

sub fisher_yates_shuffle {
    my $array = shift;
        my $i = @$array;
        while (--$i) {
        my $j = int rand($i + 1);
        @$array[$i, $j] = @$array[$j, $i];
    }
}