Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

不完全β函數反函數

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

namespace boost{ namespace math{

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

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

template <class T1, class T2, class T3, class T4>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py);

template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy&);

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

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

template <class T1, class T2, class T3, class T4>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py);

template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy&);

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

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

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

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

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

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

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

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

}} // namespaces
說明

有6個不完全β函數反函數 允許你反解出不完全β函數中的三個參數中的任意一個, 從不完全β函數beta(p)的值開始或從它的補集(q)開始。

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

[Tip] 提示

當人們談及不完全γ函數的反函數的時候,它們往往是討論關於參數x的反函數。這在這個庫中實現為 ibeta_inv 和 ibeta_inv 函數, 也是在這裡出現的所有的反函數中最高效的。

關於參數a 在統計應用中有一定的應用,但是不得不使用相當蠻力的方法來計算,也因此速度要慢幾倍。這些函數在這個庫中實現為 ibetac_inva 和 ibetac_invb 函數.

當調用參數的類型T1...TN是不同類型的時候,函數的返回值的類型使用函數返回值推導法則 來確定.

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

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

template <class T1, class T2, class T3, class T4>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py);

template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy&);

返回滿足: q = ibetac(a, b, x)的x值; 當提供了參數參數py並且py不為NULL時,設置*py = 1 - x .注意:在這個函數的內部, 函數計算x和1-x中的無論哪一個有更小的值,因此傳遞給*py的值不會有消去誤差。這意味著即使函數返回1,在*py中存儲的值也是非0的,即使非常小。

要求: a,b > 0 and 0 <= p <= 1.

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

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

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

template <class T1, class T2, class T3, class T4>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py);

template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy&);

返回滿足: q = ibetac(a, b, x)的x值; 當提供了參數參數py並且py不為NULL時,設置*py = 1 - x .注意:在這個函數的內部, 函數計算x和1-x中的無論哪一個更小的值,因此傳遞給*py的值不會有消去誤差。這意味著即使函數返回1,在*py中存儲的值也是非0的,即使非常小。

要求: a,b > 00 <= q <= 1.

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

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

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

返回滿足: p = ibeta(a, b, x)的a值;

要求: b > 0, 0 < x < 10 <= p <= 1.

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

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

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

返回滿足: q = ibetac(a, b, x)的a值;

要求: b > 0, 0 < x < 10 <= q <= 1.

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

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

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

返回滿足: p = ibeta(a, b, x)的b值;

要求: a > 0, 0 < x < 10 <= p <= 1.

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

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

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

返回滿足: q = ibetac(a, b, x)的滿足;

要求: a > 0, 0 < x < 10 <= q <= 1.

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

精確性

這些函數的精度非常接近於規則的前向不完全β函數. 然而,注意在這些函數的一部分定義域中,函數對於輸入的參數的改變非常敏感,尤其是當參數p(或者補數q)非常接近於0或1 .

測試

有兩類測試:

函數 ibeta_inv 和函數 ibetac_inv的實現

這兩個函數共享一個通用的實現.

首先計算x的初始近似值,然後使用Halley 迭代清除結果的最後幾個bit. 迭代限制被設為T類型中的bit數目中的12個,通過實驗證明,這就中以使得這些反函數的精度至少與不完全β函數的精度一樣。在極端的情況下可能需要多達5次迭代,雖然通常只需要1到2次。更進一步,迭代次數隨著a 和 b 的增長而下降(這通常構成更重要的用例

迭代的初始猜測使用下面的方法來獲得:

首先回憶一下:

我們可能希望從p或q開始,並計算x或y。此外,如果結果更有利於解決問題,那麼在任何地方我們都可以交換a和b,p和q以及x與y。

對於 a+b >= 5 初始猜測使用:

Asymptotic Inversion of the Incomplete Beta Function, by N. M. Temme. Journal of Computational and Applied Mathematics 41 (1992) 145-157.

中描述的方法。

幾乎對稱的情況(在這篇論文的第2部分)

並且涉及到首先求解誤差函數。當a = 5時這種方法的精度至少為2位十進制數字,當a =105時,這種方法的精度至少為8位十進制數字.

一般的誤差函數情形 (這篇論文的第三部分)

在這裡,又一次以誤差函數的形式來表達不完全β函數. 當a+b=5時,這種方法的精度至少為2位十進制數字,當a+ b = 105時,這種方法的精度至少為11位數字. 然而,當結果預期很小時,以及當a + b很小時,那麼它的精度變得越來越小,在這種情況下,當 p1/a < 0.0025時,最好使用下面的方法進行初始評估:

最後,對於所有的a+b > 5的情況,使用在這篇論文的第四部分中的方法。通過不完全γ函數的形式來表達不完全β函數的反函數, 因此與其它的情況相比 ,這種方法的開銷明顯更大一些。然而,當a=5時, 這種方法的精度到少為2個十進制數字,當a= 105時,精確度上升到至少10個十進制數字。這種方法局限於a>b,因此當不是這種情況時,我們需要交換a 和b,p和q以及x和y. 此外,當p接近於1,這種方法是不精確的,我們應當將y而不是x作為輸出。因此,當q很小時 (q1/p < 10-3),我們使用:

這比完整的方法更高效,並且是對q的更精確的估計。

當a和b都很小時,在文獻資料中缺少相關的如何處理這種問題的信息。我非常感謝Prof Nico Temme提供了大量耐心來提供和解釋下面的信息。下面所有的紕漏都是我自己造成的,而不是Prof Nico Temme造成的。

當a和b都小於1時,在不完全β函數中在 xs = (1 - a) / (2 - a - b)處有一個拐折點。因此,如果 p > Ix(a,b),我們交換 a 和 b ,p 和 q , 以及 x 和y ,因此,現在我們將會一直在拐折點xs下面找出一個點x,且在一條凸曲線上。使用下面的方法對x進行初始估值:

這可能在x的真值的下面:牛頓迭代 將會順利地關於x收斂而不會有因為overshooting產生的問題.

當a和b都大於1時, 但a+b太小,以至於不能使用上面提到的任何方法,我們處理的過程如下。. 觀察到在點xs = (1 - a) / (2 - a - b)處,不完全β函數具有一個拐折點 , 因此,現在我們將會一直在拐折點xs下面找出一個點x,且在一條凸曲線上。使用下面的方法對x進行初始估值:

這還可以改進為:

當 b 和x都很小(這裡使用 b < a 且 x < 0.2). 這實際上仰仗了x。這就把我們放在了使得牛頓迭代以便的x的錯誤一邊。然而,使用高階導數和Halley 迭代可以使得事情在可控的範圍之內。

最一要考慮的情況是當a和b中的一個小於或等於1,而另外一個大於1.在這時,如果b<1,我們交換 a 和 b ,p 和 q , 以及 x 和y。現在,不完全β函數的曲線是凸出的且在範圍[0,1]內沒有拐折點。對於小的p,x可以使用下面的方法進行估計:

這低估了x, 我們放在了使得牛頓迭代以便的x的正確的一邊。然而,當p很大時,這可以嚴重低估了x。當我們真的想要查找y時,這就尤其是一個問題,在這種情況下,這種方法可能是一個任意數目次數的大小的輸出,導致非常差的收斂.

當把不完全β函數看作是四分之一個歪曲的圓周時,事情可以得到改進,並且可以從下面的形式中來估算y值:

這並不保證我們可以處在x的單調收斂的正確一邊,但是,它確實使得我們足夠接近從面Halley 迭代在真值上可以快速收斂.

實現關於參數a和b的反函數

這四個函數共享同一個通用實現.

首先對a或b進行初始近似:當可能的時候針對於負的二項分佈使用 Cornish-Fisher 展開來使得結果在1附近。然而,當a或b都很小時, Cornish Fisher 展開就不可用 了,在這種情況下,選擇初始近似使得 Ix(a, b) 接近於範圍 [0,1]的中間.

初始估測值就作為一個值的根查找算法的起始值。一旦它被括號括起來,這個函數在根上收斂很快,但是將根括起來可能需要幾次迭代。一個對a和b的更好的初始近似將證明這些函數非常本質:現在,查找根大約需要調用10-20個不完全β函數.


PrevUpHomeNext