擲出 20 面骰子一回合擊敗龍的概率
在數學課上，我們被問到一個我自己無法解決的可選問題：
你正在與一條 250 點生命值的龍戰鬥，並且正在滾動一個 20 面骰子來造成傷害。龍受到的傷害等於擲出的骰子數，但有 2 個例外：如果您擲出 20，則不會造成任何傷害（龍對致命傷害免疫），如果您擲出 1，則破壞您的連擊，龍會立即殺死您。你滾動直到龍被擊敗，或者你擲出 1 並打破你的組合。打敗龍的概率是多少？
我的數學方法包括計算每卷的預期傷害$$ E_{ w }= \frac{ \sum_{ n=2 }^{ w-1 }{ n } }{ w } $$在哪裡$$ E_{ 20 }= 9.45 $$然後嘗試計算在擊敗龍所需的最少擲骰數中沒有失敗的概率，但最終導致我的概率如此之低，任何人都會認為它為 0，而且一開始似乎不正確。
我確實明白這個問題比看起來要復雜得多，因為概率隨著任何其他滾動之前的滾動數量及其結果而發生巨大變化。
通過一個小python程序的經驗測試，結果出來的概率只有25%左右，這讓我很驚訝，因為造成傷害的機會是90%。
您將如何正確解決問題，甚至將其“參數化”，以便可以更改龍的 HP 和骰子上的面數？
您解決問題的一般版本的建議是正確的。 讓我們設置一下。
die 有兩個特殊的結果：“death”，終止進程，0，沒有任何效果。我們不妨去掉 0，創建一個 19 面的“截斷骰子”。讓骰子出現一個數值的概率 $ \omega $ 是 $ p(\omega) $ 然後讓 $ p_{*} $ 是死亡的概率。 $ \Omega $ 是所有這些可能數值的集合（不包括非數值的“死亡”）。
您的目標是總共達到 $ T $ 在觀察“死亡”之前。什麼時候 $ T\le 0, $ 你已經達到了這個門檻，所以獲勝的機會是 $ 1. $ 否則，當 $ T \gt 0, $ 將“最終我贏”事件劃分為發生在其中的單獨編號的結果 $ \Omega. $ 這是一個概率公理，該事件的機會是其（非重疊）組件的機會之和。
$$ \Pr(\text{Reach } T) = \sum_{\omega\in\Omega}p(\omega) \Pr(\text{Reach } T-\omega). $$
這種遞歸可以用一個簡單形式的動態程序來執行。 除非值和概率非常特殊，否則您不能希望為解決方案編寫一個很好的封閉公式。您只需要通過計算值來執行計算 $ T=1, $ $ T=2, $ *等等，*按順序。（這在計算機科學中稱為“消除尾遞歸” 。）
該算法執行的計算次數與 $ T $ 乘以截斷骰子上唯一值的數量。這使得它對中等規模的人有效 $ T $ 和現實的骰子。
通過這樣的程序（使用雙精度浮點），我找到了達到的機會 $ T=250 $ 是 $ 0.269880432506\ldots. $
作為一個現實檢查，你預計每卷會造成大約 9.5 點傷害，這表明它需要大約 $ 250/9.5 \approx 27 $ 滾動獲勝。但是每卷都有一個 $ 1-1/20 $ 生存的機會，所以到那時你的生存機會是$$ (1-1/20)^{27} = \left[(1-1/20)^{20}\right]^{27/20} \approx \exp(-27/20) = 0.25924\ldots. $$
這與我得到的答案非常接近。
作為另一個現實檢查，這也接近您的模擬結果。實際上，我獲得了可比較的模擬結果：它們與確切答案沒有顯著差異。
我把它留給你寫程序。 它將需要一個可以存儲所有值的數據結構 $ \Pr(\text{Reach } T-\omega) $ 在公式的右側給出。考慮為此使用一個數組，由值索引 $ 0,1,2,\ldots, T. $
順便說一句，還有其他解決方法。 這個問題描述了一個馬爾可夫鏈，它的狀態是已經達到的總值（從 $ 0 $ 通過 $ T, $ 因為任何大於 $ T $ 還不如結合 $ T $ )，以及一種特殊的（吸引人的）“死亡”狀態。這個鏈可以用一個大矩陣來分析（有 $ 250+2 $ 方面）。作為一個實際問題，這個公式並不值錢，但馬爾可夫鏈理論提供了對該過程的洞察。您可以挖掘該理論以獲取有關您獲勝機會的信息，以及如果您獲勝，您可能會贏多少次。
在對該問題的評論中提出了**另一種方法：利用幾何分佈。 這是指根據你在死前將有多少卷來分析這個過程。為了處理生命值，想像一下擲骰子，將其“死亡”面與翻轉（偏向）硬幣平行移除，其功能是確定您是否死亡。（因此，在問題的情況下，骰子剩下的 19 個面中的每一個——包括 $ 0, $ 必須留在裡面——有一個 $ 1/19 $ 機會; 更一般地說，有價值的一方 $ \omega $ 有機會 $ p(\omega)/(1-p_{}). $ ）硬幣的兩面是“死”（有概率 $ p_{} $ ）和“繼續”（有概率 $ 1-p_{*} $ ）。在每個回合中，您分別滾動截斷的骰子並擲硬幣，累積生命值直到達到閾值 $ T $ 或者硬幣變成“死亡”。
**可以進行簡化，**因為很容易計算出一開始就不會翻轉“死亡”的機會 $ n=0,1,2,\ldots $ 轉：等於 $ (1-p_{})^n $ 因為所有的翻轉都是獨立的。形式上，這描述了一個隨機變量 $ N $ 其值等於 $ n $ 有概率 $ p_{}(1-p_{*})^{n}. $ （這是一個幾何分佈）。
為了模擬截斷模具的滾動，讓 $ X_1 $ 是它在第一卷中的價值， $ X_2 $ 它在第二次滾動中的值，依此類推。之後的總和 $ n $ 卷因此是 $ S_n = X_1 + X_2 + \cdots + X_n. $ （這是一個隨機遊走。）可以通過將此事件分解為對應於所需滾動數的可數無窮大可能性來計算達到閾值的機會。條件概率的基本規則告訴我們
$$ \Pr(\text{Reach }T) = \sum_{n=0}^\infty \Pr(S_N\ge T\mid N=n)\Pr(N=n).\tag{*} $$
右手邊要求我們找到每個人的這些機會 $ S_n $ 達到閾值。雖然這對於一般模具來說並沒有太大的簡化（ $ S_n $ 可以有一個非常複雜的分佈），當過程在死亡或達到閾值之前可能需要多次滾動時，它會導致一個很好的近似值。那是因為大量的總和 $ X_i $ 近似具有正態分佈（根據中心極限定理）。當截斷骰子的期望是 $ \mu $ （等於 $ 9.45/(1-0.05) $ 在問題中），它的方差是 $ \sigma^2, $ 的分佈 $ S_n $ 有一個期望 $ n\mu $ 和方差 $ n\sigma^2 $ （根據適用於自變量的期望和方差的基本規律 $ X_1,X_2,\ldots, X_n $ ）。寫作 $ \Phi(x;n\mu,n\sigma^2 $ 對於具有期望的正態分佈函數 $ n\mu $ 和方差 $ n\sigma^2, $ 我們獲得
$$ \Pr(S_N\ge T\mid N=n) \approx 1 - \Phi\left(T-\frac{1}{2}; n\mu, n\sigma^2\right). $$
將其插入 $ (*) $ 與幾何分佈規律一起產生
$$ \Pr(\text{Reach }T) \approx \sum_{n=1}^\infty \left(1 - \Phi\left(T-\frac{1}{2}; n\mu, n\sigma^2\right)\right)p_{}(1-p_{})^n. $$
實際上，我們可能會在 $ \sum_{i=n}^\infty p_{}(1-p_{})^i $ 小於一個可容忍的錯誤 $ \epsilon\gt 0, $ 因為 $ \Phi $ 係數永遠不會超過 $ 1. $ 這個上限等於 $ \log(p_{}\epsilon)/\log(1-p_{}). $ （對於問題中的情況，上限約為 418。）我們還可以計算出一個合理的值來開始求和（通過跳過非常小的初始值）。這導致了下面顯示的相對較短和簡單的代碼（用 編寫
R）。它的輸出，通過命令獲得dragon.Normal()，是 $ 0.269879\ldots, $ 同意五個有效數字的確切答案。# Use `NA` to specify the "death" sides of the die; otherwise, specify the # values on its faces. `p` gives the associated probabilities. dragon.Normal <- function(threshold=250, die=c(NA, 2:19, 0), p=rep(1/20,20), eps=1e-6) { # # Find and remove the "death" face(s) to create the truncated die. # i.death <- which(is.na(die)) p.death <- sum(p[i.death]) if (p.death <= 0) return(1) die <- die[-i.death] p <- p[-i.death] p <- p / sum(p) # # Compute the expectation and variance of the truncated die. # mu <- sum(die * p) sigma2 <- sum((die-mu)^2 * p) # # Establish limits for the sum. # N <- ceiling(log(eps * p.death) / log(1 - p.death)) if (N > 1e8) stop("Problem is too large.") Z <- qnorm(eps) n <- min(N, max(1, floor((Z * (1 - sqrt(1 + 4*threshold*mu/(Z^2*sigma2)))/(2*mu))^2 * sigma2))) # # Compute the sum. # n <- n:N sum(p.death * (1-p.death)^n * pnorm(threshold - 1/2, mu*n, sqrt(sigma2*n), lower.tail=FALSE)) }