Python SciPy integrate.qmc_quad用法及代碼示例

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

用法:

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

使用 Quasi-Monte 卡洛求積計算 N-dimensions 中的積分。

參數 ：：

func： 可調用的

被積函數。必須接受單個參數x，一個數組，指定計算 scalar-valued 被積函數的點，並返回被積函數的值。為了提高效率，函數應該被向量化以接受形狀數組(d, n_points)，其中d是變量的數量(即函數域的維數)並且n_points是正交點數，返回形狀數組(n_points,)，每個正交點的被積函數。

a, b： 類數組

一維數組分別指定每個 d 變量的積分下限和上限。

n_estimates, n_points： 整數，可選

n_estimates(默認值：8)統計上獨立的 QMC 樣本，每個樣本n_points(默認：1024)點，將由以下方式生成qrng。被積函數的總點數函數將被評估為n_points * n_estimates。詳細信息請參見注釋。

qrng： QMCEngine ，可選

從中采樣 QMC 點的 QMCEngine 實例。 QMCEngine 必須初始化為多個維度d與變量的數量相對應x1, ..., xd傳遞給函數。提供的 QMCEngine 用於生成第一個積分估計。如果n_estimates大於 1，則從第一個生成額外的 QMCEngine(如果有選項，則啟用加擾。)如果未提供 QMCEngine，則使用默認值scipy.stats.qmc.Halton將使用根據長度確定的維數進行初始化a.

log： 布爾值，默認：False

當設置為 True 時，func 返回被積函數的對數，結果對象包含積分的對數。

返回 ：：

result： 對象

具有屬性的結果對象：

不可缺少的 浮點數

積分的估計。

standard_error :

誤差估計。解釋見注釋。

注意：

每個被積函數的值n_pointsQMC 樣本的點用於生成積分的估計。該估計值是從積分的可能估計值總體中得出的，我們獲得的值取決於評估積分的特定點。我們執行這個過程n_estimates次，每次評估不同擾亂 QMC 點處的被積函數，有效地繪製獨立同分布。來自整體估計總體的隨機樣本。樣本均值\(m\)這些積分估計值是積分真實值的無偏估計量，以及平均值的標準誤差\(s\)這些估計值可用於使用 t 分布生成置信區間n_estimates - 1自由度。也許與直覺相反，增加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