Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext
負二項分佈(Negative Binomial Distribution)

#include <boost/math/distributions/negative_binomial.hpp>

namespace boost{ namespace math{ 

template <class RealType = double, 
          class Policy   = policies::policy<> >
class negative_binomial_distribution;

typedef negative_binomial_distribution<> negative_binomial;

template <class RealType, class Policy>
class negative_binomial_distribution
{
public:
   typedef RealType value_type;
   typedef Policy   policy_type;
   // 基於 successes 和 success_fraction的構造函數:
   negative_binomial_distribution(RealType r, RealType p);
   
   // 參數訪問函數(Parameter accessors):
   RealType success_fraction() const;
   RealType successes() const;
   
   // 成功分數(sucess fraction)分界:
   static RealType find_lower_bound_on_p(
      RealType trials, 
      RealType successes,
      RealType probability); // alpha
   static RealType find_upper_bound_on_p(
      RealType trials, 
      RealType successes,
      RealType probability); // alpha
      
   // 估算試驗的 最小/最大次數:
   static RealType find_minimum_number_of_trials(
      RealType k,     // 試驗失敗次數.
      RealType p,     // 成功分數(Success fraction).
      RealType probability); // 概率起始值 alpha.
   static RealType find_maximum_number_of_trials(
      RealType k,     // 試驗失敗次數.
      RealType p,     // 成功分數(Success fraction).
      RealType probability); // 概率起始值 alpha.
};

}} // namespaces

類型 negative_binomial_distribution 表示一個負二項分佈(negative_binomial distribution):當伯努利試驗(Bernoulli trial)有兩個互斥的結果時使用這個分佈。這些試驗結果標記為「成功」和「失敗」。

對於成功分數(success fraction)為p的 k + r 次伯努利試驗,負二項分佈給出在最後一次試驗成功時有k次失敗以及r次成功的概率。負二項分佈假定對於所有的(k+r)次試驗成功分數(success fraction)p是固定不變的。

[Note] 注意

負二項分佈的隨機變量是試驗的次數,(成功試驗的次數是負二項分佈的固有屬性);而對於二項分佈,隨機變量是成功試驗的次數,對於一個固定的試驗次數。

負二項分佈的PDF函數為:

下面的圖像顯示了當成功分數(success fraction)p變化時,函數PDF圖像的變化:

另一方面,下面的圖像顯示了當試驗成功次數發生變化時,函數PDF圖像的變化:

相關分佈

負二項分佈這個名字被一些人留給成功參數(success parameter)r為整數的情況。整數版本的負二項分佈也被稱作帕斯卡分佈(Pascal distribution)

這個庫的實現始終使用實數進行運算(因為負二項分佈的實現使用實值(real-valued) 版本的不完全β函數)。 實值(real-valued)版本的負二項分佈也稱作Polya 分佈。

泊松分佈(Poisson distribution) 是 帕斯卡分佈(Pascal distribution)的一般化(generalization), 其中,成功參數(success parameter) r 是一個整數:為了獲得帕斯卡分佈(Pascal distribution),你必須為r提供一個整數值,並且從返回成功次數的函數中獲取整數值(向上取整 或 向下取整)。

對於數值較大的r (成功試驗), 負二項分佈(negative binomial distribution)收斂為( converges to)泊松分佈(Poisson distribution)。

幾何分佈(geometric distribution)是成功參數(successes parameter)r = 1的特殊情況,因此僅需要一個first和成功分數。 geometric(p) = negative_binomial(1, p).

泊松分佈(Poisson distribution)是成功參數較大的特殊情況:

poisson(λ) = lim r → ∞ negative_binomial(r, r / (λ + r)))

[Caution] 注意

負二項分佈是一個離散分佈:內部函數(例如cdf和pdf)被當作它們「好像是」連續函數一樣,但實際上,僅當將整數值提供給隨機變量的時候這些函數才返回有意義的值。

分位點函數將會缺省返回一個向外捨入( rounded outwards)的整數值。也就是說,下分位點(lower quantiles)(概率小於0.5)向下捨入;上分位點(upper quantiles)(概率大於0.5)向上捨入。這種行為確保如果返回一個X%分位點值,那麼至少目標的覆蓋範圍將會在中心區域顯示,不是要求的覆蓋範圍將會在尾部(tails)顯示。

這種行為可以改變,使得分位點函數可以進行不同的捨入,或者使用策略來返回一個實值(real-valued)。在你使用二項分佈的分位點函數之前,強烈推薦你閱讀理解分佈的分位點參考文檔 描述了如何為這些分佈改變捨入策略。

成員函數
構造
negative_binomial_distribution(RealType r, RealType p);

構造函數: r 是總的成功試驗次數,p 是單次試驗成功的概率.

要求:r > 00 <= p <= 1.

訪問函數(Accessors)
RealType success_fraction() const; // successes / trials (0 <= p <= 1)

返回構造這個負二項分佈的參數p

RealType successes() const; // required successes (r > 0)

返回構造這個負二項分佈的參數r

參數p的下邊界(Lower Bound on Parameter p)
static RealType find_lower_bound_on_p(
  RealType failures, 
  RealType successes,
  RealType probability) // (0 <= alpha <= 1), 0.05 等價於 95% 的把握.

返回成功分數(success fraction)的下邊界

failures

在第r次成功試驗之前的失敗試驗的總次數。

successes

目標的成功試驗總數。

alpha

成功分數(success fraction)的真值(true value)小於 返回值的最大可接受概率。

例如,如果你想在 n = k + r 次試驗中觀察到k 次失敗和r 次成功,成功分數(success fraction)的最好估計是r/n,但是如果你想要有95%的把握:真值(true value)會 大於某個值,pmin,那麼:

pmin = negative_binomial_distribution<RealType>::find_lower_bound_on_p(
                    failures, successes, 0.05);

參看負二項分佈置信區間(confidence interval)例子.

這個函數使用Clopper-Pearson 方法來計算成功分數(success fraction)的下邊界,然而很多的資料中提到這個方法給出"精確"的結果,但在實際中,它產生一個至少是要求的覆蓋範圍的區間,對於一些失敗試驗和成功試驗的組合,這種方法可能產生悲觀的估計(pessimistic estimates)。參考:

Yong Cai and K. Krishnamoorthy, A Simple Improved Inferential Method for Some Discrete Distributions. Computational statistics and data analysis, 2005, vol. 48, no3, 605-621.

參數p的上邊界(Upper Bound on Parameter p)
static RealType find_upper_bound_on_p(
   RealType trials, 
   RealType successes,
   RealType alpha); // (0 <= alpha <= 1), 0.05 等價於 95% 的把握.

返回成功分數(success fraction)的上邊界

trials

總共的試驗次數.

successes

成功試驗的次數.

alpha

成功分數的真值(true value)大於返回值的最大可接受概率。

例如,如果你想在 n 次試驗中觀察到k 次成功,成功分數(success fraction)的最好估計是k/n,但是如果你想要有95%的把握:真值(true value)會 小於某個值, pmax,那麼:

pmax = negative_binomial_distribution<RealType>::find_upper_bound_on_p(
                    r, k, 0.05);

參看負二項分佈置信區間(confidence interval)例子.

這個函數使用Clopper-Pearson 方法來計算成功分數(success fraction)的上邊界,然而很多的資料中提到這個方法給出"精確"的結果,但在實際中,它產生一個至少是要求的覆蓋範圍的區間,對於一些失敗試驗和成功試驗的組合,這種方法可能產生悲觀的估計(pessimistic estimates)。參考:

Yong Cai and K. Krishnamoorthy, A Simple Improved Inferential Method for Some Discrete Distributions. Computational statistics and data analysis, 2005, vol. 48, no3, 605-621.

估計至少發生特定次數的失敗試驗所需要的總試驗次數
static RealType find_minimum_number_of_trials(
   RealType k,     // 失敗試驗次數.
   RealType p,     // 成功分數(success fraction).
   RealType alpha); // 概率起始值 (0.05 等價於 95%).

這個函數用於估算為了達到一個特定的概率( 觀察到超過k次失敗試驗)所需要的試驗總數。

k

最大的失敗試驗總數.

p

每次試驗的成功概率.

alpha

觀察到的失敗試驗總數為k次或少於k次的最大可接受風險值.

例如:
negative_binomial_distribution<RealType>::find_minimum_number_of_trials(10, 0.5, 0.05);

返回具有95%把握的以1/2頻率發生的10次失敗試驗最小的試驗次數。

可運行的例子.

這個函數使用負二項分佈的數值反演(numeric inversion)來獲取結果:結果的另一種解釋是,查找試驗的次數(成功+失敗)使得發生k次或更少的失敗試驗的概率為alpha

估算試驗的總數來確保失敗試驗的最大或更少的次數(Estimating Number of Trials to Ensure a Maximum Number of Failures or Less)
static RealType find_maximum_number_of_trials(
   RealType k,     // 失敗試驗數.
   RealType p,     // 成功分數(success fraction).
   RealType alpha); // 概率起始值 (0.05 等價於 95%).

這個函數用於估算達到k次或更少的失敗試驗次數的概率所需要的試驗最大次數。

k

最大的失敗試驗次數.

p

每次試驗成功的概率.

alpha

發生多於k次失敗試驗的最大可接受風險值。

 

例如:

negative_binomial_distribution<RealType>::find_maximum_number_of_trials(0, 1.0-1.0/1000000, 0.05);

返回最大的試驗次數,並且仍然有95%的可能性以百萬分之一的頻率不發生任何失敗試驗。

這個函數使用負二項分佈的數值反演(numeric inversion)來獲取結果:結果的另一種解釋是,查找試驗的次數(成功+失敗)使得發生超過k失敗試驗的概率為alpha

非成員訪問函數(Non-member Accessors)

支持所有的分佈都通用的 常見的非成員訪問函數累積分佈函數(Cumulative Distribution Function)概率密度函數(Probability Density Function)分位點(Quantile)故障率函數(Hazard Function)累積危險函數(Cumulative Hazard Function)均值(mean)中位數(median)眾數(mode)方差(variance)標準差(standard deviation)偏斜(skewness)峰態(kurtosis)峰態超越(kurtosis_excess)值域(range) 以及 支持(support)

值得花一時間來定義在負二項分佈中的這些訪問函數的實際含義:

表2.非成員訪問函數的含義.

函數

含義

概率密度函數(Probability Density Function)

在成功分數(success fraction)為p的k+r次試驗中產生剛好k次失敗試驗 且最後一次試驗成功的概率。例如:

pdf(negative_binomial(r, p), k)

累積分佈函數(Cumulative Distribution Function)

在成功分數(success fraction)為p的k+r次試驗中產生k次或更少次數的失敗試驗且最後一次試驗成功的概率。例如:

cdf(negative_binomial(r, p), k)

累積分佈函數的補集(Complement of the Cumulative Distribution Function)

在成功分數(success fraction)為p的k+r次試驗中產生超過k次失敗試驗 且最後一次試驗成功的概率。例如:

cdf(complement(negative_binomial(r, p), k))

分位點(Quantile)

在成功分數(success fraction)為p,概率為P的k+r次試驗中可以產生的最大失敗試驗次數k。注意:返回值是一個實數而不是一個整數。依賴於使用情況,你可能想要向上(ceiling)捨入或是向下(floor)捨入的結果值。例如:

quantile(negative_binomial(r, p), P)

概率補集的分位點(Quantile from the complement of the probability)

在成功分數(success fraction)為p,概率為P的k+r次試驗中可以產生的最小失敗試驗次數k。注意:返回值是一個實數而不是一個整數。依賴於使用情況,你可能想要向上(ceiling)捨入或是向下(floor)捨入的結果值。例如:

quantile(complement(negative_binomial(r, p), P))


精確度

負二項分佈使用不完全β函數ibetaibetac來實現,請參看這些函數瞭解關於精度的信息。

實現

在下面的表中,p 是單次試驗的成功概率(成功分數),r 是成功試驗的次數,k是失敗試驗的次數, p是概率,且 q = 1-p

函數

實現註解

pdf

pdf = exp(lgamma(r + k) - lgamma(r) - lgamma(k+1)) * pow(p, r) * pow((1-p), k)

使用ibeta_derivative函數實現:

(p/(r + k)) * ibeta_derivative(r, static_cast<RealType>(k+1), p) 在這裡使用函數 ibeta_derivative ,因為它已經針對最低可能誤差( lowest possible error)進行了優化- 實際上,這僅僅是對不完全β函數的部分內部實現進行了一個簡單的包裝。

cdf

使用下面的關係:

cdf = Ip(r, k+1) = ibeta(r, k+1, p)

= ibeta(r, static_cast<RealType>(k+1), p)

cdf 補集(complement)

使用下面的關係:

1 - cdf = Ip(k+1, r)

= ibetac(r, static_cast<RealType>(k+1), p)

分位點(quantile)

ibeta_invb(r, p, P) - 1

補集的分位點(quantile from the complement )

ibetac_invb(r, p, Q) -1)

均值(mean)

r(1-p)/p

方差(variance)

r (1-p) / p * p

眾數(mode)

floor((r-1) * (1 - p)/p)

偏斜(skewness)

(2 - p) / sqrt(r * (1 - p))

峰態(kurtosis)

6 / r + (p * p) / r * (1 - p )

峰態超越(kurtosis excess)

6 / r + (p * p) / r * (1 - p ) -3

參數估算成員函數(parameter estimation member functions)

find_lower_bound_on_p

ibeta_inv(successes, failures + 1, alpha)

find_upper_bound_on_p

ibetac_inv(successes, failures, alpha) plus see comments in code.

find_minimum_number_of_trials

ibeta_inva(k + 1, p, alpha)

find_maximum_number_of_trials

ibetac_inva(k + 1, p, alpha)

實現註解:


PrevUpHomeNext