前言:

前一篇, 我們利用霍納, 秦九韶的算法寫出了可以快速計算各種多項式函數的 FORTH 程式碼. 這一篇, 我們就來進階應用吧! 原來啊, 多項式函數還有個很有用的用途, 一些超越函數或是特殊函數, 有些是可以用多項式函數來做個多階的逼近跟取值的, (冪級數power series, 無窮階的多項式級數展開...)

而我們的 FORTH系統雖然在 ANS FORTH 已經定義了一些重要的數學函數. 但其實, 有一些數學函數其實也很重要, 例如鼎鼎大名的高斯分佈機率函數! 因為在工程的領域, 或是實驗, 統計的計算, 分佈的估計, 假設檢定. 樣樣都離不開這個重要的函數的!

所以來建構這個連 Excel 都內建的常態分佈機率及其機率密度函數來提供給我們的 FORTH 系統使用吧! 這樣就更如虎添翼囉, 有這個函數的加持, 把 FORTH 當統計軟體, 或統計語言來使用, 更是一點都不奇怪, 整個 FORTH 系統更是威上加威囉, 會更方便來解決任何科學/統計/工程上的任何問題.

這篇建構完畢整個常態分佈函數群後, 會再提供個簡單的示範, 示範如何用我們所建構的函數來計算估計並解決我們所遇到的工程問題. (水果商對於蘋果出貨分類的分佈跟預估!) 讓大家瞧瞧看, 如何把 FORTH 當 R 語言使用!

 

常態分佈(高斯分佈)機率密度函數及其機率分佈:

這個由平均值跟標準差兩個參數所標定的常態分佈真的是所有理工科的學生都再熟悉不過了, 它的重要性緣由於於數學上的中央極限定理. 可以證明, 今天在一個未知現象未知分佈的一個母群體中, 即使它本身的分佈不是常態分佈. 但只要我們對它多次的取樣, 最後取樣出來每次所計算的平均值所形成的分佈, 隨著取樣數目的增加, 會慢慢地變成是一個常態分佈. (平均值樣本數目超過30個以上) 

也因此常態分佈廣泛地出現在所有的自然科學現象, 物品的生產, 製造之中.

這裡有些精闢的演示跟討論, 常態分佈本質. 根據他的論點, 因為抽樣的過程, 本身就是一種二項式分佈(巴斯卡三角形). 純數學上可以證明, 當項數目很大(無窮大)的時候, 二項式分佈會等於常態分佈. 所以中央極限定理的出現也就不足為奇囉!

這裡有篇用 R 語言演示中央極限定理的範例, 可以看得很清楚. 各種不同的分佈, 透過抽樣之後, 它們的樣本平均值的分佈會漸漸地變成常態分佈.

 

製程 SPC 控制:

也因為如此, 當我們工廠生產時, 雖然母群體的分佈我們是不知道的. 但是因為我們都是透過抽樣來做量測. 所以抽樣樣本的平均值, 如果抽樣的數目夠多(大於30以上)的話, 由中央極限定理, 它們會遵守常態分佈.

所以 SPC 的制定方式就是這樣, 選定30組以上的樣本數, 透過量測最後把這個樣本平均值所需要遵守的常態分佈決定出來. 一個常態分佈由樣本的平均值跟標準差決定, 所以算出這些樣本平均值的平均值跟標準差, 這個常態分佈就決定了.

根據常態分佈的密度函數, 

以平均值為中心, 樣本落於  +/- 1 sigma (正負一個標準差) 內的機率為: 0.6827 (68.27%)

以平均值為中心, 樣本落於  +/- 2 sigma (正負兩個標準差) 內的機率為: 0.9545 (95.45%), 一般製程工程師喜歡用這個來制定 SPC

以平均值為中心, 樣本落於  +/- 3 sigma (正負三個標準差) 內的機率為: 0.9973 (99.73%)

這邊可以看到, +/- 2 sigma 下, 你量測的樣本平均值是有 95.45% 的機率落在裡面的. 所以正常應該幾乎不太可能跑出去 (只有 4.55% 的很小機率), 現在假如你做 SPC 控制的量測, 結果發現量測值跑出去了, 那可能代表製程已經失控了, 已經不是原來的那個常態分佈.

這時候只要再以相同製程重新製作監控樣本後再量測一次, 如果量測數值還是跑出 +/- 2 sigma 外, 那就更肯定製程是出了問題, 因為這樣的機率只有 4.55% x 4.55% = 0.207% 非常微小的機率. 然後就是機台開始被停機, 設備跟製程開始 Trouble Shooting, 找原因囉!

 

這裡我們要建構這個常態分佈的函數, 第一個是可以算出機率密度, 會給你一個漂漂亮亮的鐘型高斯常態分佈曲線, 第二個是可以給出常態分佈的機率, 就是給定範圍然後會給你這個範圍下, 這個常態分佈的機率值. (積分面積)

 

I. 常態分佈機率密度函數

就是高斯分佈囉.

這個有精確的機率密度函數分布公式囉, 有分兩個.

(1) 第一個是比較簡單的, 平均值mean =0, 標準差sigma=1 的標準化常態機率密度函數.

精確公式如下

gauss formula Z.png

這個正規化的機率密度函數, 不需要任何參數的, 分佈的圖形如下

normal Z.png

 

(2) 第二個是廣義化的, 任意平均值 mean, 及標準差 sigma 下的常態機率密度函數!

精確公式如下

gauss formula.png

這個就需要兩個參數囉, 平均值 mean, 跟標準差 sigma. 不同的值標定不同的常態分佈.

分佈的圖形如下

normal.png

 

FORTH 程式

這兩個高斯分佈函數, 因為有公式囉, 所以很簡單的直接套一下就寫出來囉

(1) 平均值mean =0, 標準差sigma=1 的標準化常態機率密度函數: NormDistZ

: NormDistZ ( F: z -- F: f)
   fdup f*  -2.0E0  f/ fexp   
   [ 2.0E0 pi f* fsqrt ] fliteral f/
;

 

(2) 廣義化的, 任意平均值 mean, 及標準差 sigma 下的常態機率密度函數: NormDist 

利用標準化常態機率密度函數: NormDistZ, 很容易就把它廣義化, 推廣到任意的標準差跟平均值下的分佈.

因為 NormDist(x, mean, sigma) = (1/sigma)*NormDistZ((x-mean)/sigma )

 

: NormDist ( F: x mean sigma -- F: f)
    frot frot f-   ( F: sigma x-mean)
    fover f/       ( F: sigma {x-mean}/sigma)
    NormDistZ
    fswap f/
;

 

II. 常態分佈機率函數

再來這個就是重點囉, 常態分佈的機率函數. (Cumulative distribution function, 或簡稱 CDF)

機率密度函數, 給定範圍對它積分我們就可以得到這個給定範圍內的 "機率"! 為了方便, 這裡的範圍一律從負無窮大積分到特定x值. 如下圖所示,

norm Prob.png

這個積分, 可以真的用數值積分去做啦. 可是, 可能太慢的關係, 常看到數學家們會先引進一個幾乎一樣的 erf(x) 特殊積分函數, 再用多項式展開對它求值.

這裡, 我們要採用一個被稱之為 ACM Algorithm #209 的演算法, 來做這個計算. 這個演算法裡面的核心也是多項式展開, 就是用兩個多項式函數來對 z=0..1, 跟 z=1..3 之間分開做計算!

當 z=0..1,

f(z) = (((((((0.000124818987 * z

        – 0.001075204047) * z + 0.005198775019) * z
        – 0.019198292004) * z + 0.059054035642) * z
        – 0.151968751364) * z + 0.319152932694) * z
        – 0.531923007300) * z + 0.797884560593

當 z=1..3

f(z) = (((((((((((((-0.000045255659 * z
        + 0.000152529290) * z – 0.000019538132) * z
        – 0.000676904986) * z + 0.001390604284) * z
        – 0.000794620820) * z – 0.002034254874) * z
        + 0.006549791214) * z – 0.010557625006) * z
        + 0.011630447319) * z – 0.009279453341) * z
        + 0.005353579108) * z – 0.002141268741) * z
        + 0.000535310849) * z + 0.999936657524

演算法的細節, 大家就自己去翻閱文獻囉. 這裡就直接秀出程式實作的結果囉!

 

FORTH 程式實作

這兩個多項式剛好可以利用我們前一篇所說的的多項式函數定義來定義它們囉.

PolyDef Poly:0..1  ( F: z -- dist)
   0.000124818987E0  Coef,
  -0.001075204047E0  Coef,
   0.005198775019E0  Coef,
  -0.019198292004E0  Coef,
   0.059054035642E0  Coef,
  -0.151968751364E0  Coef,
   0.319152932694E0  Coef,
  -0.531923007300E0  Coef,
   0.797884560593E0  Coef,
PolyEnd

PolyDef Poly:1..3  ( F: z -- dist)
  -0.000045255659E0  Coef,
   0.000152529290E0  Coef,
  -0.000019538132E0  Coef,
  -0.000676904986E0  Coef,
   0.001390604284E0  Coef,
  -0.000794620820E0  Coef,
  -0.002034254874E0  Coef,
   0.006549791214E0  Coef,
  -0.010557625006E0  Coef,
   0.011630447319E0  Coef,
  -0.009279453341E0  Coef,
   0.005353579108E0  Coef,
  -0.002141268741E0  Coef,
   0.000535310849E0  Coef,
   0.999936657524E0  Coef,
PolyEnd

 

剩下就照文獻裡面演算法裡所描述的數學公式,  把他們組合起來就是囉.

: Dist:0..1  ( F: z -- F: dist)
    fdup fdup f* Poly:0..1   ( F: z dist)
    f* 2e0 f*
;

: Dist:1..3  ( F: z -- F: dist)
    2e0 f-
    Poly:1..3
;

 

最終的結果, 組合這兩個多項式函數成一個, 給定 z, 就會計算 負無窮大到 z 之間的, 正規化(mean=0, sigma=1)的常態分佈的機率給你.

: NormProbZ  ( F: Z — F: p)
    fdup f0>  Z大於零嗎? 先判斷, 旗號暫放

    fabs  2e0 f/

    fdup 3e0 f>    if   fdrop 1e0       else  |Z/2| 大於3
    fdup 1e0 f>    if   Dist:1..3       else  |Z/2| 介於1..3之間
                        Dist:0..1       then then  |Z/2| 介於0..1之間

    剛剛的, Z大於零嗎?

    if       1e0  f+                是就直接加1
    else     1e0  fswap f-     then 不是就1減
    2e0  f/
;

 

廣義的常態分布機率函數.

剛剛那個 NormProbZ 是正規化的, mean=0, sigma=1 的常態分佈機率函數.

現在一樣, 我們來推廣一下, 推廣到任意平均值 mean, 跟任意標準差 sigma 下都可以使用的機率函數. 這樣使用起來才會更方便.

很簡單的, 因為

NormProb(x, mean, sigma) = NormProbZ((x-mean)/sigma)

 

所以

: NormProb ( F: x mean sigma -- F: p)
    frot frot f- fswap f/
    NormProbZ
;

 

數值驗證

這幾個函數其實就是 Excel 裡的公式 NormDist(x,mean,sigma,CDF),  CDF是個旗號, FALSE 時為高斯分佈機率密度函數, TRUE 時為高斯分佈機率函數 CDF.

對應到我們所建造的函數如下

NormDistZ(x) = NormDist(x, 0, 1, FALSE)

NormDist(x, mean, sigma) = NormDist(x, mean, sigma, FALSE)

NormProbZ(x) = NormDist(x, 0, 1, TRUE)

NormProb(x, mean, sigma) = NormDist(x, mean, sigma, TRUE)

Excel 所跑出來四個函數的數據, (這裡廣義的函數, 一蓋以 mean=0.45, Sigma=0.78 來測試)

NormDist.png

 

我們寫個小程式, 列印一下四個函數所計算的結果

result1.png

result2.png

結果是一模一樣的! 驗證完畢, 這四個函數是可以正確使用的!

 

如何新增 FORTH 編譯器行為

在下一個範例前, 簡單介紹如何自己創造一個新的 FORTH 編譯器指令, 來新增 FORTH 編譯器的編譯行為.

在 ANS FORTH 標準裡面, 對於浮點數字的辨識輸入是靠 E這個字元加上一個數字.

例如要輸入  3.1416 這個浮點數字, 你必須輸入 3.1416E0 這樣的字串, 系統才會認知它是個浮點數字, 把它轉換成浮點數放上浮點堆疊.

為什麼會這麼不自然呢?? 原來啊, 小數點在最早的 FORTH 標準裡面, 已經被當作是雙整數的識別字元了. 例如, 假如你直接輸入 3.1416 這樣的數字, 系統會認為它是一個 31416 的一個雙整數, 而且會放入數字(參數)堆疊, 而不是浮點堆疊.

這樣, 某些情況下, 程式碼裡頭的浮點數可讀性就變差了哩! 一堆如 120E0  110E0 4.37E0 ... 看得頭昏眼花的, 科學記號的數字, 真的是好不自然喲!

假如能用個提示字元, 讓系統知道後面接過來的數字是浮點數, 那可讀性就會大大的提升. 例如 F 120  F 110  F 4.37 ... 用個提示字元 'F', 當系統看到 F 時, 就知道後面接著的數字是個浮點數, 就自動轉換它並放上浮點堆疊, 這樣會更方便. 就來實作這樣的 FORTH 編譯器行為的擴充吧!

這裡要介紹一下, FORTH 系統裡面有兩個司令官, 一個是翻譯器, 另一個是編譯器. 我們平常跟 FORTH 系統交談時, 是在翻譯器中. 一旦使用冒號指令時, 系統就會被切到編譯器裡, 接受文字輸入進行編譯的動作直到分號指令 ; 時才會離開. 同時系統裡也提供一個狀態變數 state, 當這個變數裡面為 true 的時候, 代表正在編譯的狀態. 為 false 的時候代表正在翻譯的狀態.

FORTH 的編譯器行為是非常簡單的, 就只是重複性的查字典把指令的位址找出來, 編到冒號定義裡面. 這樣的動作, 是非常不足以應付複雜的編譯行為的, 所以祖師爺 Chuck Moore 很聰明地留了一個後門, 每個FORTH指令有一個 bit 來標示這個指令是不是立即詞, 當這個 bit 為真, 就是 FORTH 術語裡面的立即詞. 所謂立即詞就是編譯器字, 編譯時永遠不會被編譯入冒號定義裡面,  而是作為編譯器編譯行為的一種擴充.

FORTH 編譯器每找到一個指令要編入記憶體前會先檢查這個 bit, 假如發現它是立即詞時, 這個指令並不會被編譯入記憶體, 反而是內建的編譯器會直接的將控制權交給這個指令立即執行. 透過這個指令來延伸額外的編譯器編譯動作. 所以透過一些特定的編譯器指令來觸發額外的編譯器行為, 這一直都是 FORTH 語言裡的拿手好戲. (例如: BEGIN, UNTIL, IF THEN... 這些都是編譯器字)

這個立即詞 bit 的設定很簡單, 只需透過 immediate 指令. 每當你定義完一個新指令時, 後面只要再接上 immediate 指令, 新指令的立即詞 bit 屬性就會被設成1, 就完成了立即詞的設定動作. 也就是你新增了一個新的不會被編譯, 反而是編譯時會被執行的編譯器指令.

這裡再來我們要利用一個很重要的, ANS FORTH 裡面字串轉換成浮點數的 >float 指令. 它的參數變化如下.

>float ( c-addr u -- true | false) ; ( F: -- r)

c-addr 是要轉換的字串起始位址跟長度u, >float 會直接嘗試將這個字串轉換成浮點數, 然後這個轉換是更有彈性的, 轉換過程中只要是遇到 'E', 'e', 'D', 'd', '+', '-' 這幾個符號, 後面的數字就視為 10 的次方數字. 假如轉換成功的話, 會在參數堆疊留下 true 的旗號跟轉換後的浮點數於浮點堆疊. 轉換失敗的話會留下一個 false 旗號, 浮點堆疊什麼都不留.

 

所以這個 F 的編譯器立即詞指令如下

: F ( <number> -- F: f)
  BL word count  從字串緩衝區裡以空白BL為分界取出下個指令的字串

  >float  嘗試轉換成浮點數
  invert abort" failed to convert to float point..." 檢查是否失敗,是的話就 abort
  
  state @  現在是在編譯狀態嗎?
  if postpone fliteral then  是的話就把剛轉換完的浮點數編入字典吧!
; immediate

這裡因為 fliteral 是個立即詞, 它會從浮點堆疊抓一個數字編譯入冒號定義中, 我們要借助它來編譯浮點數字, 但因為是立即詞, 無法被編譯. 所以必須利用 postpone 這個指令, 這個指令可以強迫一個立即詞被編譯入冒號定義中. 這樣就可以在冒號定義裡頭被執行囉.

(舊的標準是用 [compile] 這個指令, 新 ANS 標準改用 postpone 這個指令. 新標準的 postpone 比較聰明, 用在立即詞上面, 它等於舊標準的 [compile] 這個指令; 用在非立即詞上面, 它等於舊標準的 compile 這個延遲編譯指令.)

 

來測試一下吧

問, +/- 1 sigma 下, 常態分佈的信賴區間是多少啊??

F 1.0 NormProbZ   F -1.0 NormProbZ   f-  f.    結果是 0.682689 = 68.27%

再問, +/- 2 sigma 下, 常態分佈的信賴區間是多少呀??

F 2.0 NormProbZ   F -2.0 NormProbZ   f-  f.    結果是 0.9545 = 95.45% 

再問, +/- 3 sigma 下, 常態分佈的信賴區間是多少呀??

F 3.0 NormProbZ   F -3.0 NormProbZ   f-  f.    結果是 0.9973 = 99.73%

result3.png

 

 

最後, 來個更複雜的範例計算吧!

假定今天有台收獲蘋果果實, 並進行重量分類的機器. 這個機器啊, 可以將蘋果依據重量分類成 A1, A2, A3, ... A10 總共10個分類. 其規格依序為:

A1:150g - 160g,  A2:160g - 170g,   A3:170g - 180g, ... A10:240g - 250g

某次的收穫中, 透過這台機器, 得到如下的分佈:

A1: 0%, A2: 1%, A3: 1.5%, A4: 13.5%, A5: 21.5%, A6: 32.4%, A7: 18.9%, A8: 11%, A9: 0.2%, A10: 0%

這台機器啊, 因為分類的速度太快了, 最後被發現有個致命的缺點. 也就是說, 不是那麼精準. 例如: A3 分類內的, 理論上應該是只有 170g - 180g 之間重量的蘋果. 但是實際拿更精準的秤來測量後, 蘋果的重量會是以 175g為中心, 8.11g 的標準差的常態分佈散佈出去. 也就是說會有一些分佈跑出 170g - 180g.

現在有了一家蘋果汁工廠,  買了這批頻果, 

(1) 那請根據收穫時 Bin 分佈的比例, 估計出這批正確的重量分佈出來.

(2) 蘋果汁工廠決策單位想知道, 這批蘋果大於 210g 的數量比例到底有多少呢??

 

這樣的題目, 其實用 Excel 配上Excel內建的 NormDist(x,mean,sigma,CDF) 的統計函數, 很快的可以建模, 進行計算後得到答案.

但我們已經在我們的 FORTH 裡建立起這樣的函數囉, 所以也來計算一下囉.

計算的方式也很簡單啦, 因為真實的 Bin 為裡面, 除了自己的分佈外, 還會有其他 Bin 別所貢獻過來的分佈. 所以必須把每個 Bin 跑過來的比例給加進來, 這樣就是真實的分佈囉.

currentBin 這個變數指向目前想要計算的 Bin 位

以目前的這個 currentBin Bin 位,溢出去的分佈為何?

: CurrentBinDist ( F: x -- F: p-dist)
   currentBin @ BinCenterWeight@   SigmaInBin f@  NormProb
;

根據這個 currentBin Bin 位所溢出去的分佈, 計算出在位置 bin 所接收到溢出分佈的數量.

: CurrentRealRatioAtBin ( bin -- F: CurrentRealRatio@bin)
   BinLabel f@  
   fdup  BinRange f+   CurrentBinDist
              fswap    CurrentBinDist 
   f-       currentBin @ BinRatio@    f*
;

 

某特定 bin 位置, 將不同 currentBin 位置所溢出到這個特定 bin 位置的數目加一加, 就是最後這個特定 bin 的真實分佈囉.

: AllRealRatioAtBin ( bin -- F: AllRealRatio@bin)
   0E0
   TotalBin# 0
   do  i currentBin !        dup CurrentRealRatioAtBin
       f+
   loop  drop
;

 

剩下就是分佈畫一畫, 大於多少的比例給它全部加一加. 就是最終的答案囉!

result4.png

 

 

 

原始程式列表

 

 

\
\  Normal Cumulative Distribution
\            Frank Lin 2018.06.04
\

\
\ Polynomial by horner's method 
\

-414159 constant CheckFlag

: PolyDef
   create  CheckFlag  [ CheckFlag s>d d>f ] fliteral
           here      0 ,  ( num)
   does>   ( addr F: x — result)
      dup cell+
      swap @ 0     0e0          ( addr’ F: x b )
      do   fover f*             ( addr' F: x b*x)
           dup f@               ( addr’ F: x b*x Coef)
           f+                   ( addr’ F: x b*x+Coef)
           
           float+  ( step to next coef)
      loop
      drop fswap fdrop
;

: Coef,
   f,  1 over +!
;

: PolyEnd
   @ 1 < ( no data)
   swap CheckFlag <>                         ( Checkflag disappear)
   [ CheckFlag s>d d>f ] fliteral f= invert  ( Checkflag disappear on f stack)
   or or
   abort" Errror... Something went wrong..."
;

\
\ facility words
\

: F ( " -- F: f)
  BL word count >float
  invert abort" failed to convert to float point..."
  
  state @
  if postpone fliteral then
; immediate

\
\ Normal Distributions function and Cumulative distribution
\

PolyDef Poly:0..1  ( F: z -- dist)
   0.000124818987E0  Coef,
  -0.001075204047E0  Coef,
   0.005198775019E0  Coef,
  -0.019198292004E0  Coef,
   0.059054035642E0  Coef,
  -0.151968751364E0  Coef,
   0.319152932694E0  Coef,
  -0.531923007300E0  Coef,
   0.797884560593E0  Coef,
PolyEnd

PolyDef Poly:1..3  ( F: z -- dist)
  -0.000045255659E0  Coef,
   0.000152529290E0  Coef,
  -0.000019538132E0  Coef,
  -0.000676904986E0  Coef,
   0.001390604284E0  Coef,
  -0.000794620820E0  Coef,
  -0.002034254874E0  Coef,
   0.006549791214E0  Coef,
  -0.010557625006E0  Coef,
   0.011630447319E0  Coef,
  -0.009279453341E0  Coef,
   0.005353579108E0  Coef,
  -0.002141268741E0  Coef,
   0.000535310849E0  Coef,
   0.999936657524E0  Coef,
PolyEnd


: Dist:0..1  ( F: z -- F: dist)
    fdup fdup f* Poly:0..1   ( F: z dist)
    f* 2e0 f*
;

: Dist:1..3  ( F: z -- F: dist)
    2e0 f-

    Poly:1..3
;

: NormProbZ  ( F: Z — F: p)
    fdup f0>

    fabs  2e0 f/

    fdup 3e0 f>    if   fdrop 1e0       else
    fdup 1e0 f>    if   Dist:1..3       else
                        Dist:0..1       then then

    if       1e0  f+  
    else     1e0  fswap f-                  then
    2e0  f/
;

: NormProb ( F: x mean sigma -- F: p)
    frot frot f- fswap f/
    NormProbZ
;


false    \ if your system doesn't have pi const, set this to true
[if]

3.1415926535897932384626E0 fconstant pi

[then]


: NormDistZ ( F: z -- F: f)
   fdup f*  -2.0E0  f/ fexp   
   [ 2.0E0 pi f* fsqrt ] fliteral f/
;


: NormDist ( F: x mean sigma -- F: f)
    frot frot f-   ( F: sigma x-mean)
    fover f/       ( F: sigma {x-mean}/sigma)
    NormDistZ
    fswap f/
;

 

false
[if]


\
\ verification 
\

: s>f  s>d d>f ;

: .NormProbZ ( --)

   cr cr
   ." Z     NormProbZ"
   cr cr

   31 -30
   do  i s>f 10E0 f/
       fdup f.   4 spaces   NormProbZ f.  cr
   5
   +loop
;

0.45E0 fconstant  mean
0.78E0 fconstant sigma

: .NormProb ( --)

   cr cr
   ." x     NormProb (mean: 0.45 sigma: 0.78)"
   cr cr

   31 -30
   do  i s>f 10E0 f/
       fdup f.   4 spaces   mean sigma NormProb f.  cr
   5
   +loop
;


: .NormDistZ ( --)

   cr cr
   ." Z     NormDistZ"
   cr cr

   31 -30
   do  i s>f 10E0 f/
       fdup f.   4 spaces   NormDistZ f.  cr
   5
   +loop
;

: .NormDist ( --)

   cr cr
   ." x     NormDist (mean: 0.45 sigma: 0.78)"
   cr cr

   31 -30
   do  i s>f 10E0 f/
       fdup f.   4 spaces   mean sigma NormDist f.  cr
   5
   +loop
;

 


\
\ Example
\


\
\  Assuming there is one Apple fruit inspection sorting tool
\  it will sort apple fruits into 10 Bins as
\  A1 A2 A3 A4 ... A10, each bins with 10 grams weight scale
\            A1  A2  A3  A4  A5  A6  A7  A8  A9  A10
\  Weight:   150 160 170 180 190 200 210 220 230 240-250
\  Because the limitation of measurement, some error always will be introduced
\  during weighting. 
\  Even in the same bin, the real weight will have about sigma = 8.11 gram nature
\  spread.
\
\  In one apple fruit harvest, the distribution after sorting are
\   A1  A2  A3    A4    A5    A6     A7     A8   A9    A10
\   0%  1%  1.5%  13.5% 21.5% 32.4%  18.9%  11%  0.2%  0%
\
\  Now, one customer wants to order the apples with their weight over 210 gram,
\  Estimate the real distribution and the ratio of the Apple over 210 gram
\

 

: Data
   create 
   does> ( index addr -- addr')
      swap floats +
;

 

F 10 fconstant BinRange
  10  constant TotalBin#

Data BinRatio ( index -- addr)
\ A1     A2      A3        A4        A5        A6        A7        A8       A9     A10
F 0 f, F 1 f, F 1.5 f, F 13.5 f, F 21.5 f, F 32.4 f, F 18.9 f, F 11.0 f, F 0.2 f, F 0 f,

: BinRatio@  ( index -- F: ratio%)
   BinRatio f@    100E0 f/
;

Data BinLabel  ( index -- addr)
\  A1       A2       A3       A4       A5       A6       A7       A8       A9      A10
F 150 f, F 160 f, F 170 f, F 180 f, F 190 f, F 200 f, F 210 f, F 220 f, F 230 f, F 240 f,
    
: BinCenterWeight@ ( index -- F: gram)
   BinLabel f@  BinRange  2E0 f/  f+
;
    

 variable currentBin
fvariable SigmaInBin    F 8.11 SigmaInBin f!


: CurrentBinDist ( F: x -- F: p-dist)
   currentBin @ BinCenterWeight@   SigmaInBin f@  NormProb
;

: CurrentRealRatioAtBin ( bin -- F: CurrentRealRatio@bin)
   BinLabel f@  
   fdup  BinRange f+   CurrentBinDist
              fswap    CurrentBinDist 
   f-       currentBin @ BinRatio@    f*
;

: AllRealRatioAtBin ( bin -- F: AllRealRatio@bin)
   0E0
   TotalBin# 0
   do  i currentBin !        dup CurrentRealRatioAtBin
       f+
   loop  drop
;

: f.%  ( F: f --)
   F 100. f* f.  ." %"
;

: .dist
    cr
    F 0.  ( for summation)
    TotalBin# 0
    do  ." A" i 1+ .  ." (" i BinLabel f@  fdup f.  ." - "  BinRange f+ f.  ." ) :  "
        i AllRealRatioAtBin  fdup f.%  cr
        f+
    loop
    cr
    ." Total:"  f.%
    cr
;

: RatioOverBin  ( bin -- F: summation)
    F 0.
    TotalBin# swap
    ?do  i AllRealRatioAtBin  f+
    loop
;

: .Answer
     cr
    .dist cr
    ." The Ratio of the Apples with their weight over 210 gram (A7) is" cr
    6 RatioOverBin f.%  cr
;

6 set-precision       


[then]

 

 

xxx

 

 

 

 

arrow
arrow

    ohiyooo2 發表在 痞客邦 留言(0) 人氣()