Python SciPy stats.jarque_bera用法及代碼示例

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

用法:

scipy.stats.jarque_bera(x, *, axis=None, nan_policy='propagate', keepdims=False)#

對樣本數據執行 Jarque-Bera 擬合優度檢驗。

Jarque-Bera 檢驗測試樣本數據是否具有與正態分布匹配的偏度和峰度。

請注意，此檢驗僅適用於足夠多的數據樣本 (>2000)，因為檢驗統計量漸近具有 2 個自由度的卡方分布。

參數 ：：

x： array_like

隨機變量的觀察。

axis： int 或無，默認：無

如果是 int，則計算統計量的輸入軸。輸入的每個axis-slice(例如行)的統計信息將出現在輸出的相應元素中。如果 None ，輸入將在計算統計數據之前被分解。

nan_policy： {‘propagate’, ‘omit’, ‘raise’}

定義如何處理輸入 NaN。

• propagate ：如果計算統計數據的軸切片(例如行)中存在NaN，則輸出的相應條目將為 NaN。

• omit : 計算時將省略NaNs。如果計算統計數據的軸切片中剩餘的數據不足，則輸出的相應條目將為 NaN。

• raise ：如果存在 NaN，則會引發 ValueError。

keepdims： 布爾值，默認值：假

如果將其設置為 True，則縮小的軸將作為尺寸為 1 的尺寸留在結果中。使用此選項，結果將針對輸入數組正確廣播。

返回 ：：

result： SignificanceResult

具有以下屬性的對象：

統計 浮點數

檢驗統計量。

p值 浮點數

假設檢驗的 p 值。

注意：

從 SciPy 1.9 開始，np.matrix 輸入(不建議用於新代碼)在執行計算之前轉換為 np.ndarray。在這種情況下，輸出將是標量或適當形狀的 np.ndarray 而不是 2D np.matrix 。同樣，雖然屏蔽數組的屏蔽元素被忽略，但輸出將是標量或 np.ndarray 而不是帶有 mask=False 的屏蔽數組。

參考：

[1]

Jarque, C. 和 Bera, A. (1980) “回歸殘差的正態性、同方差性和序列獨立性的有效檢驗”，6 Econometric Letters 255-259。

[2]

夏皮羅，S.S. 和威爾克，M.B. (1965)。正態性方差檢驗分析(完整樣本)。生物計量學，52(3/4)，591-611。

[3]

B. Phipson 和 G. K. Smyth。 “排列 P 值不應該為零：隨機抽取排列時計算精確的 P 值。”遺傳學和分子生物學中的統計應用 9.1 (2010)。

[4]

帕納吉奧塔科斯，D.B. (2008)。 p 值在生物醫學研究中的值。開放心血管醫學雜誌，2, 97。

例子：

假設我們希望從測量中推斷醫學研究中成年男性的體重是否不服從正態分布 [2]。重量(磅)記錄在下麵的數組x中。

>>> import numpy as np
>>> x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])

Jarque-Bera 測試首先根據樣本偏度和峰度計算統計量。

>>> from scipy import stats
>>> res = stats.jarque_bera(x)
>>> res.statistic
6.982848237344646

由於正態分布具有零偏度和零(“excess” 或 “Fisher”)峰度，因此對於從正態分布中抽取的樣本，此統計量的值往往較低。

該檢驗是通過將統計量的觀測值與零分布進行比較來執行的：零分布是在權重從正態分布中得出的零假設下得出的統計值的分布。對於Jarque-Bera 檢驗，非常大的樣本的零分布是具有兩個自由度的卡方分布。

>>> import matplotlib.pyplot as plt
>>> dist = stats.chi2(df=2)
>>> jb_val = np.linspace(0, 11, 100)
>>> pdf = dist.pdf(jb_val)
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> def jb_plot(ax):  # we'll re-use this
...     ax.plot(jb_val, pdf)
...     ax.set_title("Jarque-Bera Null Distribution")
...     ax.set_xlabel("statistic")
...     ax.set_ylabel("probability density")
>>> jb_plot(ax)
>>> plt.show()

比較通過 p 值進行量化：零分布中大於或等於統計觀測值的值的比例。

>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> jb_plot(ax)
>>> pvalue = dist.sf(res.statistic)
>>> annotation = (f'p-value={pvalue:.6f}\n(shaded area)')
>>> props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
>>> _ = ax.annotate(annotation, (7.5, 0.01), (8, 0.05), arrowprops=props)
>>> i = jb_val >= res.statistic  # indices of more extreme statistic values
>>> ax.fill_between(jb_val[i], y1=0, y2=pdf[i])
>>> ax.set_xlim(0, 11)
>>> ax.set_ylim(0, 0.3)
>>> plt.show()

>>> res.pvalue
0.03045746622458189

如果 p 值為 “small” - 也就是說，如果從正態分布總體中采樣數據產生統計數據的極值的概率較低 - 這可以作為反對零假設的證據另一種選擇：權重不是從正態分布中得出的。注意：

• 反之則不成立；也就是說，檢驗不用於為原假設提供證據。

• 將被視為 “small” 的值的閾值是在分析數據之前應做出的選擇 [3]，同時考慮誤報(錯誤地拒絕原假設)和漏報(未能拒絕假設)的風險。錯誤的原假設)。

請注意，卡方分布提供了零分布的漸近近似；它僅對於具有許多觀測值的樣本是準確的。對於像我們這樣的小樣本， scipy.stats.monte_carlo_test 可以提供更準確的(盡管是隨機的)精確 p 值的近似值。

>>> def statistic(x, axis):
...     # underlying calculation of the Jarque Bera statistic
...     s = stats.skew(x, axis=axis)
...     k = stats.kurtosis(x, axis=axis)
...     return x.shape[axis]/6 * (s**2 + k**2/4)
>>> res = stats.monte_carlo_test(x, stats.norm.rvs, statistic,
...                              alternative='greater')
>>> fig, ax = plt.subplots(figsize=(8, 5))
>>> jb_plot(ax)
>>> ax.hist(res.null_distribution, np.linspace(0, 10, 50),
...         density=True)
>>> ax.legend(['aymptotic approximation (many observations)',
...            'Monte Carlo approximation (11 observations)'])
>>> plt.show()

>>> res.pvalue
0.0097  # may vary

此外，盡管具有隨機性，以這種方式計算的 p 值可用於精確控製原假設的錯誤拒絕率 [4]。