Monte Carlo 模拟抛硬币获得特定模式

Question

受这篇文章的启发：Statistics of Coin-Toss Patterns，我进行了一个Monte Carlo模拟，使用Excel [=31]来确定抛硬币的预期次数以获得某种模式=].下面的代码是Monte Carlo模拟抛一个公平的硬币得到模式HTH，其中H是头（1），T是尾（0）。

Sub Tossing_Coin()
    Dim Toss(1000000) As Double, NToss(1000000) As Double, AVToss(1000000) As Double
    t0 = Timer
    Sheet2.Cells.Clear
    a = 0

    For j = 1 To 1000000

        p1 = Rnd()
        If p1 <= 0.5 Then
            Toss(1) = 1
        Else
            Toss(1) = 0
        End If

        p2 = Rnd()
        If p2 <= 0.5 Then
            Toss(2) = 1
        Else
            Toss(2) = 0
        End If

        i = 2
        Do
            p3 = Rnd()
            If p3 <= 0.5 Then
                Toss(i + 1) = 1
            Else
                Toss(i + 1) = 0
            End If

            i = i + 1
        Loop Until Toss(i - 2) = 1 And Toss(i - 1) = 0 And Toss(i) = 1

        NToss(j) = i
        a = a + NToss(j)
        AVToss(j) = a / j
        b = AVToss(j)
    Next j

    MsgBox "The expected number of tossing is " & b & "." _
        & vbNewLine & "The running time of simulation is " & Round(Timer - t0, 2) & " s."
End Sub

程序输出如下图：

与文章中显示的结果一致。其他用于抛硬币的模式也是匹配的。尽管有结果，我仍然不确定我写的程序是否正确。当硬币不公平时，我会产生疑问，这意味着 p1、p2 和 p3 不等于 0.5，因为我没有任何信息来检查其准确性。我还想知道如何在 VBA Excel 或 R 中编写一个高效的程序来执行上面的模拟，以获得更长的模式，如 TTHTHTHTHT、THTTHHTHTTH 等，它的循环超过 1,000,000（可能是 100,000,000 或1,000,000,000）但仍然很快？有什么想法吗？

Answer 1

您需要一个字符串来存储您要查找的模式。

然后在每次投掷后将最新结果附加到结果字符串的末尾。

然后检查结果字符串的最后 n 位是否 == 模式，其中 n = 模式的长度。

如果匹配则记录投掷次数和空白结果字符串并再次...

您大概可以在大约 20 行代码内完成！希望对您有所帮助。

Answer 2

为了提高效率，您可以通过抛出分配一个位来使用变量的位。然后对于每次投掷，旋转左侧的位并将投掷结果添加到第一个位置，直到右侧的位与模式匹配：

pattern "HTH"  : 101
mask for "XXX" : 111

1 toss "H" :       1 And 111 = 001
2 toss "T" :      10 And 111 = 010
3 toss "T" :     100 And 111 = 100
4 toss "H" :    1001 And 111 = 001
5 toss "H" :   10011 And 111 = 011
6 toss "T" :  100110 And 111 = 110
7 toss "H" : 1001101 And 111 = 101  : "HTH" matches the first 3 bits

注意 VBA 没有位移运算符，但向左移动 1 位与乘以 2 相同：

  decimal  9 =   1001 in bits
 9 +  9 = 18 =  10010 in bits
18 + 18 = 36 = 100100 in bits

下面是一个获取匹配序列的平均投掷次数的示例：

Sub UsageExample()
    Const sequence = "HTH"
    Const samples = 100000

    MsgBox "Average: " & GetTossingAverage(sequence, samples)
End Sub

Function GetTossingAverage(sequence As String, samples As Long) As Double
    Dim expected&, tosses&, mask&, tossCount#, i&
    Randomize ' Initialize the random generator. '

    ' convert the [TH] sequence to a sequence of bits. Ex: HTH -> 00000101 '
    For i = 1 To Len(sequence)
        expected = expected + expected - (Mid$(sequence, i, 1) = "T")
    Next

    ' generate the mask for the rotation of the bits. Ex: HTH -> 01110111 '
    mask = (2 ^ (Len(sequence) * 2 + 1)) - (2 ^ Len(sequence)) - 1

    ' iterate the samples '
    For i = 1 To samples
        tosses = mask

        ' generate a new toss until we get the expected sequence '
        Do
            tossCount = tossCount + 1
            ' rotate the bits on the left and rand a new bit at position 1 '
            tosses = (tosses + tosses - (Rnd < 0.5)) And mask
        Loop Until tosses = expected
    Next

    GetTossingAverage = tossCount / samples
End Function

Monte Carlo 模拟抛硬币获得特定模式

Monte Carlo simulation for tossing a coin to get a certain pattern

excel

vba

montecarlo

probability-theory