Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

不完全β函數

概要

#include <boost/math/special_functions/gamma.hpp>

namespace boost{ namespace math{

template <class T1, class T2, class T3>
calculated-result-type ibeta(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta(T1 a, T2 b, T3 x, const Policy&);

template <class T1, class T2, class T3>
calculated-result-type ibetac(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac(T1 a, T2 b, T3 x, const Policy&);

template <class T1, class T2, class T3>
calculated-result-type beta(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type beta(T1 a, T2 b, T3 x, const Policy&);

template <class T1, class T2, class T3>
calculated-result-type betac(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type betac(T1 a, T2 b, T3 x, const Policy&);

}} // namespaces
說明

有4個不完全β函數 : 兩個標準版本( normalised versions) (也稱作規則 β函數 ) 返回值的區間為 [0, 1],以及兩個非標準版本( non-normalised )且返回值的區間為[0, beta(a, b)]. 對於統計應用感興趣的用戶應當使用標準版本 (或規則版本 ) (ibeta 和 ibetac).

所有的這些函數都要求: a > 0, b > 00 <= x <= 1.

當T1,T2和T3是不同類型的時候,函數的返回值的類型使用函數返回值推導法則 來確定。.

最後一個策略 參數是可選的並且可以用來控制函數的行為: 如何處理錯誤, 使用哪種層次的精度等等. 參見策略文檔瞭解更多信息。

template <class T1, class T2, class T3>
calculated-result-type ibeta(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta(T1 a, T2 b, T3 x, const Policy&);

返回 a, b 和 x的標準不完全β函數。

template <class T1, class T2, class T3>
calculated-result-type ibetac(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac(T1 a, T2 b, T3 x, const Policy&);

返回a,b和x的完整的標準完全和標準不完全β函數:

template <class T1, class T2, class T3>
calculated-result-type beta(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type beta(T1 a, T2 b, T3 x, const Policy&);

返回a,b和x的 (非標準) 不完全β函數:

template <class T1, class T2, class T3>
calculated-result-type betac(T1 a, T2 b, T3 x);

template <class T1, class T2, class T3, class Policy>
calculated-result-type betac(T1 a, T2 b, T3 x, const Policy&);

返回a,b和x的 (不標準) 不完全β函數的補集:

精確度

下面的表顯示了輸入參數a, b 和 x的不同定義域的峰值誤差, 以及與GSL-1.9 庫和Cephes 庫的比較. 注意:在一些平台上使用寬浮點類型來表示窄浮點類型的結果具有有效的零誤差.

注意:對於 80 和 128-bit long double 的精度明顯要比double 高:這是因為更寬的指數範圍允許更測試更極端的測試情況.例如,在double精度預期為0的結果在指數範圍更寬的long double類型中是有限的,但是非常小。

表3.函數 ibeta(a,b,x) 的出錯率

有效數字位數

平台和編譯器

0 < a,b < 10

0 < x < 1

0 < a,b < 100

0 < x < 1

1x10-5 < a,b < 1x105

0 < x < 1

53

Win32, Visual C++ 8

峰值=42.3 均值=2.9

(GSL 峰值=682 均值=32.5)

(Cephes 峰值=42.7 均值=7.0)

峰值=108 均值=16.6

(GSL 峰值=690 均值=151)

(Cephes 峰值=1545 均值=218)

峰值=4x103 均值=203

(GSL 峰值~3x105 均值~2x104)

(Cephes 峰值~5x105 均值~2x104)

64

Redhat Linux IA32, gcc-3.4.4

峰值=21.9 均值=3.1

峰值=270.7 均值=26.8

峰值~5x104 均值=3x103

64

Redhat Linux IA64, gcc-3.4.4

峰值=15.4 均值=3.0

峰值=112.9 均值=14.3

峰值~5x104 均值=3x103

113

HPUX IA64, aCC A.06.06

峰值=20.9 均值=2.6

峰值=88.1 均值=14.3

峰值~2x104 均值=1x103


表4.函數 ibetac(a,b,x) 的出錯率

有效數字位數

平台和編譯器

0 < a,b < 10

0 < x < 1

0 < a,b < 100

0 < x < 1

1x10-5 < a,b < 1x105

0 < x < 1

53

Win32, Visual C++ 8

峰值=13.9 均值=2.0

峰值=56.2 均值=14

峰值=3x103 均值=159

64

Redhat Linux IA32, gcc-3.4.4

峰值=21.1 均值=3.6

峰值=221.7 均值=25.8

峰值~9x104 均值=3x103

64

Redhat Linux IA64, gcc-3.4.4

峰值=10.6 均值=2.2

峰值=73.9 均值=11.9

峰值~9x104 均值=3x103

113

HPUX IA64, aCC A.06.06

峰值=9.9 均值=2.6

峰值=117.7 均值=15.1

峰值~3x104 均值=1x103


表5.函數 beta(a, b, x)的出錯率

有效數字位數

平台和編譯器

0 < a,b < 10

0 < x < 1

0 < a,b < 100

0 < x < 1

1x10-5 < a,b < 1x105

0 < x < 1

53

Win32, Visual C++ 8

峰值=39 均值=2.9

峰值=91 均值=12.7

峰值=635 均值=25

64

Redhat Linux IA32, gcc-3.4.4

峰值=26 均值=3.6

峰值=180.7 均值=30.1

峰值~7x104 均值=3x103

64

Redhat Linux IA64, gcc-3.4.4

峰值=13 均值=2.4

峰值=67.1 均值=13.4

峰值~7x104 均值=3x103

113

HPUX IA64, aCC A.06.06

峰值=27.3 均值=3.6

峰值=49.8 均值=9.1

峰值~6x104 均值=3x103


表6.函數 betac(a,b,x)的出錯率

有效數字位數

平台和編譯器

0 < a,b < 10

0 < x < 1

0 < a,b < 100

0 < x < 1

1x10-5 < a,b < 1x105

0 < x < 1

53

Win32, Visual C++ 8

峰值=12.0 均值=2.4

峰值=91 均值=15

峰值=4x103 均值=113

64

Redhat Linux IA32, gcc-3.4.4

峰值=19.8 均值=3.8

峰值=295.1 均值=33.9

峰值~1x105 均值=5x103

64

Redhat Linux IA64, gcc-3.4.4

峰值=11.2 均值=2.4

峰值=63.5 均值=13.6

峰值~1x105 均值=5x103

113

HPUX IA64, aCC A.06.06

峰值=15.6 均值=3.5

峰值=39.8 均值=8.9

峰值~9x104 均值=5x103


測試

有兩種類別的測試: 抽樣測試(spot tests)比較數據取自於Mathworld's 在線函數求值 和這個庫: 它們對這個庫的實現提供了一個基本的 "合理性測試(sanity check)" , 對第一種實現都有一個抽取測試(spot-test ) (參見下面的實現說明).

精度測試使用精度非常高的測試數據 (NTL RR 類 設置為 1000-bit 精度), 使用"教書(textbook)"連分數表示( continued fraction representation )(參見下面的實現討論中的第一個連分數).注意:連分數並沒有用在這個庫的實現中, 因此我們的測試數據完全獨立於這些代碼。

實現

這個庫的實現緊密地基於"Algorithm 708; Significant digit computation of the incomplete beta function ratios", DiDonato and Morris, ACM, 1992.

所有的這四個函數共享一個通用的實現: 傳遞 x 和 y, 且返回 p 或 q 。這些變量通過下面的等式等起來:

因此,如果結果處於一個更好的情況,我們可以在任何位置交換 a 和 b, x 和 y 以及 p 和 q . 通常都會進行這些交換,這使得我們永遠都是計算一個小於0.9的值,當需要的時候,這可能從1中減去而不會有過度的消去誤差(cancellation error)

下面的連分數表示可以在很多教科書中找到,但是卻沒有用在這個庫的實現中-與其它的實現相比,它不但更慢而且精度更低 - 然而它被用來生成測試數據:

下面的連分數出自於Didonato and Morris, 當a 和 b 都大於1的時候,下面的等式用在這個庫的實現中:

對於較小的b和x可以使用級數表示:

當 b << a時,從0到1的轉換的出現非常接近於 x = 1 並且這就需要一些額外的處理來接管這個計算方法,在那種情況下使用下面的級數表示:

其中 Q(a,x)是一個不完全γ函數. 注意:這種方法依賴於保持一個已計算的 pn 值的表, 這就限制了這種方法的精度,依賴於所使用的表的大小。

ab 都是較小的整數時, 那麼我們可以將不完全的β函數和二項分佈聯繫起來並使用下面的有限和:

最後, 我們可以迴避迴避的區域, 或者轉移到一個更高效的計算方式,通過使用下面的倍角公式:

各種計算方法所使用的a,b和x的定義域與Didonato and Morris TOMS 708 paper.所描述的是等價的。


PrevUpHomeNext