excel - 用于抛硬币以获得特定模式的蒙特卡罗模拟

Question

受本文启发:Statistics of Coin-Toss Patterns ，我进行了蒙特卡罗模拟，以确定通过使用 Excel VBA 获得特定模式所需的抛硬币次数。下面的代码是蒙特卡洛模拟，用于抛掷一枚公平的硬币得到模式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

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

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

下面是一个获取匹配序列的平均 throw 次数的示例:

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