做品管的，大都使用過管制圖。P管制圖和NP管制圖係依據二項分布之假設，而 C管制圖和U管制圖則依據卜瓦松分布之假設繪製而成（關於管制圖與這些統計分布的關聯性描述，以後有時間再開另篇文章談，這須用到數學公式說明）。二項分布和卜瓦松分布在特定條件下可藉由常態分布來近似，故將平均數+3個標準差作為管制圖管制上限及平均數-3個標準差作為管制圖管制下限。至於為何是3個標準差，這議題日後有機會再開另一篇文章來描述。
這篇文章的談論主角是「卜瓦松分布」(順便為日後的卜瓦松迴歸模型議題留伏筆)，既然要談論它，它的知識來源，包含伯努利分布和二項分布就一起在此篇文章說明。這些分布描述的對象都屬離散型資料（所以能見度很高的幾何分布也拉進來一起談），例如藥物不良反應件數指標，本月23件，上月為18件。18～23間無法無限制分割，例如分割線落在21.56件(該月是幾件就幾件，不會有21.56件這種帶小數位的件數)；相對來說，本次監測血壓值為88.5mmHg，上次監測為92.4mmHg，88.5～92.4間可無限制分割，例如分割線落在91.026mmHg（測出來的血壓值可以是這樣，如果機器夠靈敏的話），這種可無限制分割下去的資料型態就叫連續型資料。平日小弟傳授管制圖課程，不會帶入以上機率分布概念(學員會倒一片，且不學不會防礙管制圖類別選擇，但因無機率分配概念，故小弟提的判圖原則須謹記於心)，至於離散型和連續型資料判別，則會於課程中加以說明(這有助於管制圖類別選擇)。而下表內的統計公式，只是本文章的裝潢物件(這樣看起來很有學問，是不是？)，不是本篇文章的說明重點。
好，開始切入主題(小弟會儘量不提公式)。
伯努利分布 (Bernoulli distribution)
只有2種可能(p機率+q機率＝1)結果的單次隨機試驗稱為伯努利試驗。
例如拋10元硬幣，人頭朝上(假設機率p=0.5)記為1，人頭朝下(假設機率q=0.5)記為0，接著就可說拋1次硬幣就是一個p=q=1/2的伯努利試驗。，它是最簡單的離散分布，也是構成其它複雜離散分布的基礎。
但現實生活中，p和q不太會相等。接下來用SAS迴圈模擬出10萬筆伯努利試驗後資料(p=0.6)。
/*伯努利分布*/
data bern;
do n=1 to 100000;
p=0.6;
x=rand("bernoulli",p);
output;end;
keep x;
run;
proc sgplot;
title h=2 "伯努利分布";
histogram x / binwidth=0.5 binstart=0;
density x / type=normal legendlabel="常態分布" lineattrs=(pattern=solid color=red);
run;
proc means mean;
run;
二項分布(Binomial distribution)
N次重複獨立的伯努利試驗中，獲得n次成功的離散機率分布。
例如1次拋50枚硬幣(視為重複50次，但實際上此50次視為1次試驗)，來看人頭朝上的枚數，此時已產生1個人頭朝上的率值，此率值就是二項分布。第1次拋有23枚人頭朝上、第2次拋有19枚人頭朝上、第3次拋有23枚人頭朝上、第4次拋有28枚人頭朝上、第5次拋有15枚人頭朝上、..X軸依序為人頭朝上的次數，由1到50，Y軸為出現次數，此直方圖即為N=50, p=0.5的二項分布。
所以可以說伯努利分布是二項分布N=1的特例，所以二項分布的期望值和變異數都是伯努利分布的期望值和變異數的N倍(因為重複N次)，當重複次數N=1時，服從二項分布X~Bin(1,p)，它等於伯努利分布X~Bern(p)，也等於區間[0,1]上的離散均勻分布Unif{0,1}(有沒有覺得用了一些鬼符號讓一般人看不懂就覺得很專業？小弟會遵守諾言「儘量不提公式」並「適可而止」)。
用SAS迴圈模擬出1000次試驗，每一次試驗重複20次，p=0.5的二項分布資料。
/*二項分布*/
data bin;
do t=1 to 1000;
p=0.5; n=20;
x=rand("binomial",p,n);
output;end;
keep x;
run;
proc sgplot;
title h=2 "二項分布";
histogram x / binwidth=1 binstart=1;
xaxis grid minorgrid values=(1 to 20 by 2);
yaxis min=1 max=20 integer grid minorgrid;
run;
proc means mean;
run;
proc means mean;
run;
當重複的次數N很大時，且p不接近1或0，則de Moivre - Laplace中央極限定理(為中央極限定理的初始版)證明，該二項分布會近似於常態分布，一般當N>=30時，逼近常態的效果較佳。
來驗證上面橘色內容：
用SAS迴圈模擬出10萬次試驗，每一次試驗重複1萬次，p=0.5的二項分布資料。
/*二項分布逼近常態分布*/
data bin_to_normal;
do t=1 to 100000; /*試驗次數設10萬次*/
n=10000; p=0.5; /*每一試驗含有1萬次重複，這樣的重複次數N很大了。p為0.5,這樣不接近0或1*/
x=rand("binomial",p,n);
output;end;
keep x;
run;
proc sgplot;
title h=2 "二項分布逼近常態分布";
histogram x ;
density x / type=normal legendlabel="常態分布" lineattrs=(pattern=solid color=red);
density x / type=kernel legendlabel="二項分布,每一試驗重複1萬次, p=0.5" lineattrs=(pattern=solid color=green);
xaxis display=(nolabel) ;
run;
紅綠線幾乎重疊了，不用再驗證了。
卜瓦松分布(Poisson distribution)
終於輪朕上場了。先不說明，承橘色內容，但p接近1或0時，該二項分布會近似於卜瓦松分布。
用SAS迴圈模擬出10萬次試驗，每一次試驗重複1萬次，p=0.05的二項分布資料。
/*二項分布逼近卜瓦松分布*/
data bin_to_poi;
do t=1 to 100000; /*試驗次數設10萬次*/
n=10000; p=0.05; /*每一試驗含有1萬次重複，這樣的重複次數N很大了。p為0.05也很小了*/
x=rand("binomial",p,n);
m=n*p;
x2=rand("poisson",m);
output;end;
keep x x2;
run;
proc sgplot;
title h=2 "二項分布逼近卜瓦松分布";
density x / type=kernel legendlabel="二項分布,每一試驗重複1萬次, p=0.05" lineattrs=(pattern=solid color=red);
density x2 / type=kernel legendlabel="卜瓦松分布,m=np=50" lineattrs=(pattern=solid color=green);
xaxis display=(nolabel) ;
run;
「現主時」，二項分布Bin(n,p)逼近lamda=np (先不管它是什麼鬼東西)的卜瓦松分布。
好，開始解釋卜瓦松分布。
它是單位時間內或範圍內獨立事件發生次數的機率分布(而指數分布則是獨立事件發生時間間隔的機率分布，常用於排隊統計)，它的前提是事件間不能有關聯。
日常生活中，有大量事件發生在單位時間內是有固定次數規律的。例如某醫院平均每1小時接生3個嬰兒、醫院每一台自助繳費機在一定時間內服務的民眾次數或機器出現故障次數等。這種有固定發生頻率的事件就服從卜瓦松分布。固然可以估計在某個時段內發生的次數，但永遠無法知道發生的具體時間。
卜瓦松的公式來自於二項分布公式，當預期成功的試驗次數當作自變量，而不再是特定機率p下的試驗次數N，且當試驗次數足夠大(微積分後)時，其分布機率公式如下：
來手算數學，例如高速公路每5秒過2部車，則lamda=2 (即每單位時間平均有2部車)
(1)5 秒內恰有一車通過的機率？
5秒，所以1個時間單位，所以t=1
一車通過，所以n=1
(2)10 秒內恰有一車通過的機率？
10秒，所以2個時間單位，所以t=2
一車通過，所以n=1
10秒，所以2個時間單位，所以t=2
沒車通過，所以n=0
(4) 一分鐘平均有幾輛車通過？
1分鐘有12個5秒，即t=12
/*卜瓦松分布平均件數*/
data poi_n;
do x=0 to 25;
p1=pdf("poisson",x,1);
p4=pdf("poisson",x,4);
p10=pdf("poisson",x,10);
output;end;
run;
proc sgplot;
title h=2 "卜瓦松分布平均件數";
series x=x y=p1 / markers markerattrs=(color=red) lineattrs=(pattern=solid color=red) legendlabel="lamda=1";
series x=x y=p4 / markers markerattrs=(color=green) lineattrs=(pattern=solid color=green) legendlabel="lamda=4";
series x=x y=p10 / markers markerattrs=(color=blue) lineattrs=(pattern=solid color=blue) legendlabel="lamda=10";
keylegend / location=inside position=topright across=1;
yaxis label="機率值";
run;
轉成累積型機率分布後
/*卜瓦松分布平均件數_累積型*/
data poi_n_cum;
do x=0 to 25;
p1=cdf("poisson",x,1);
p4=cdf("poisson",x,4);
p10=cdf("poisson",x,10);
output;end;
run;
proc sgplot;
title h=2 "卜瓦松分布平均件數";
series x=x y=p1 / markers markerattrs=(color=red) lineattrs=(pattern=solid color=red) legendlabel="lamda=1";
series x=x y=p4 / markers markerattrs=(color=green) lineattrs=(pattern=solid color=green) legendlabel="lamda=4";
series x=x y=p10 / markers markerattrs=(color=blue) lineattrs=(pattern=solid color=blue) legendlabel="lamda=10";
keylegend / location=inside position=topright across=1;
yaxis label="機率值";
run;
前面提到「某醫院平均每1小時接生3個嬰兒」例子，換用SAS來計算，這次不手工計算。
(1)1小時內生產0個嬰兒的機率p0
(2)1小時內生產1個嬰兒的機率p1
(3)1小時內至少生產2個以上嬰兒的機率p
(3)的解法為
該機率為 1-（1小時內接生1人的機率p1）－（1小時內接生0人的機率p0）
/*某醫院平均每1小時接生3個嬰兒專案*/
data occur;
lamda=3;
t1=1;
m1=lamda*t1;
n1=1;
p1=pdf("poisson",n1,m1); /*1小時內至少生產1個嬰兒的機率*/
t0=1;
m0=lamda*t0;
n0=0;
p0=pdf("poisson",n0,m0); /*1小時內沒有接生嬰兒的機率*/
p=1-p1-p0; /*1小時內至少生產2個嬰兒的機率*/
put p=;
if (p>0.9973) then put "Note: 極可能發生(>99.7％)";
else if (p>0.954) then put "Note: 很可能發生(>95％)";
else if (p>0.683) then put "Note: 有可能發生(>68％)";
else if (p<0.003) then put "Note: 不可能發生(<0.3％)";
run;
Note: 卜瓦松迴歸模型
卜瓦松機率曲線圖也可透過統計模型來適配(就像線性迴歸統計模型那樣)，並評估此卜瓦松迴歸模型(為廣義線性模型Generalized linear model項下的其中一分支)內，每一自變項對依變項的影響程度(此議題以後可以另寫文章來介紹)。SAS的Proc genmod可以執行此類分析，且分析前會對卜瓦松分布進行評估。一旦發現過度離散現象Overdispersion (原因常為數據間不獨立、尚有影響依變項的重要自變項未入模型中、數據存在異常值。卜瓦松分布有一重要特徵，其平均數等於或接近變異數，若出現該現象，此特徵通常不存在)，則適配前須對該現象進行校正，校正的方法有多種，或改對數據執行「負二項迴歸分析」。
而且，它只是視訊課程（不過很有誠意的連負二項、零膨脹模型都講了）。
幾何分布(Poisson distribution)
它是唯一的離散型無記憶隨機分布，可視為指數分布的離散模式。它是N次重複的伯努利試驗中第一次試驗成功之前的失敗次數，或者獲得1次成功所需要的試驗次數所服從的統計機率分布，它也是負二項分布的基礎。
例如一直投骰子，直到點數為6的投擲次數，則該投擲次數服從p=1/6的幾何分布。
用SAS迴圈模擬成功投到點數6的次數達1萬次，統計每一次成功前的失敗次數。
/*幾何分布*/
data geo;
n=10000;
t=0;
do i=1 to n;
x=floor(rand("uniform")*6)+1;
if x=6 then do; x=t; t=0;
output;end;
t=t+1;end;
keep x;
run;
proc sgplot;
title h=2 "幾何分布";
histogram x ;
density x / type=normal legendlabel="常態分布" lineattrs=(pattern=solid color=red);
density x / type=kernel legendlabel="幾何分布" lineattrs=(pattern=solid color=green);
xaxis display=(nolabel) values=(-20 to 60 by 10) ;
run;
還是有些人要投20次以上才投的出點數6，可謂「天選之人」。
這篇文章充份展示了各種統計機率分布之間在某些條件下相互轉化的數學規律。
沒有留言:
張貼留言