Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext
二項拋擲硬幣(Binomial Coin-Flipping)實例

伯努利過程(Bernoulli process) 的一個例子是硬幣拋擲(coin flipping)。在這樣一個序列中的一個變量被稱作伯努利變量(Bernoulli variable)。

這個例子顯示了使用二項分佈(binomial distribution)來預測當拋擲一個硬幣時正面和反面(heads and tails)的概率。

正確答案的數量 (比如說,硬幣正面), X,的分佈為一個試驗次數n=10且結果為硬幣正面的概率(成功分數)為 p = 0.5 (一個 公平('fair)' 的硬幣)的二項分佈的變量的分佈。

(我們所使用的硬幣假定是公平的(fair),但是我們可以很容易地將成功分數(success fraction)的值從0.5更改為其它的值來模擬一個不公平的(unfair)硬幣,比如說對於一個背面粘有口香糖的硬幣,落下的時候更可能為反面朝下而正面朝上,將其成功分數改為0.6).

首先我們需要一些include和using 語句來使得可以使用二項分佈(binomial distribution),一些標準輸入和輸出,然後我們可以開始了:

#include <boost/math/distributions/binomial.hpp>
  using boost::math::binomial;

#include <iostream>
  using std::cout;  using std::endl;  using std::left;
#include <iomanip>
  using std::setw;

int main()
{
  cout << "Using Binomial distribution to predict how many heads and tails." << endl;
  try
  {

參考使用try catch 瞭解為什麼一直使用try catch是一個好主意。

首先,以拋擲硬幣的次數和成功分數(success fraction)為參數來構造一個二項分佈:

const double success_fraction = 0.5; // = 50% = 1/2 for a 'fair' coin.
int flips = 10;
binomial flip(flips, success_fraction);

cout.precision(4);

然後是一些使用二項矩(Binomial moments)的例子 (並顯示參數)。

cout << "From " << flips << " one can expect to get on average "
  << mean(flip) << " heads (or tails)." << endl;
cout << "Mode is " << mode(flip) << endl;
cout << "Standard deviation is " << standard_deviation(flip) << endl;
cout << "So about 2/3 will lie within 1 standard deviation and get between "
  <<  ceil(mean(flip) - standard_deviation(flip))  << " and "
  << floor(mean(flip) + standard_deviation(flip)) << " correct." << endl;
cout << "Skewness is " << skewness(flip) << endl;
// 二項分佈的偏斜(Skewness)為0(對稱的)
// 如果 success_fraction 如果是1/2,
// 例如, 當拋擲一個'公平的(fair)' 硬幣.
cout << "Skewness if success_fraction is " << flip.success_fraction()
  << " is " << skewness(flip) << endl << endl; // 對於一個'公平的(fair)' 硬幣,預期結果為0.

現在我們顯示多個正面概率的預測:

cout << "For " << flip.trials() << " coin flips: " << endl;
cout << "Probability of getting no heads is " << pdf(flip, 0) << endl;
cout << "Probability of getting at least one head is " << 1. - pdf(flip, 0) << endl;

當我們想要計算某個範圍內的概率時,我們可以將函數PDF的值求和:

cout << "Probability of getting 0 or 1 heads is "
  << pdf(flip, 0) + pdf(flip, 1) << endl; // sum of exactly == probabilities

或者我們可以使用函數CDF。

cout << "Probability of getting 0 or 1 (<= 1) heads is " << cdf(flip, 1) << endl;
cout << "Probability of getting 9 or 10 heads is " << pdf(flip, 9) + pdf(flip, 10) << endl;

注意,使用

cout << "Probability of getting 9 or 10 heads is " << 1. - cdf(flip, 8) << endl;

比使用補集(complement)的精度稍低一些。

cout << "Probability of getting 9 or 10 heads is " << cdf(complement(flip, 8)) << endl;

因為減法可能涉及消去錯誤(cancellation error),又因為函數cdf(complement(flip, 8)) 內部並不使用減法,因此不會產生這個問題。

為了獲得在某個範圍內顯示硬幣正面的概率,我們可以將每個顯示硬幣正面的pdf函數的值累加起來:

cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "
  //  P(X == 4) + P(X == 5) + P(X == 6)
  << pdf(flip, 4) + pdf(flip, 5) + pdf(flip, 6) << endl;

但這與使用函數cdf相比,效率要低一些。

cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "
  // P(X <= 6) - P(X <= 3) == P(X < 4)
  << cdf(flip, 6) - cdf(flip, 3) << endl;

當然,對於一個大的範圍,類似於 3 到 7

cout << "Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is "
  // P(X <= 7) - P(X <= 2) == P(X < 3)
  << cdf(flip, 7) - cdf(flip, 2) << endl;
cout << endl;

最後,打印兩個表的概率結果:剛好 顯示某一數目的硬幣正面的概率和至少 顯示某一數目的硬幣正面的概率。

// 打印剛好顯示某一數目的硬幣正面的概率表
cout << "Probability of getting exactly (==) heads" << endl;
for (int successes = 0; successes <= flips; successes++)
{ // 假設成功的意思是說顯示一次硬幣的正面 (或者等價的,成功為顯示一次硬幣的反面).
  double probability = pdf(flip, successes);
  cout << left << setw(2) << successes << "     " << setw(10)
    << probability << " or 1 in " << 1. / probability
    << ", or " << probability * 100. << "%" << endl;
} // 對於 i
cout << endl;

// 將顯示0次正面與顯示多達10次正面的概率製表。
cout << "Probability of getting upto (<=) heads" << endl;
for (int successes = 0; successes <= flips; successes++)
{ // 假設成功的意思是說顯示一次硬幣的正面
  // (或者等價的,成功為顯示一次硬幣的反面).
  double probability = cdf(flip, successes); // P(X <= heads)
  cout << setw(2) << successes << "        " << setw(10) << left
    << probability << " or 1 in " << 1. / probability << ", or "
    << probability * 100. << "%"<< endl;
} // 對於 i

最後 (0 to 10 heads)的概率必須,當然, 是 100% 的概率。

}
catch(const std::exception& e)
{
  //

包含 try & catch 永遠都是最基本的,因為對於超出範圍的定義域或類似於數值溢出產生的錯誤,缺省的策略是拋出異常。

缺少 try & catch ,程序將會退出(abort),然而下面的從異常中顯示的消息將會是查找問題原因的一個很好的線索。

  std::cout <<
    "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
}

參考binomial_coinflip_example.cpp 查看完整代碼,程序的輸出類似於下面:

Using Binomial distribution to predict how many heads and tails.
From 10 one can expect to get on average 5 heads (or tails).
Mode is 5
Standard deviation is 1.581
So about 2/3 will lie within 1 standard deviation and get between 4 and 6 correct.
Skewness is 0
Skewness if success_fraction is 0.5 is 0

For 10 coin flips:
Probability of getting no heads is 0.0009766
Probability of getting at least one head is 0.999
Probability of getting 0 or 1 heads is 0.01074
Probability of getting 0 or 1 (<= 1) heads is 0.01074
Probability of getting 9 or 10 heads is 0.01074
Probability of getting 9 or 10 heads is 0.01074
Probability of getting 9 or 10 heads is 0.01074
Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6562
Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6563
Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is 0.8906

Probability of getting exactly (=) heads
0      0.0009766  or 1 in 1024, or 0.09766%
1      0.009766   or 1 in 102.4, or 0.9766%
2      0.04395    or 1 in 22.76, or 4.395%
3      0.1172     or 1 in 8.533, or 11.72%
4      0.2051     or 1 in 4.876, or 20.51%
5      0.2461     or 1 in 4.063, or 24.61%
6      0.2051     or 1 in 4.876, or 20.51%
7      0.1172     or 1 in 8.533, or 11.72%
8      0.04395    or 1 in 22.76, or 4.395%
9      0.009766   or 1 in 102.4, or 0.9766%
10     0.0009766  or 1 in 1024, or 0.09766%

Probability of getting upto (<) heads
0         0.0009766  or 1 in 1024, or 0.09766%
1         0.01074    or 1 in 93.09, or 1.074%
2         0.05469    or 1 in 18.29, or 5.469%
3         0.1719     or 1 in 5.818, or 17.19%
4         0.377      or 1 in 2.653, or 37.7%
5         0.623      or 1 in 1.605, or 62.3%
6         0.8281     or 1 in 1.208, or 82.81%
7         0.9453     or 1 in 1.058, or 94.53%
8         0.9893     or 1 in 1.011, or 98.93%
9         0.999      or 1 in 1.001, or 99.9%
10        1          or 1 in 1, or 100%


PrevUpHomeNext