"random" --- 生成偽隨機數
*************************

**原始碼：**Lib/random.py

======================================================================

本章中所提及的 module（模組）用來實現各種分佈的虛擬隨機數產生器。

對於整數，可以從範圍中進行均勻選擇。對於序列，有一個隨機元素的均勻選擇
，一個用來原地 (in-place) 產生隨機排列清單的函式，以及一個用來隨機採樣
不替換的函式。

在實數線上，有一些函式用於處理均勻分佈、常態分佈（高斯分佈）、對數常態
分佈、負指數分佈、gamma 分佈和 Beta 分佈。對於生成角度分佈，可以使用馮
·米塞斯分佈 (von Mises distribution)。

幾乎所有 module 函式都相依於基本函式 "random()"，此函式在半開放範圍
"0.0 <= X < 1.0" 內均勻地生成一個隨機 float（浮點數）。Python 使用
Mersenne Twister（梅森旋轉演算法）作為核心的產生器，它產生 53 位元精度
float，其週期為 2**19937-1，透過 C 語言進行底層的實作既快速又支援執行
緒安全 (threadsafe)。Mersenne Twister 是現存最廣泛被驗證的隨機數產生器
之一，但是基於完全確定性，它並不適合所有目的，並且完全不適合加密目的。

該 module 提供的函式實際上是 "random.Random" class（類別）中一個隱藏實
例的綁定方法 (bound method)。你可以實例化自己的 "Random" 實例，以得到
不共享狀態的產生器。

如果你想使用你自己設計的基本產生器，"Random" 也可以進行子類別化
(subclass)。有關詳細資訊，請參考該類別的文件。

"random" module 也提供了 "SystemRandom" class，使用系統函式
"os.urandom()" 從作業系統提供的來源產生隨機數。

警告:

  本章所提及的虛擬隨機數產生器不應該使用於安全目的。有關安全性或加密用
  途，請參考 "secrets" module。

也參考:

  M. Matsumoto and T. Nishimura, "Mersenne Twister: A
  623-dimensionally equidistributed uniform pseudorandom number
  generator", ACM Transactions on Modeling and Computer Simulation
  Vol. 8, No. 1, January pp.3--30 1998.

  進位互補乘法 (Complementary-Multiply-with-Carry) 用法，可作為隨機數
  產生器的一個可相容替代方案，具有較長的週期和相對簡單的更新操作。


簿記函式 (bookkeeping functions)
================================

random.seed(a=None, version=2)

   初始化隨機數產生器。

   如果 *a* 被省略或為 "None"，則使用當前系統時間。如果隨機來源由作業
   系統提供，則使用它們而不是系統時間（有關可用性的詳細資訊，請參考
   "os.urandom()" 函式）。

   如果 *a* 為 int（整數），則直接使用它。

   如使用版本 2（預設值），"str"、"bytes" 或 "bytearray" 物件將轉換為
   "int"，並使用其所有位元。

   若使用版本 1（為復現於舊版本 Python 中產生隨機序列而提供），"str"
   和 "bytes" 的演算法會生成範圍更窄的種子 (seed)。

   在 3.2 版的變更: 移至版本 2 方案，該方案使用字串種子中的所有位元。

   在 3.11 版的變更: *seed* 必須是以下型別之一："None"、"int"、"float"
   、"str"、"bytes"、"bytearray"。

random.getstate()

   回傳一個物件，捕獲產生器的當前內部狀態。此物件可以傳遞給
   "setstate()" 以恢復狀態。

random.setstate(state)

   *state* 應該要從之前對 "getstate()" 的呼叫中獲得，並且以
   "setstate()" 將產生器的內部狀態恢復到呼叫 "getstate()" 時的狀態。


回傳位元組的函式
================

random.randbytes(n)

   產生 *n* 個隨機位元組。

   此方法不應使用於產生安全性權杖 (Token)。請改用
   "secrets.token_bytes()"。

   在 3.9 版被加入.


回傳整數的函式
==============

random.randrange(stop)
random.randrange(start, stop[, step])

   傳回從 "range(start, stop, step)" 中隨機選擇的元素。

   這大致相當於 "choice(range(start, stop, step))"，但支援任意大的範圍
   ，並針對常見情況進行了最佳化。

   位置引數模式與 "range()" 函式相符。

   不應使用關鍵字引數，因為它們可能會以意想不到的方式被直譯。例如
   "randrange(start=100)" 會被直譯為 "randrange(0, 100, 1)"。

   在 3.2 版的變更: "randrange()" 在產生均勻分佈的值方面更為複雜。以前
   ，它使用像 "int(random()*n)" 這樣的樣式，這可能會產生稍微不均勻的分
   佈。

   在 3.12 版的變更: 已經不再支援非整數類型到等效整數的自動轉換。像是
   "randrange(10.0)" 和 "randrange(Fraction(10, 1))" 的呼叫將會引發
   "TypeError"。

random.randint(a, b)

   回傳一個隨機整數 *N*，使得 "a <= N <= b"。為 "randrange(a, b+1)" 的
   別名。

random.getrandbits(k)

   回傳一個具有 *k* 個隨機位元的非負 Python 整數。此方法會隨 Mersenne
   Twister 產生器一起提供，一些其他的產生器也可能將其作為 API 的可選部
   分。如果可用，"getrandbits()" 使 "randrange()" 能夠處理任意大的範圍
   。

   在 3.9 版的變更: 此方法現在接受 *k* 為零。


回傳序列的函式
==============

random.choice(seq)

   從非空序列 *seq* 回傳一個隨機元素。如果 *seq* 為空，則引發
   "IndexError"。

random.choices(population, weights=None, *, cum_weights=None, k=1)

   回傳從 *population* 中重置取樣出的一個大小為 *k* 的元素 list。如果
   *population* 為空，則引發 "IndexError"。

   如果指定了 *weights* 序列，則根據相對權重進行選擇。另外，如果給定
   *cum_weights* 序列，則根據累積權重進行選擇（可能使用
   "itertools.accumulate()" 計算）。例如，相對權重 "[10, 5, 30, 5]" 等
   同於累積權重 "[10, 15, 45, 50]"。在內部，相對權重在進行選擇之前會轉
   換為累積權重，因此提供累積權重可以節省工作。

   如果既未指定 *weights* 也未指定 *cum_weights*，則以相等的機率進行選
   擇。如果提供了加權序列，則該序列的長度必須與 *population* 序列的長
   度相同。它是一個 "TypeError" 來指定 *weights* 和 *cum_weights*。

   *weights* 或 *cum_weights* 可以使用任何與 "random()" 所回傳的
   "float" 值（包括整數、float 和分數，但不包括小數）交互操作
   (interoperates) 的數值類型。權重假定為非負數和有限的。如果所有權重
   均為零，則引發 "ValueError"。

   對於給定的種子，具有相等權重的 "choices()" 函式通常產生與重複呼叫
   "choice()" 不同的序列。"choices()" 使用的演算法使用浮點運算來實現內
   部一致性和速度。"choice()" 使用的演算法預設為整數運算和重複選擇，以
   避免捨入誤差產生的小偏差。

   在 3.6 版被加入.

   在 3.9 版的變更: 如果所有權重均為零，則引發 "ValueError"。

random.shuffle(x)

   將序列 *x* 原地 (in place) 隨機打亂位置。

   要打亂一個不可變的序列並回傳一個新的被打亂的 list（串列），請使用
   "sample(x, k=len(x))"。

   請注意，即使對於較小的 "len(x)"，*x* 的置換總數也會快速成長到大於大
   多數隨機數產生器的週期。這意味著長序列的大多數置換永遠無法產生。例
   如，長度為 2080 的序列是 Mersenne Twister 隨機數產生器週期內可以容
   納的最大序列。

   在 3.11 版的變更: 移除可選參數 *random*。

random.sample(population, k, *, counts=None)

   回傳從母體序列中選擇出的一個包含獨特元素、長度為 *k* 的 list。用於
   不重置的隨機取樣。

   回傳包含母體元素的新清單，同時保持原始母體不變。產生的清單按選擇順
   序排列，因此所有子切片也會是有效的隨機樣本。這允許抽獎獲獎者（樣本
   ）分為大獎和第二名獲獎者（子切片）。

   母體成員不必是 *hashable* 或唯一的。如果母體包含重複項，則每次出現
   都是樣本中可能出現的一個選擇。

   可以一次指定一個重複元素，也可以使用可選的僅關鍵字 *counts* 參數指
   定重複元素。例如 "sample(['red', 'blue'], counts=[4, 2], k=5)" 等同
   於 "sample(['red', 'red', 'red', 'red', 'blue', 'blue'], k=5)"。

   若要從整數範圍中選擇範例，請使用 "range()" 物件作為引數。這對於從大
   型母體中取樣特別快速且節省空間："sample(range(10000000), k=60)"。

   如果樣本大小大於母體大小，"ValueError" 會被引發。

   在 3.9 版的變更: 新增 *counts* 參數。

   在 3.11 版的變更: *population* 必須是一個序列。不再支援將 set 自動
   轉換為 list。


離散分布
========

以下函式產生離散分佈。

random.binomialvariate(n=1, p=0.5)

   二項分佈 (Binomial distribution)。回傳 *n* 個獨立試驗的成功次數，每
   個試驗的成功機率為 *p*：

   數學上等價於：

      sum(random() < p for i in range(n))

   試驗次數 *n* 應為非負整數。成功的機率 *p* 應在 "0.0 <= p <= 1.0" 之
   間。結果是 "0 <= X <= n" 範圍內的整數。

   在 3.12 版被加入.


實數分布
========

以下函式產生特定的實數分佈。函式參數以分佈方程中的對應變數命名，如常見
的數學實踐所示；這些方程式中的大多數都可以在任意統計文本中找到。

random.random()

   回傳範圍 "0.0 <= X < 1.0" 中的下一個隨機浮點數

random.uniform(a, b)

   回傳一個隨機浮點數 *N*，當 "a <= b" 時確保 N 為 "a <= N <= b" 、"b
   < a" 時確保 N 為 "b <= N <= a"。

   終點值 "b" 可能包含在範圍內，也可能不包含在範圍內，取決於運算式 "a
   + (b-a) * random()" 中的浮點捨入。

random.triangular(low, high, mode)

   回傳一個隨機浮點數 *N*，使得 "low <= N <= high"，並在這些邊界之間具
   有指定的 *mode*。*low* 和 *high* 邊界預設為零和一。*mode* 引數預設
   為邊界之間的中點，從而給出對稱分佈。

random.betavariate(alpha, beta)

   Beta（貝它）分布。參數的條件為 "alpha > 0" 和 "beta > 0"。回傳值的
   範圍介於 0 和 1 之間。

random.expovariate(lambd=1.0)

   指數分佈。*lambd* 為 1.0 除以所需的平均數。它應該不為零。（該參數將
   被稱為 "lambda"，但這是 Python 中的保留字）如果 *lambd* 為正，則回
   傳值的範圍從 0 到正無窮大；如果 *lambd* 為負，則回傳值的範圍從負無
   窮大到 0。

   在 3.12 版的變更: 新增 "lambd" 的預設值。

random.gammavariate(alpha, beta)

   Gamma（伽瑪）分佈。（*不是* Gamma 函式！）。形狀 (shape) 和比例
   (scale) 參數 *alpha* 和 *beta* 必須具有正值。（根據呼叫習慣不同，部
   分來源會將 'beta' 定義為比例的倒數）。

   Probability distribution function（機率密度函式）是：

                x ** (alpha - 1) * math.exp(-x / beta)
      pdf(x) =  --------------------------------------
                  math.gamma(alpha) * beta ** alpha

random.gauss(mu=0.0, sigma=1.0)

   常態分佈，也稱為高斯分佈。*mu* 是平均數，*sigma* 是標準差。這比下面
   定義的 "normalvariate()" 函式快一點。

   多執行緒須注意：當兩個執行緒同時呼叫此函式時，它們可能會收到相同的
   傳回值。這可以透過三種方式避免。1）讓每個執行緒使用隨機數產生器的不
   同實例。2）在所有呼叫周圍加鎖。3）使用較慢但執行緒安全的
   "normalvariate()" 函式代替。

   在 3.11 版的變更: *mu* 和 *sigma* 現在有預設引數。

random.lognormvariate(mu, sigma)

   對數常態分佈。如果你取此分佈的自然對數，你將獲得一個具有平均數 *mu*
   和標準差 *sigma* 的常態分佈。*mu* 可以為任何值，並且 *sigma* 必須大
   於零。

random.normalvariate(mu=0.0, sigma=1.0)

   常態分佈。*mu* 是平均數，*sigma* 是標準差。

   在 3.11 版的變更: *mu* 和 *sigma* 現在有預設引數。

random.vonmisesvariate(mu, kappa)

   *mu* 是平均角度，以 0 到 2**pi* 之間的弧度表示，*kappa* 是濃度參數
   ，必須大於或等於零。如果 *kappa* 等於零，則此分佈在 0 到 2**pi* 的
   範圍內將減小為均勻的隨機角度。

random.paretovariate(alpha)

   Pareto distribution（柏拉圖分佈）。*alpha* 是形狀參數。

random.weibullvariate(alpha, beta)

   Weibull distribution（韋伯分佈）。*alpha* 是比例參數，*beta* 是形狀
   參數。


替代產生器
==========

class random.Random([seed])

   實現 "random" 模組使用的預設偽隨機數產生器的 class。

   在 3.11 版的變更: 過去 *seed* 可以是任何可雜湊物件，但現在必須是以
   下類型之一："None"、"int"、"float"、"str"、"bytes"、"bytearray"。

   如果 "Random" 的子類別希望使用不同的基礎產生器，則應該覆寫以下方法
   ：

   seed(a=None, version=2)

      在子類別中覆寫此方法以自訂 "Random" 實例的 "seed()" 行為。

   getstate()

      在子類別中覆寫此方法以自訂 "Random" 實例的 "getstate()" 行為。

   setstate(state)

      在子類別中覆寫此方法以自訂 "Random" 實例的 "setstate()" 行為。

   random()

      在子類別中覆寫此方法以自訂 "Random" 實例的 "random()" 行為。

   或者，自訂產生器子類別還可以提供以下方法：

   getrandbits(k)

      在子類別中覆寫此方法以自訂 "Random" 實例的 "getrandbits()" 行為
      。

class random.SystemRandom([seed])

   使用 "os.urandom()" 函式從作業系統提供的來源產生隨機數的 Class。並
   非在所有系統上都可用。不依賴於軟體狀態，並且序列不可復現。因此
   "seed()" 方法沒有效果且被忽略。如果呼叫 "getstate()" 和
   "setstate()" 方法會引發 "NotImplementedError"。


關於 Reproducibility（復現性）的注意事項
========================================

有時，能夠重現偽隨機數產生器給出的序列很有用。只要多執行緒未運行，透過
重複使用種子值，同一序列就應該可以被復現。

大多數隨機 module 的演算法和 seed 設定函式在 Python 版本中可能會發生變
化，但可以保證兩個方面不會改變：

* 如果增加了新的 seed 設定函式，則將提供向後相容的播種器 (seeder)。

* 當相容的播種器被賦予相同的種子時，產生器的 "random()" 方法將持續產生
  相同的序列。


範例
====

基礎範例：

   >>> random()                          # Random float:  0.0 <= x < 1.0
   0.37444887175646646

   >>> uniform(2.5, 10.0)                # Random float:  2.5 <= x <= 10.0
   3.1800146073117523

   >>> expovariate(1 / 5)                # Interval between arrivals averaging 5 seconds
   5.148957571865031

   >>> randrange(10)                     # Integer from 0 to 9 inclusive
   7

   >>> randrange(0, 101, 2)              # Even integer from 0 to 100 inclusive
   26

   >>> choice(['win', 'lose', 'draw'])   # Single random element from a sequence
   'draw'

   >>> deck = 'ace two three four'.split()
   >>> shuffle(deck)                     # Shuffle a list
   >>> deck
   ['four', 'two', 'ace', 'three']

   >>> sample([10, 20, 30, 40, 50], k=4) # Four samples without replacement
   [40, 10, 50, 30]

模擬：

   >>> # Six roulette wheel spins (weighted sampling with replacement)
   >>> choices(['red', 'black', 'green'], [18, 18, 2], k=6)
   ['red', 'green', 'black', 'black', 'red', 'black']

   >>> # Deal 20 cards without replacement from a deck
   >>> # of 52 playing cards, and determine the proportion of cards
   >>> # with a ten-value:  ten, jack, queen, or king.
   >>> deal = sample(['tens', 'low cards'], counts=[16, 36], k=20)
   >>> deal.count('tens') / 20
   0.15

   >>> # Estimate the probability of getting 5 or more heads from 7 spins
   >>> # of a biased coin that settles on heads 60% of the time.
   >>> sum(binomialvariate(n=7, p=0.6) >= 5 for i in range(10_000)) / 10_000
   0.4169

   >>> # Probability of the median of 5 samples being in middle two quartiles
   >>> def trial():
   ...     return 2_500 <= sorted(choices(range(10_000), k=5))[2] < 7_500
   ...
   >>> sum(trial() for i in range(10_000)) / 10_000
   0.7958

統計 bootstrapping（自助法）的範例，使用有重置的重新取樣來估計樣本平均
數的信賴區間：

   # https://www.thoughtco.com/example-of-bootstrapping-3126155
   from statistics import fmean as mean
   from random import choices

   data = [41, 50, 29, 37, 81, 30, 73, 63, 20, 35, 68, 22, 60, 31, 95]
   means = sorted(mean(choices(data, k=len(data))) for i in range(100))
   print(f'The sample mean of {mean(data):.1f} has a 90% confidence '
         f'interval from {means[5]:.1f} to {means[94]:.1f}')

重新取樣排列測試的範例，來確定觀察到的藥物與安慰劑之間差異的統計學意義
或 p 值：

   # Example from "Statistics is Easy" by Dennis Shasha and Manda Wilson
   from statistics import fmean as mean
   from random import shuffle

   drug = [54, 73, 53, 70, 73, 68, 52, 65, 65]
   placebo = [54, 51, 58, 44, 55, 52, 42, 47, 58, 46]
   observed_diff = mean(drug) - mean(placebo)

   n = 10_000
   count = 0
   combined = drug + placebo
   for i in range(n):
       shuffle(combined)
       new_diff = mean(combined[:len(drug)]) - mean(combined[len(drug):])
       count += (new_diff >= observed_diff)

   print(f'{n} label reshufflings produced only {count} instances with a difference')
   print(f'at least as extreme as the observed difference of {observed_diff:.1f}.')
   print(f'The one-sided p-value of {count / n:.4f} leads us to reject the null')
   print(f'hypothesis that there is no difference between the drug and the placebo.')

模擬多伺服器佇列 (queue) 的到達時間與服務交付：

   from heapq import heapify, heapreplace
   from random import expovariate, gauss
   from statistics import mean, quantiles

   average_arrival_interval = 5.6
   average_service_time = 15.0
   stdev_service_time = 3.5
   num_servers = 3

   waits = []
   arrival_time = 0.0
   servers = [0.0] * num_servers  # time when each server becomes available
   heapify(servers)
   for i in range(1_000_000):
       arrival_time += expovariate(1.0 / average_arrival_interval)
       next_server_available = servers[0]
       wait = max(0.0, next_server_available - arrival_time)
       waits.append(wait)
       service_duration = max(0.0, gauss(average_service_time, stdev_service_time))
       service_completed = arrival_time + wait + service_duration
       heapreplace(servers, service_completed)

   print(f'Mean wait: {mean(waits):.1f}   Max wait: {max(waits):.1f}')
   print('Quartiles:', [round(q, 1) for q in quantiles(waits)])

也參考:

  Statistics for Hackers 是由 Jake Vanderplas 製作的教學影片，僅使用幾
  個基本概念（包括模擬、取樣、洗牌、交叉驗證）進行統計分析。

  Economics Simulation是由 Peter Norvig 對市場進行的模擬，顯示了該模組
  提供的許多工具和分佈（高斯、均勻、樣本、 beta 變數、選擇，三角形、隨
  機數）的有效使用。

  機率的具體介紹（使用Python）為 Peter Norvig 的教學課程，涵蓋了機率理
  論的基礎知識與如何模擬以及使用 Python 執行數據分析。


使用方案
========

這些使用方案展示了如何有效地從 "itertools" 模組的組合疊代器
(combinatoric iterators) 中進行隨機選擇：

   def random_product(*args, repeat=1):
       "Random selection from itertools.product(*args, **kwds)"
       pools = [tuple(pool) for pool in args] * repeat
       return tuple(map(random.choice, pools))

   def random_permutation(iterable, r=None):
       "Random selection from itertools.permutations(iterable, r)"
       pool = tuple(iterable)
       r = len(pool) if r is None else r
       return tuple(random.sample(pool, r))

   def random_combination(iterable, r):
       "Random selection from itertools.combinations(iterable, r)"
       pool = tuple(iterable)
       n = len(pool)
       indices = sorted(random.sample(range(n), r))
       return tuple(pool[i] for i in indices)

   def random_combination_with_replacement(iterable, r):
       "Choose r elements with replacement.  Order the result to match the iterable."
       # Result will be in set(itertools.combinations_with_replacement(iterable, r)).
       pool = tuple(iterable)
       n = len(pool)
       indices = sorted(random.choices(range(n), k=r))
       return tuple(pool[i] for i in indices)

預設的 "random()" 回傳 *0.0 ≤ x < 1.0* 範圍內 2⁻⁵³ 的倍數。所有數字都
是均勻分佈的，並且可以完全表示為 Python float。但是，該間隔中的許多其
他可表示的 float 不是可能的選擇。 例如 "0.05954861408025609" 不是 2⁻⁵³
的整數倍。

以下範例採用不同的方法。間隔中的所有 float 都是可能的選擇。尾數來自
*2⁵² ≤ 尾數 < 2⁵³* 範圍內的整數均勻分佈。指數來自幾何分佈，其中小於
*-53* 的指數的出現頻率是下一個較大指數的一半。

   from random import Random
   from math import ldexp

   class FullRandom(Random):

       def random(self):
           mantissa = 0x10_0000_0000_0000 | self.getrandbits(52)
           exponent = -53
           x = 0
           while not x:
               x = self.getrandbits(32)
               exponent += x.bit_length() - 32
           return ldexp(mantissa, exponent)

Class 中的所有實數分佈都將使用新方法：

   >>> fr = FullRandom()
   >>> fr.random()
   0.05954861408025609
   >>> fr.expovariate(0.25)
   8.87925541791544

該範例在概念上等效於一種演算法，該演算法從 *0.0 ≤ x < 1.0* 範圍內
2⁻¹⁰⁷⁴ 的所有倍數中進行選擇。這些數字都是均勻分佈的，但大多數必須向下
捨入到最接近的可表示的 Python float。（2⁻¹⁰⁷⁴ 是最小為正的非正規化
float，等於 "math.ulp(0.0)"）

也參考:

  產生偽隨機浮點值  Allen B. Downey 的一篇論文描述了產生比通常由
  "random()" 產生的 float 更 fine-grained（細粒的）的方法。
