Ctrl+K
scipy.integrate.

qmc_quad#

scipy.integrate.qmc_quad(func, a, b, *, n_estimates=8, n_points=1024, qrng=None, log=False)[原始碼]#

使用準蒙地卡羅積分法計算 N 維積分。

參數:
funccallable

被積分函數。必須接受單一引數 x，一個指定要評估純量值被積分函數的點的陣列，並傳回被積分函數的值。為了提高效率，該函數應該向量化以接受形狀為 (d, n_points) 的陣列，其中 d 是變數的數量（即函數域的維度），而 n_points 是積分點的數量，並傳回形狀為 (n_points,) 的陣列，即每個積分點的被積分函數。

a, barray-like

一維陣列，分別指定每個 d 變數的積分下限和上限。

n_estimates, n_pointsint, optional

n_estimates (預設值：8) 統計上獨立的 QMC 樣本，每個樣本有 n_points (預設值：1024) 個點，將由 qrng 產生。被積分函數 func 將被評估的總點數為 n_points * n_estimates。詳情請參閱「注意事項」。

qrngQMCEngine, optional

QMCEngine 的實例，用於从中採樣 QMC 點。QMCEngine 必須初始化為維度數量 d，該數量與傳遞給 func 的變數 x1, ..., xd 的數量相對應。提供的 QMCEngine 用於產生第一個積分估計值。如果 n_estimates 大於 1，則從第一個 QMCEngine 衍生出額外的 QMCEngine（如果選項允許，則啟用擾亂）。如果未提供 QMCEngine，則預設的 scipy.stats.qmc.Halton 將使用從 a 的長度確定的維度數量進行初始化。

logboolean, default: False

當設定為 True 時，func 傳回被積分函數的對數，且結果物件包含積分的對數。

傳回值:
resultobject

具有屬性的結果物件

integralfloat

積分的估計值。

standard_error

誤差估計值。有關解釋，請參閱「注意事項」。

注意事項

在 QMC 樣本的每個 n_points 點處的被積分函數值用於產生積分的估計值。此估計值取自積分的可能估計值的母體，我們獲得的值取決於積分被評估的特定點。我們執行此過程 n_estimates 次，每次在不同的擾亂 QMC 點評估被積分函數，有效地從積分估計值的母體中抽取 i.i.d. 隨機樣本。這些積分估計值的樣本平均值 \(m\) 是積分真實值的無偏估計量，而這些估計值的平均值 \(s\) 的標準誤差可用於使用自由度為 n_estimates - 1 的 t 分佈產生信賴區間。也許違反直覺的是，增加 n_points 同時保持函數評估點的總數 n_points * n_estimates 固定往往會減少實際誤差，而增加 n_estimates 往往會減少誤差估計值。

範例

QMC 積分法特別適用於計算較高維度的積分。一個範例被積分函數是多變數常態分佈的機率密度函數。

>>> import numpy as np
>>> from scipy import stats
>>> dim = 8
>>> mean = np.zeros(dim)
>>> cov = np.eye(dim)
>>> def func(x):
...     # `multivariate_normal` expects the _last_ axis to correspond with
...     # the dimensionality of the space, so `x` must be transposed
...     return stats.multivariate_normal.pdf(x.T, mean, cov)

若要計算單位超立方體上的積分

>>> from scipy.integrate import qmc_quad
>>> a = np.zeros(dim)
>>> b = np.ones(dim)
>>> rng = np.random.default_rng()
>>> qrng = stats.qmc.Halton(d=dim, seed=rng)
>>> n_estimates = 8
>>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
>>> res.integral, res.standard_error
(0.00018429555666024108, 1.0389431116001344e-07)

雙邊 99% 積分信賴區間可以估計為

>>> t = stats.t(df=n_estimates-1, loc=res.integral,
...             scale=res.standard_error)
>>> t.interval(0.99)
(0.0001839319802536469, 0.00018465913306683527)

事實上，scipy.stats.multivariate_normal 報告的值在此範圍內。

>>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
0.00018430867675187443
在本頁中