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

Question

Inspired by this article: Statistics of Coin-Toss Patterns, I have conducted a Monte Carlo simulation for determining the expected number of tossing a coin to get a certain pattern by using Excel VBA. The following code is the Monte Carlo simulation for tossing a fair coin to get pattern HTH, where H is head (1) and T is tail (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

The output of the program is as shown below:

which agrees with the result as shown in the article. Other patterns for tossing a fair coin are also match. Despite the results, I'm still feeling uncertain whether the program I wrote is correct or not. My doubt arises when the coin is unfair, meaning p1, p2, and p3 are not equal to 0.5, since I don't have any information to check its accuracy. I also want to know how to write an efficient program in VBA Excel or R to perform simulation above for a longer pattern like THTHTHTHT, THTTHHTHTTH, etc. and its looping is more than 1,000,000 (maybe 100,000,000 or 1,000,000,000) but still pretty fast? Any idea?

Answer 1

To make it more efficient, you could use the bits of a variable by assigning a bit with a toss. Then for each toss, rotate the bits on the left and add the toss result at the first position until the bits on the right are a match with the pattern :

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

Note that VBA doesn't have a bit shift operator, but shifting 1 bit on the left is the same as multiplying by 2 :

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

Here is an example to get the average number of toss to match a sequence:

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