Visual Basic for Applications/VBA 的 PRNG

摘要

偽隨機數生成器 (PRNG) 如果執行足夠長的時間，會生成一個基於其演算法的特徵序列。這個序列永遠重複，而且是不變的。VBA 的 Rnd() 函式，如果在迴圈中使用，不帶引數，也不使用 Randomize()，將生成 16,777,216 個介於零和一之間的值，然後從頭開始，生成一個週期長度為 16,777,216 的重複序列。使用者唯一的選擇是選擇該序列中的起點。這是透過選擇一個開始值或種子來完成的。在 Rnd() 的情況下，開始值透過兩種方式選擇：預設情況下，使用系統計時器，或者使用使用者設定的數字。同樣，使用者設定的開始引數不會生成新的序列，只會決定該序列中的哪一部分將被使用。線性同餘生成器 (LCG)，即 Microsoft 的 Rnd() 函式使用的型別，在 線性同餘生成器 中有詳細描述。
單級 LCG 的最大週期長度等於其模數。對於組合生成器，最大週期長度等於各個生成器的週期長度的最小公倍數。設計良好的生成器將具有最大週期長度，並且在其整個序列中只包含唯一的值，但並非所有生成器都是設計良好的。上面的連結描述了設計 LCG 所需的值，以便其在所有起始值上都具有最大週期長度。
下面的程式碼模組包含 VBA 中的 Wichmann-Hill (1982) CLCG (組合 LCG)，並且完全可用。它被稱為 RndX()，並與它自己的 RandomizeX() 配合使用。它比 Microsoft 的 Rnd() 具有更長的重複週期。下面給出了 RndX() 的最有用設定的摘要，以及那些需要更多細節的人在下拉框中的額外資訊。遺憾的是，作者缺乏進行任何嚴肅的數字生成器測試的工具和知識，因此以下內容可能只對初學者感興趣。
長週期生成器在 Excel 中難以研究。Microsoft 的 Rnd() 和使用者函式 RndX() 的週期都太長，無法將完整週期寫入單個工作表列中。解決方案是隻列出長流的一部分，或者建立一個週期足夠短以供完整列出的數字生成器。以這種方式在單個列中列出，可以確認重複週期的長度，然後在刪除重複項後修剪行，然後計算刪除重複項後的行數，可以向懷疑者確認所有值都是唯一的。模組中包含了一些過程，用於列出 Wichmann-Hill 實現的一部分，它適合大約 30269 行左右，以及另一個非常簡單的生成器，用於進一步測試，它只適合 43 行。

微軟的 Rnd() 演算法

微軟的 Visual Basic for Applications (VBA)，目前在 Rnd() 函式中使用線性同餘生成器 (LCG) 進行偽隨機數生成。由於溢位，嘗試在 VBA 中實現 Microsoft 演算法失敗。以下是其基本演算法。

      x1 = ( x0 * a + c ) MOD m
  and;       
      Rnd() = x1/m
  where:
      Rnd() = returned value
      m = modulus = (2^24)
      x1 = new value
      x0 = previous value (initial value 327680)
      a = 16598013
      c = 12820163
      Repeat length = m = (2^24) = 16,777,216

可以注意到 Microsoft 的 Rnd() 與下面由 Wichmann-Hill (1982) 描述的演算法之間的相似之處，其中三個 LCG 表示式的總和用於生成每個輸出數字。表示式的組合給出了 RndX()，其編碼值具有其大大改善的迴圈長度，為

      Cycle length = least_common_multiple(30268, 30306, 30322) = 30268 * 30306 * 30322 / 4 = 6,953,607,871,644

VBA 程式碼 - Wichmann-Hill (1982)

關於模組級變數的提醒可能是有序的。模組級變數在過程執行之間保持其值。實際上，它們將保留值，直到 VBA 不再使用或程式碼被編輯。程式碼中已加入這些變數的重置，以確保從預期的值開始，而不是從之前的頂級過程執行的舊儲存值開始。

注意：雖然該演算法比駐留的 Rnd() 具有改進的屬性，但執行這些生成器的應用程式本身並不特別安全。還要考慮，如果起始值已知，則所有 LCG 編碼的輸出都是完全可預測的。實際上，如果已知該流的任何部分，那麼那些有意傷害的人可以透過將其與儲存的值進行比較來找到整個流。這些事實加在一起限制了此類 VBA 實現的使用範圍，僅限於學習或非關鍵應用程式。

話雖如此，以下可能是最實用的引數配置：在每種情況下，RandomizeX() 應該只調用一次，在呼叫 RndX() 之前，並且在任何包含 RndX() 的生成器迴圈之外。此建議也適用於 Microsoft 函式 Rnd() 及其配套函式 Randomize()。

要生成具有不可預測起始點的輸出，以及每次執行時不同的起始點
- 在呼叫 RndX 之前，呼叫 RandomizeX，不帶任何引數。這使用了系統計時器。
要從大量可重複的起始點生成輸出，並且由使用者引數選擇
- 在呼叫 RndX 之前，呼叫 RandomizeX，使用任何數字引數。更改 RandomizeX 引數值會導致標準演算法流的不同起始點。
要生成不可預測的單個值，每次執行時都不同
- 在呼叫 RndX 之前，呼叫 RandomizeX，不帶任何引數，並使用引數為零。這使用了系統計時器。
要生成可重複的單個值，與使用者引數相關，並由使用者引數選擇
- 在呼叫 RndX 之前，呼叫 RandomizeX，使用任何數字引數，並使用引數為零。更改 RandomizeX 引數值會導致不同的值，這些值對每個引數都是特有的。
請參考下面的下拉框，獲取引數設定及其結果的完整列表。

PRNG RndX() 和 RandomizeX() 引數詳細資訊

RndX() 和 RandomizeX() 引數詳細資訊
RandomizeX() 引數	RndX() 引數	函式的行為 (假設編碼是用於生成序列)
無。	無	PRNG 流由計算機系統計時器的執行時取樣決定。流不確定。
無	正數	PRNG 流由計算機系統計時器的執行時取樣決定。流不確定。RndX() 的正引數對其沒有影響。
無	負數	一個數字，可重複，每個數字都不同，並且取決於RndX() 引數的值。例如：RndX(-3) 導致 0.05079271
無	零	一個數字，可重複，由計算機系統計時器的執行時取樣決定；例如：序列為 0.1741…，01741…
數字²	無	PRNG 流，可重複，每個數字都不同，並且取決於RandomizeX() 引數的值。
數字	正數	PRNG 流，可重複，每個數字都不同，並且取決於RandomizeX() 引數的值。RndX() 的正引數對其沒有影響。
數字	負數	一個數字，可重複，每個數字都不同，並且取決於 RndX() 引數的值。RandomizeX() 引數值沒有任何影響。例如：RndX(-51) 導致 0.8634…
數字	零	一個數字，可重複，每個數字都不同，並且取決於RandomizeX() 引數的值。例如：RandomizeX(2346) 導致 0.2322…
函式未使用	無	預設 PRNG 流，可重複，始終相同。例如：序列為 0.8952…，0.1114…，0.9395…
函式未使用	正數	預設 PRNG 流，可重複，始終相同。例如：序列為 0.8952…，0.1114…，0.9395…
函式未使用¹	負數或零	一個數字，可重複，每個數字都不同，並且取決於RndX() 引數的值。例如：RndX(0) = 0.8694...: -5 = 0.0846…

1. 術語函式未使用是指函式不是由使用者在程式碼中專門呼叫。在某些情況下，例如這種情況，RandomizeX() 函式仍然需要在程式碼中可用，以便RndX() 函式的內部呼叫。
2. 數字項是可以轉換為數字的項。RandomizeX() 函式使用其變體引數中給定的種子值生成一個正整數。它對字串的任何前導部分也是如此，直到第一個無法識別為數字的字元為止。

本節中的程式碼應另存為 Excel 中的單獨標準模組。

Option Explicit
Dim nSamples As Long
Dim nX As Long, nY As Long, nZ As Long

Sub TestRndX()
    'run this to obtain RndX() samples
    'Wichmann, Brian; Hill, David (1982), Algorithm AS183:
    'An Efficient and Portable Pseudo-Random Number Generator,
    'Journal of the Royal Statistical Society. Series C
    Dim n As Long
   
    'reset module variables
    nX = 0: nY = 0: nZ = 0
    
    RandomizeX
    For n = 1 To 10
        Debug.Print RndX()
        MsgBox RndX()
    Next n
   
    'reset module variables
    nX = 0: nY = 0: nZ = 0

End Sub

Sub TestScatterChartOfPRNG()
    'run this to make a point scatter chart
    'using samples from RndX
    
    Dim vA As Variant, n As Long
    Dim nS As Long, nR As Double
    
    'remove any other charts
    'DeleteAllCharts
    
    'reset module variables
    nX = 0: nY = 0: nZ = 0
    
    'set number of samples here
    nSamples = 1000
    ReDim vA(1 To 2, 1 To nSamples) 'dimension array
        
    'load array with PRNG samples
    RandomizeX
    For n = 1 To nSamples
        nR = RndX()
        vA(1, n) = n  'x axis data - sample numbers
        vA(2, n) = nR 'y axis data - prng values
    Next n
    
    'make scatter point chart from array
    ChartScatterPoints vA, 1, 2, nSamples & " Samples of RndX()", _
                "Sample Numbers", "PRNG Values [0,1]"
    
    'reset module work variables
    nX = 0: nY = 0: nZ = 0

End Sub

Sub RandomizeX(Optional ByVal nSeed As Variant)
   'sets variables for PRNG procedure RndX()
      
   Const MaxLong As Double = 2 ^ 31 - 1
   Dim nS As Long
   Dim nN As Double
   
   'make multiplier
   If IsMissing(nSeed) Then
      nS = Timer * 60
   Else
      nN = Abs(Int(Val(nSeed)))
      If nN > MaxLong Then 'no overflow
         nN = nN - Int(nN / MaxLong) * MaxLong
      End If
      nS = nN
   End If
   
   'update variables
   nX = (nS Mod 30269)
   nY = (nS Mod 30307)
   nZ = (nS Mod 30323)
   
   'avoid zero state
   If nX = 0 Then nX = 171
   If nY = 0 Then nY = 172
   If nZ = 0 Then nZ = 170

End Sub

Function RndX(Optional ByVal nSeed As Long = 1) As Double
   'PRNG - gets pseudo random number - use with RandomizeX
   'Wichmann-Hill algorithm of 1982
   
   Dim nResult As Double
   
   'initialize variables
   If nX = 0 Then
      nX = 171
      nY = 172
      nZ = 170
   End If
   
   'first update variables
   If nSeed <> 0 Then
      If nSeed < 0 Then RandomizeX (nSeed)
      nX = (171 * nX) Mod 30269
      nY = (172 * nY) Mod 30307
      nZ = (170 * nZ) Mod 30323
   End If
   
   'use variables to calculate output
   nResult = nX / 30269# + nY / 30307# + nZ / 30323#
   RndX = nResult - Int(nResult)

End Function

Sub ChartScatterPoints(ByVal vA As Variant, RowX As Long, RowY As Long, _
                     Optional sTitle As String = "", Optional sXAxis As String, _
                     Optional sYAxis As String)
    
    'array input must contain two data rows for x and y data
    'parameters for user title, x axis and y axis labels
    'makes a simple point scatter chart
    
    Dim LBC As Long, UBC As Long, LBR As Long, UBR As Long, n As Long, bOptLim As Boolean
    Dim X As Variant, Y As Variant, sX As String, sY As String, sT As String, oC As Chart
    
    LBR = LBound(vA, 1): UBR = UBound(vA, 1)
    LBC = LBound(vA, 2): UBC = UBound(vA, 2)
    ReDim X(LBC To UBC)
    ReDim Y(LBC To UBC)
    
    'labels for specific charts
    If sTitle = "" Then sT = "Title Goes Here" Else sT = sTitle
    If sXAxis = "" Then sX = "X Axis Label Goes Here" Else sX = sXAxis
    If sYAxis = "" Then sY = "Y Axis Label Goes Here" Else sY = sYAxis
    
    If RowX < LBR Or RowX > UBR Or RowY < LBC Or RowY > UBC Then
        MsgBox "Parameter data rows out of range in ChartColumns - closing"
        Exit Sub
    End If
    
    'transfer data to chart arrays
    For n = LBC To UBC
        X(n) = vA(RowX, n) 'x axis data
        Y(n) = vA(RowY, n) 'y axis data
    Next n
    
    'make chart
    Charts.Add
    
    'set chart type
    ActiveChart.ChartType = xlXYScatter 'point scatter chart
        
    'remove unwanted series
    With ActiveChart
        Do Until .SeriesCollection.Count = 0
            .SeriesCollection(1).Delete
        Loop
    End With
    
    
    'assign the data and labels to a series
    With ActiveChart.SeriesCollection
        If .Count = 0 Then .NewSeries
            If Val(Application.Version) >= 12 Then
                .Item(1).Values = Y
                .Item(1).XValues = X
            Else
                .Item(1).Select
                Names.Add "_", X
                ExecuteExcel4Macro "series.x(!_)"
                Names.Add "_", Y
                ExecuteExcel4Macro "series.y(,!_)"
                Names("_").Delete
            End If
    End With
        
    'apply title string, x and y axis strings, and delete legend
    With ActiveChart
        .HasTitle = True
        .ChartTitle.Text = sT
        .SetElement (msoElementPrimaryCategoryAxisTitleAdjacentToAxis) 'X
        .Axes(xlCategory).AxisTitle.Text = sX
        .SetElement (msoElementPrimaryValueAxisTitleRotated) 'Y
        .Axes(xlValue).AxisTitle.Text = sY
        .Legend.Delete
    End With
    
    'trim axes to suit
    With ActiveChart
    'X Axis
        .Axes(xlCategory).Select
        .Axes(xlCategory).MinimumScale = 0
        .Axes(xlCategory).MaximumScale = nSamples
        .Axes(xlCategory).MajorUnit = 500
        .Axes(xlCategory).MinorUnit = 100
        Selection.TickLabelPosition = xlLow
        
    'Y Axis
        .Axes(xlValue).Select
        .Axes(xlValue).MinimumScale = -0.2
        .Axes(xlValue).MaximumScale = 1.2
        .Axes(xlValue).MajorUnit = 0.1
        .Axes(xlValue).MinorUnit = 0.05
    End With
    
    
    ActiveChart.ChartArea.Select
    
    Set oC = Nothing

End Sub

Sub DeleteAllCharts5()
    'run this to delete all ThisWorkbook charts
    
    Dim oC
       
    Application.DisplayAlerts = False
    
    For Each oC In ThisWorkbook.Charts
        oC.Delete
    Next oC
    
    Application.DisplayAlerts = True
    
End Sub

PRNG 的簡單測試

下面的程式碼模組包含 Wichmann-Hill (1982) 演算法的簡化版本，實際上只使用了其三個計算部分中的第一個。它將使用不同的起始值在執行它的工作簿的Sheet1 上生成多個完整的數值流。請注意，第一個值在第 30269 行都重複出現，如果擴充套件，整個流也會重複出現。生成列表後，使用電子表格的列排序和刪除重複項函式來檢視每列是否包含適當數量的唯一條目。還包含了一個更簡單的生成器，其重複週期僅為 43，這可能使研究更易於管理，並且可以透過執行TestMSRnd 來檢視 Microsoft 的Rnd() 的週期在 16777216 (+1) 處重複。

本節中的程式碼應另存為 Excel 中的單獨標準模組。

Option Explicit

Private ix2 As Long

Sub TestWHRnd30269()
    'makes five columns for complete output streams
    'each with a different start point
    'runs a simplified LCNG with mod 30269
        
    Dim sht As Worksheet, nS As Double, nSamp As Long
    Dim c As Long, r As Long, nSeed As Long
    
    'set seed value for Rnd2()
    nSeed = 327680 'WH initial seed
    
    'set number of random samples to make
    nSamp = 30275 '30269 plus say, 6
    
    'set initial value of carry variable
    ix2 = nSeed
    
    Set sht = ThisWorkbook.Worksheets("Sheet1")
    
    'clear the worksheet
    sht.Cells.Cells.ClearContents
    
    'load sheet with set of samples
    For c = 1 To 5                'number of runs
        ix2 = nSeed + c           'change start value
        For r = 1 To nSamp        'number of samples
            nS = WHRnd30269()     'get a sample
            sht.Cells(r, c) = nS  'write to sheet
        Next r
    Next c

    sht.Cells(1, 1).Select

End Sub

Function WHRnd30269() As Double
   'first part of Wichmann-Hill tripple.
   'When started with seed ix2 = 171,
   'full sequence repeats from n = 30269
   'without any repeated values before.
   
   Dim r As Double
   
   'ix2 cannot be 0.
   If ix2 = 0 Then
      ix2 = 171
   End If
   
   'calculate Xn+1 from Xn
   ix2 = (171 * ix2) Mod 30269
   
   'make an output value
   r = ix2 / 30269#
   WHRnd30269 = r - Int(r)

End Function

Sub TestSimpleRnd43()
    'makes five columns for complete output streams
    'each with a different start point
    'runs a very simple LCNG with mod 43
        
    Dim sht As Worksheet, nS As Double, nSamp As Long
    Dim c As Long, r As Long, nSeed As Long
    
    'set seed value for Rnd2()
    nSeed = 17 'initial seed
    
    'set number of random samples to make
    nSamp = 45 '43 plus say, 2
    
    'set initial value of carry variable
    ix2 = nSeed
    
    Set sht = ThisWorkbook.Worksheets("Sheet1")
    
    'clear the worksheet
    sht.Cells.Cells.ClearContents
    
    'load sheet with set of samples
    For c = 1 To 5                'number of runs
        ix2 = nSeed + c           'change start value
        For r = 1 To nSamp        'number of samples
            nS = SimpleRnd43()    'get a sample
            sht.Cells(r, c) = nS  'write to sheet
        Next r
    Next c

    sht.Cells(1, 1).Select

End Sub

Function SimpleRnd43() As Double
   'simple Lehmer style LCNG to show repeat streams
   'produces one sequence of 42 unique values - then repeats entire sequence
   'start value decides only where the predictable sequence begins
   
   Dim r As Double
   
   'Note; Makes 42 unique values before sequence repeats
   'Modulus = 43: Multiplier = 5: Initial Seed = 17
   '43 is prime
   '5 is primitive root mod 43
   '17 is coprime to 43
   
   'ix2 cannot be 0.
   If ix2 = 0 Then
      ix2 = 17
   End If
   
   'calculate a new carry variable
   ix2 = (5 * ix2) Mod 43
   
   'make an output value
   r = ix2 / 43#
   SimpleRnd43 = r - Int(r)

End Function

Sub TestMSRnd()
    'makes two sets of single data using MS Rnd
    'the first 10 samples of Rnd() values
    'followed by values around sample 16777216
    'confirms sequence probably re-starts at M+1 = 16777217
    
    Dim sht As Worksheet, nS As Double
    Dim c As Long, r As Long, nMod As Long
    
    'note modulus
    nMod = 16777216
    
    Set sht = ThisWorkbook.Worksheets("Sheet1")
    
    'clear the worksheet
    sht.Cells.Cells.ClearContents
    
    'load sheet with set of samples
    
        For r = 1 To nMod + 20   'number of samples
            nS = Rnd()            'get a sample
            Select Case r
                Case 1 To 10
                    sht.Cells(r, 1) = r
                    sht.Cells(r, 2) = nS
                Case (nMod - 4) To (nMod + 5)
                    sht.Cells(r - 16777211 + 10, 1) = r
                    sht.Cells(r - 16777211 + 10, 2) = nS
            End Select
        Next r
    
    sht.Cells(1, 1).Select

End Sub

參考文獻

Wichmann, Brian; Hill, David (1982), Algorithm AS183: An Efficient and Portable Pseudo-Random Number Generator, Journal of the Royal Statistical Society. Series C

另請參閱

Wichmann-Hill CLCG: 關於該特定組合線性同餘生成器的維基百科文章。
線性同餘生成器: 維基百科對完全週期條件的良好描述。

外部連結

Visual Basic 如何為 RND 函式生成偽隨機數: Microsoft kb231847 知識庫條目