Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext
查找均值和標準差實例(Find mean and standard deviation example)

首先我們包含一些頭文件來訪問正態分佈,用來查找位置(location)和尺度(scale)的函數(以及標準輸出)。

#include <boost/math/distributions/normal.hpp> // normal_distribution
  using boost::math::normal; // typedef 提供的缺省類型為 double.
#include <boost/math/distributions/cauchy.hpp> //  cauchy_distribution
  using boost::math::cauchy; // typedef 提供的缺省類型為 double.
#include <boost/math/distributions/find_location.hpp>
  using boost::math::find_location;
#include <boost/math/distributions/find_scale.hpp>
  using boost::math::find_scale;
  using boost::math::complement;
  using boost::math::policies::policy;

#include <iostream>
  using std::cout; using std::endl; using std::left; using std::showpoint; using std::noshowpoint;
#include <iomanip>
  using std::setw; using std::setprecision;
#include <limits>
  using std::numeric_limits;

使用函數 find_location 和函數 find_scale 來滿足分配和測量規範( dispensing and measurement specifications)

考慮來自於 K Krishnamoorthy, Handbook of Statistical Distributions with Applications, ISBN 1-58488-635-8, (2006) p 126, example 10.3.7的一個例子。

"一台機器用於在每個盒子中裝3 kg牛肉(ground beef)。經過很長一段時間發現每個盒子中的牛肉平均重量為3kg且標準差為0.1kg。假定將牛肉裝進例子的過程為正態分佈。"

我們通過用給定的參數構造一個正態分佈來開始:

double mean = 3.; // kg
double standard_deviation = 0.1; // kg
normal packs(mean, standard_deviation);

我們可以找到重量大於3.1kg的盒子的百分數。

double max_weight = 3.1; // kg
cout << "Percentage of packs > " << max_weight << " is "
<< cdf(complement(packs, max_weight)) * 100. << endl; // P(X > 3.1)

我們可能想要確定95%以上的盒子的重量超過最小的給定重量,那麼我們想均值滿足 : P(X < 2.9) = 0.05。

使用均值為 3 kg,我們可以估算沒有滿足指定的2.9kg的盒子的百分數。

double minimum_weight = 2.9;
cout <<"Fraction of packs <= " << minimum_weight << " with a mean of " << mean
  << " is " << cdf(complement(packs, minimum_weight)) << endl;
// fraction of packs <= 2.9 with a mean of 3 is 0.841345

結果為 0.84 - 小於目標百分數 0.95。如果我們想要95% 的盒子超過最小重量,均值應當設為多少?

使用這本書中提供的 KK StatCalc 程序,在第126頁給出的方法計算的值為 3.06449。

我們可以通過使用一個安全的邊緣均值(safety margin mean)3.06449來構造一個新的稱之為『'xpacks'』的分佈來證實這一點:

double over_mean = 3.06449;
normal xpacks(over_mean, standard_deviation);
cout << "Fraction of packs >= " << minimum_weight
<< " with a mean of " << xpacks.mean()
  << " is " << cdf(complement(xpacks, minimum_weight)) << endl;
// fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005

使用數學工具,我們可以直接計算要求的均值:

double under_fraction = 0.05;  //  95% 的盒子超過最小重量均值 - sd = 2.9
double low_limit = standard_deviation;
double offset = mean - low_limit - quantile(packs, under_fraction);
double nominal_mean = mean + offset;
// mean + (mean - low_limit - quantile(packs, under_fraction));

normal nominal_packs(nominal_mean, standard_deviation);
cout << "Setting the packer to " << nominal_mean << " will mean that "
  << "fraction of packs >= " << minimum_weight
  << " is " << cdf(complement(nominal_packs, minimum_weight)) << endl;
// 設置盒子重量為 3.06449 將意味著 >= 2.9 的盒子的百分比為 0.95

這種計算一般化為稱為find_location的自由函數。

為了使用這個函數我們需要

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

然後使用函數 find_location來查找 safe_mean,並構造一個新的正態分佈稱之為 'goodpacks'。

double safe_mean = find_location<normal>(minimum_weight, under_fraction, standard_deviation);
normal good_packs(safe_mean, standard_deviation);

使用與上面相同的證實(confirmation):

cout << "Setting the packer to " << nominal_mean << " will mean that "
  << "fraction of packs >= " << minimum_weight
  << " is " << cdf(complement(good_packs, minimum_weight)) << endl;
// 設置盒子重量為 3.06449 將意味著 >= 2.9 的盒子的百分比為 0.95

使用柯西-洛倫茨分佈取代正態分佈(Using Cauchy-Lorentz instead of normal distribution)

在檢驗了大量的盒子重量的分佈之後 ,我們可能會決定,畢竟,正態分佈的假設是不公平的。我們可能會覺得這種情況更適合使用柯西分佈(Cauchy Distribution)。這個分佈有更寬的範圍('wings'),所以,儘管大部分的數值與正常值(normal value)相比更接近於均值(mean value),仍然有更多的值與正常值(normal value)相比離均值更遠。

這可能是因為一塊比正常重量更重的牛肉要麼被裝入盒子中要麼沒有被裝入盒子中。

我們首先使用原始的均值和標準差來生成一個柯西分佈(Cauchy Distribution) ,並估算少於指定的最小重量的盒子的百分數

cauchy cpacks(mean, standard_deviation);
cout << "Cauchy Setting the packer to " << mean << " will mean that "
  << "fraction of packs >= " << minimum_weight
  << " is " << cdf(complement(cpacks, minimum_weight)) << endl;
// 設置盒子重量為 3 將意味著 >= 2.9 的盒子的百分比為 0.75

注意:只有更少的盒子滿足這個指定的標準,僅是 75% 而不是 95%。現在我們可以使用柯西分佈( cauchy distribution)來取代正態分佈(normal distribution)作為模板參數來重複調用函數 find_location。

double lc = find_location<cauchy>(minimum_weight, under_fraction, standard_deviation);
cout << "find_location<cauchy>(minimum_weight, over fraction, standard_deviation); " << lc << endl;
// find_location<cauchy>(minimum_weight, over fraction, packs.standard_deviation()); 3.53138

注意:現在的safe_mean值需要設置得更高一些, 3.53138 而不是 3.06449,因此我們將獲得更少的利益。

再一次證實滿足指定要求的百分數與預期的是一樣的。

cauchy goodcpacks(lc, standard_deviation);
cout << "Cauchy Setting the packer to " << lc << " will mean that "
  << "fraction of packs >= " << minimum_weight
  << " is " << cdf(complement(goodcpacks, minimum_weight)) << endl;
// 設置盒子重量為 3.53138 將意味著 >= 2.9 的盒子的百分比為 0.95

最後,我們可以估算一個更緊(tighter)的規範的效果,99%的盒子滿足這個規範。

cout << "Cauchy Setting the packer to "
  << find_location<cauchy>(minimum_weight, 0.99, standard_deviation)
  << " will mean that "
  << "fraction of packs >= " << minimum_weight
  << " is " << cdf(complement(goodcpacks, minimum_weight)) << endl;

將盒子的重量設為 3.13263 將意味著重量大於2.9的盒子的百分數為0.099。但均值損失將會由每盒0.0644 增大到 0.133 kg而加倍。

當然,這種計算並不局限於裝牛肉的盒子( packs of meat),它可以用於計量(dispensing)任何東西,就像任何的測量一樣,這種計算還可以應用於「虛擬的」物質。

唯一要注意的是:這種計算假定標準差(尺度)是已知的且不確定性相當低,在實際中這很難確定。並且這個分佈是定義良好的正態分佈(Normal Distribution)柯西分佈(Cauchy Distribution),或者其它的分佈。

如果某人只是簡單地測量大量的盒子,那麼就可能測量成百上千的盒子的重量。在有一個合理的「自由度」(degree of freedom)的情況下,標準差的置信區間就不會太寬,對於大量的觀測,典型的情況是大約 + 以及 - 10% 。

對於其它的應用,當做大量觀測很難或是很昂貴的時候,置信區間就會很寬。

參考標準差的置信區間(Confidence Intervals on the standard deviation) 來查看一個可運行的例子,chi_square_std_dev_test.cpp 估算這些區間。

改變尺度或標準差(Changing the scale or standard deviation)

另一方面,我們可以投資一個標準差更低的更好的(精確度更高)的打包機。

這可能會花費更多,但是會減少為了達到規範而不得不「贈送(give away)」的數量。

為了估算新的打包機到底有多好(標準差到底有多小),我們必須使得5%的分位點位於 under_weight limit, 2.9 處。

double p = 0.05; // wanted p th quantile.
cout << "Quantile of " << p << " = " << quantile(packs, p)
  << ", mean = " << packs.mean() << ", sd = " << packs.standard_deviation() << endl;

Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1

對於當前的打包機 (mean = 3, sd = 0.1), 5% 分位點在 2.8551 kg處,比我們2.9kg的目標稍低一些。所以,標準差將會變小。

讓我們通過假設將其對分(現在為0.1)為標準差為0.05kg開始。

normal pack05(mean, 0.05);
cout << "Quantile of " << p << " = " << quantile(pack05, p)
  << ", mean = " << pack05.mean() << ", sd = " << pack05.standard_deviation() << endl;
// Quantile of 0.05 = 2.91776, mean = 3, sd = 0.05

cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
  << " and standard deviation of " << pack05.standard_deviation()
  << " is " << cdf(complement(pack05, minimum_weight)) << endl;
// 重量大於 2.9 ,均值為 3,標準差為0.05的包的百分比為 0.97725

所以0.05 是一個非常好的估計,但對於2.9的目標稍高了一點,所以標準差會稍大一點。所以,我們可以進行更多的估計以使得結果更接近一些,假設把標準差提高到 0.06 kg,構造另一個稱之為 pack06的新的分佈。

normal pack06(mean, 0.06);
cout << "Quantile of " << p << " = " << quantile(pack06, p)
  << ", mean = " << pack06.mean() << ", sd = " << pack06.standard_deviation() << endl;
// Quantile of 0.05 = 2.90131, mean = 3, sd = 0.06

cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
  << " and standard deviation of " << pack06.standard_deviation()
  << " is " << cdf(complement(pack06, minimum_weight)) << endl;
// 重量大於 2.9 ,均值為 3,標準差為0.06的包的百分比為 0.95221

現在我們的計算結果真的非常接近了,為了更好的完成這個工作,我們可能需要使用根查找算法( root finding method),例如已經提供的工具,參考 不使用導數查找根(Root Finding Without Derivatives)

但在這個分佈(正態分佈)情形中,我們應當更靈活並進行一個直接的計算。

我們要求的界限是 minimum_weight = 2.9 kg,通常稱作隨機變量z。對於一個標準正態分佈,那麼概率 p = N((minimum_weight - mean) / sd)。

我們想要查找滿足界限的標準差,使得第p個分位點位於點z處(minimum_weight)。在這種情況下,0.05 (5%) 分位點位於 2.9 kg 處,當均值為 3 kg,確保 0.95 (95%) 的盒子的重量在最小重量之上。

重新整理,我們可以直接計算要求的標準差:

normal N01; // 均值為0且標準差為1的標準正態分佈.
p = 0.05;
double qp = quantile(N01, p);
double sd95 = (minimum_weight - mean) / qp;

cout << "For the "<< p << "th quantile to be located at "
  << minimum_weight << ", would need a standard deviation of " << sd95 << endl;
// 為了使 0.05th 分位點位於 2.9處, 需要的標準差為 0.0607957

現在我們可以為這個更好的打包機構造一個更好的分佈pack95,並且檢查我們的分佈將會滿足規範。

normal pack95(mean, sd95);
cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
  << " and standard deviation of " << pack95.standard_deviation()
  << " is " << cdf(complement(pack95, minimum_weight)) << endl;
// 重量大於 2.9 ,均值為 3,百分比為 0.95

這種計算被泛化為函數 find_scale,如下面顯示的那樣,給出相同的標準差。

double ss = find_scale<normal>(minimum_weight, under_fraction, packs.mean());
cout << "find_scale<normal>(minimum_weight, under_fraction, packs.mean()); " << ss << endl;
// find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957

如果我們已經定義了一個 over_fraction,或者必須超過規範的百分比:

double over_fraction = 0.95;

並 (錯誤地) 寫作

double sso = find_scale<normal>(minimum_weight, over_fraction, packs.mean());

基於缺省的策略,我們將會得到下面的出錯消息:

Message from thrown exception was:
   Error in function boost::math::find_scale<Dist, Policy>(double, double, double, Policy):
   Computed scale (-0.060795683191176959) is <= 0! Was the complement intended?

這會返回一個負的 標準差 - 顯然不可能。 概率應當為 1 - over_fraction,而不是 over_fraction,因此:

double ss1o = find_scale<normal>(minimum_weight, 1 - over_fraction, packs.mean());
cout << "find_scale<normal>(minimum_weight, under_fraction, packs.mean()); " << ss1o << endl;
// find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957

但請注意使用 '1 - over_fraction' - 將會導致 精度損失,尤其是當函數非常接近1. 。在這種情況下(非常常見的情況),我們將會使用補集(complements) 來替代,給出最精確的結果。。

double ssc = find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean()));
cout << "find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean())); " << ssc << endl;
// find_scale<normal>(complement(minimum_weight, over_fraction, packs.mean())); 0.0607957

注意,我們的猜測值非常接近於精確值 0.060795683191176959。

我們可以再一次證實我們的預測:

normal pack95c(mean, ssc);
cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " << mean
  << " and standard deviation of " << pack95c.standard_deviation()
  << " is " << cdf(complement(pack95c, minimum_weight)) << endl;
// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.0607957 is 0.95

注意這兩個簡單的問題:

且 / 或

是實際上非常常見的問題。

測量牛肉重量可以被其它的任何東西的測量替代,藥片含量,阿波羅登陸火箭點火(Apollo landing rocket firing),X放射線治療......

在分配(dispensing)或測量的不確定性中,尺度(scale)將會是一個變化量。

參考find_mean_and_sd_normal.cpp 查看完整代碼& 附加的程序輸出。


PrevUpHomeNext