蒙特卡洛模拟抛硬币得到一定的模式
受到这篇文章的启发： 投币模式的统计 ，我已经进行了一个蒙特卡罗模拟，以确定抛硬币的预期数量，通过使用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 程序的输出如下所示：
这与文中所显示的结果是一致的。 抛硬币的其他模式也是匹配的。 尽pipe如此，我仍然不确定我写的程序是否正确。 当硬币是不公平的时候，我怀疑是什么意思， p1 ， p2和p3不等于0.5，因为我没有任何信息去检查它的准确性。 我也想知道如何用VBA Excel或R编写一个高效的程序来执行上面的模拟，比如THTHTHTHT，THTTHHTHTTH等更长的模式，其循环超过100万（可能是1亿或10亿），但是还是相当快的？ 任何想法？
为了使它更有效率，可以通过分配一个位来使用variables的位。 然后，对于每一次抛球，旋转左侧的比特，并在第一个位置添加折腾结果，直到右侧的比特与图案匹配：
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
您将需要一个string来存储您想要查找的模式。
然后在每次折腾后将最新结果追加到结果string的结尾。
然后检查结果string的最后n个数字==模式，其中n =模式的长度。
如果匹配，然后logging扔数和空白结果string，然后再去…
你大概可以做20行左右的代码！ 希望有所帮助。