Python SciPy sampling.DiscreteGuideTable用法及代碼示例

本文簡要介紹 python 語言中 scipy.stats.sampling.DiscreteGuideTable 的用法。

用法:

class  scipy.stats.sampling.DiscreteGuideTable(dist, *, domain=None, guide_factor=1, random_state=None)#

離散導表法。

Discrete Guide Table 方法從任意但有限的概率向量中采樣。它使用大小為\(N\) 的概率向量或具有有限支持的概率質量函數從分布中生成隨機數。 Discrete Guide Table 的設置速度非常慢(與矢量長度成線性關係)，但采樣速度非常快。

參數 ：：

dist： 數組 或對象，可選

分布的概率向量 (PV)。如果 PV 不可用，則需要一個具有 pmf 方法的類的實例。 PMF 的簽名預計為：def pmf(self, k: int) -> float。即它應該接受一個 Python 整數並返回一個 Python 浮點數。

domain： 整數，可選

支持 PMF。如果概率向量 (pv) 不可用，則必須給出有限域。即 PMF 必須有一個有限的支持。默認為 None 。當 None ：

• 如果一個support方法由分發對象提供距離，用於設置分布的域。

• 否則，假定支持為 (0, len(pv)) 。當此參數與概率向量一起傳遞時，domain[0] 用於將分布從 (0, len(pv)) 重新定位到 (domain[0], domain[0]+len(pv))，而 domain[1] 將被忽略。有關更詳細的說明，請參閱注釋和教程。

guide_factor: int, optional：

相對於 PV 長度的導向台尺寸。更大的引導表會導致更快的生成時間，但需要更昂貴的設置。不推薦大於 3 的尺寸。如果相對大小設置為 0，則使用順序搜索。默認值為 1。

random_state： {無，整數， numpy.random.Generator ，

numpy.random.RandomState}, optional

NumPy 隨機數生成器或底層NumPy 隨機數生成器的種子，用於生成統一隨機數流。如果random_state是無(或np.random)， 這numpy.random.RandomState使用單例。如果random_state是一個 int，一個新的RandomState使用實例，播種random_state.如果random_state已經是一個Generator或者RandomState實例然後使用該實例。

注意：

當有限概率向量可用或分布的 PMF 可用時，此方法有效。如果隻有 PMF 可用，則還必須給出 PMF 的有限支持(域)。建議首先通過評估支持中每個點的 PMF 來獲得概率向量，然後再使用它。

DGT 樣本來自任意但有限的概率向量。隨機數是通過反演方法生成的，即

1. 生成一個隨機數 U ~ U(0,1)。

2. 找到最小整數 I 使得 F(I) = P(X<=I) >= U。

步驟(2)是關鍵步驟。使用順序搜索需要 O(E(X)) 比較，其中 E(X) 是分布的期望值。然而，索引搜索使用引導表跳轉到靠近 I 的一些 I' <= I 以在恒定時間內找到 X。實際上，當引導表與概率向量具有相同大小(這是默認值)時，預期的比較次數減少到 2。對於較大的引導表，此數字會變小(但始終大於 1)，對於較小的表，它會變大。對於表大小為 1 的極限情況，該算法隻是進行順序搜索。

另一方麵，引導表的設置時間為 O(N)，其中 N 表示概率向量的長度(對於大小 1，不需要預處理)。此外，對於非常大的引導表，內存效應甚至可能會降低算法的速度。所以我們不建議使用比給定概率向量大三倍以上的引導表。如果隻需要生成幾個隨機數，(很多)較小的表大小會更好。引導表相對於給定概率向量長度的大小可以通過guide_factor參數設置。

如果給定一個概率向量，它必須是一個非負浮點數的一維數組，沒有任何 inf 或 nan 值。此外，必須至少有一個非零條目，否則會引發異常。

默認情況下，概率向量從 0 開始索引。但是，這可以通過傳遞 domain 參數來更改。當 domain 與 PV 一起給出時，它具有將分布從 (0, len(pv)) 重新定位到 (domain[0], domain[0] + len(pv)) 的效果。 domain[1] 在這種情況下被忽略。

參考：

[1]

UNU.RAN 參考手冊，第 5.8.4 節，“DGT -(離散)引導表方法(索引搜索)”https://statmath.wu.ac.at/unuran/doc/unuran.html#DGT

[2]

HC陳和 Y. Asau (1974)。關於從經驗分布中生成隨機變量，AIIE Trans。 6，第 163-166 頁。

例子：

>>> from scipy.stats.sampling import DiscreteGuideTable
>>> import numpy as np

要使用概率向量創建隨機數生成器，請使用：

>>> pv = [0.1, 0.3, 0.6]
>>> urng = np.random.default_rng()
>>> rng = DiscreteGuideTable(pv, random_state=urng)

RNG 已設置完畢。現在，我們現在可以使用 rvs 方法從分布中生成樣本：

>>> rvs = rng.rvs(size=1000)

為了驗證隨機變量是否遵循給定分布，我們可以使用卡方檢驗(作為goodness-of-fit 的度量)：

>>> from scipy.stats import chisquare
>>> _, freqs = np.unique(rvs, return_counts=True)
>>> freqs = freqs / np.sum(freqs)
>>> freqs
array([0.092, 0.355, 0.553])
>>> chisquare(freqs, pv).pvalue
0.9987382966178464

由於 p 值非常高，我們無法拒絕觀察到的頻率與預期頻率相同的零假設。因此，我們可以安全地假設變量是根據給定的分布生成的。請注意，這隻是給出了算法的正確性，而不是樣本的質量。

如果 PV 不可用，也可以傳遞具有 PMF 方法和有限域的類的實例。

>>> urng = np.random.default_rng()
>>> from scipy.stats import binom
>>> n, p = 10, 0.2
>>> dist = binom(n, p)
>>> rng = DiscreteGuideTable(dist, random_state=urng)

現在，我們可以使用 rvs 方法從分布中采樣，並測量樣本的 goodness-of-fit：

>>> rvs = rng.rvs(1000)
>>> _, freqs = np.unique(rvs, return_counts=True)
>>> freqs = freqs / np.sum(freqs)
>>> obs_freqs = np.zeros(11)  # some frequencies may be zero.
>>> obs_freqs[:freqs.size] = freqs
>>> pv = [dist.pmf(i) for i in range(0, 11)]
>>> pv = np.asarray(pv) / np.sum(pv)
>>> chisquare(obs_freqs, pv).pvalue
0.9999999999999989

為了檢查樣本是否來自正確的分布，我們可以可視化樣本的直方圖：

>>> import matplotlib.pyplot as plt
>>> rvs = rng.rvs(1000)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> x = np.arange(0, n+1)
>>> fx = dist.pmf(x)
>>> fx = fx / fx.sum()
>>> ax.plot(x, fx, 'bo', label='true distribution')
>>> ax.vlines(x, 0, fx, lw=2)
>>> ax.hist(rvs, bins=np.r_[x, n+1]-0.5, density=True, alpha=0.5,
...         color='r', label='samples')
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('PMF(x)')
>>> ax.set_title('Discrete Guide Table Samples')
>>> plt.legend()
>>> plt.show()

要設置引導表的大小，請使用 guide_factor 關鍵字參數。這設置了引導表相對於概率向量的大小

>>> rng = DiscreteGuideTable(pv, guide_factor=1, random_state=urng)

要使用 \(n=4\) 和 \(p=0.1\) 計算二項分布的 PPF：我們可以設置如下指導表：

>>> n, p = 4, 0.1
>>> dist = binom(n, p)
>>> rng = DiscreteGuideTable(dist, random_state=42)
>>> rng.ppf(0.5)
0.0