Lesson 5: 機率分配與圖形（Probability distributions and the Shapes）

(Objective:

學習繪製離散型與連續型隨機變數的機率密度函數圖（PMF, PDF）與累積機密度函數圖（CDF），譬如，Binomial, Poisson, Normal, T, $\chi^2, \beta$ … 等分配.
熟悉繪圖技巧與分配函數的觀念，以便繪製出能代表分配特色的圖形。
透過繪製 PMF, PDF 與 CDF 圖，了解並熟悉各種分配函數的「長相」與參數間的關係。

Prerequisite:

熟悉前面單元的繪圖技巧。
熟悉 Python 基本的資料型態，特別是矩陣資料。

範例 1：Normal Distribution

繪製常態分配 $N(\mu, \sigma^2)$ 的機率密度函數圖 PDF。其中 $N(\mu, \sigma^2)$ 自選。

注意事項：

統計相關的計算都涵蓋在 scipy 套件裡。
本範例也呈現了 95% 信賴區間的範圍並塗上顏色（patched with color）。
本範例刻意採用財經界著名的 538 圖型 plt.style.use(‘fivethirtyeight’)。看看與一般的圖型有何差別。另外 matplotlib.pyplot 還有其他的 style，如用於 R 語言的 plt.style.use(‘ggplot’) ，更多選擇請詳見手冊說明。另，當採用某種 style 之後，後續的程式都會沿用，如欲回復，可使用 plt.style.use(‘default’)
熟悉機率分配常用的兩種指令，以常態為例， norm.pdf and norm.ppf，分別為計算 PDF 與 inverse CDF。
學習繪製累積密度函數 CDF 圖。
繪圖區的格線（grid）採虛線。

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# print(norm.__doc__) # see basic information
mu, sigma = 0, 1
xlim = [mu - 5 * sigma, mu + 5 * sigma]
x = np.linspace(xlim[0], xlim[1], 1000)
y = norm.pdf(x, mu, sigma)

plt.plot(x, y)
plt.style.use('fivethirtyeight') # 538 style

# patch the area of confidence interval
conf = 0.95
ci = norm.ppf([(1-conf)/2, (1+conf)/2]) # inverse CDF　function
x_ci = np.linspace(ci[0], ci[1], 1000)
y_ci = norm.pdf(x_ci, mu, sigma)

plt.fill_between(x_ci, y_ci, 0, color = 'g', alpha = 0.5)
plt.grid(True, linestyle='--', which='major')
plt.title('Normal({}, {})'.format(mu, sigma))
plt.show()

練習：Normal Distribution 兼 Python Broadcasting 練習

學習機率分配的參數與圖型分布的關係，最好的方式便是寫支程式，一次呈現多個分配圖，方能看出參數對形狀的影響。請繪製常態分配 $N(\mu, \sigma^2)$ 的 PDF 圖，其中 $\mu = 0$ ，而 $\sigma = 1, 2, 3, 4, 5$ , 即繪製下左圖。

注意事項：

想對每個分配函數圖有清晰的觀察視野，必須謹慎的選擇 x 值的範圍。譬如，下左圖考慮了當 $\sigma = 5$ 時，需要比較大的範圍，才能呈現該有的鐘形。
想要畫條數條函數線時（有規律的參數），第一個想法通常是用迴圈，此時最好先畫一條函數線成功後，再加上迴圈，免得太多不熟悉的東西加在一起，一旦出錯，不容易檢查出來。
接著，改換參數為固定 $\sigma =1$ ，即 $N(\mu, 1)$ ，改變 $\mu = 1, 2, 3, 4, 5$ 。
只會用迴圈，不會 broadcasting，便浪費了 Python 的優勢。拿掉迴圈，嘗試 broadcasting 怎麼完成任務。下列程式碼提供參考。關鍵在指令 Y = norm.pdf(x.reshape(-1, 1), mu, s) 將 x 豎起來成為一行向量，與 s 的列向量相乘形成一矩陣。接著 plt.plot 同時對 Y 矩陣的每一行繪製一條線，在此以大寫變數名稱 Y 區分矩陣與向量之別。另外，加入 legend 的做法也一併示範，這個做法稱為 list comprehension (loop in bracket)。
練習使用分割圖 subplot，將前述兩種情況放在一起比較，如下右圖。
本範例展示 norm.pdf()，而累積機率分配函數則是 norm.cdf()。其他分配也比照這個命名原則。

# Broadcasting 示範
fig, ax = plt.subplots()
xlim = [-15, 15]
mu = 0
s = np.arange(1, 6)
x = np.linspace(xlim[0], xlim[1], 1000)
Y = norm.pdf(x.reshape(-1, 1), mu, s)
# print(y.shape)
label = ["$\sigma$={}".format(i) for i in s]
plt.plot(x, y, label=label)
plt.legend()
plt.show()

範例 2： $\chi^2$ Distribution

繪製 $\chi^2$ 分配的 PDF 圖，其中自由度 $\nu = 4:2:32$ , 如下圖。

注意事項：

scipy 的 $\chi^2$ PDF 函數指令為 scipy.stats.chi2.pdf(x, $\nu$ ).
繪製 $\chi^2$ PDF 函數可以使用時間暫留製作動畫效果（plt.pause(0.5)），以展現隨著自由度增加，圖形逐漸像右移動，且右傾現象逐漸趨緩。當然，程式必須必須以 py 為副檔名。

import numpy as np
from scipy.stats import chi2
import matplotlib.pyplot as plt

xlim = [0, 50]
x = np.linspace(xlim[0], xlim[1], 1000)

# df 
df = np.arange(4, 32, 2)
# fix xlim before animation
plt.figure()
plt.axis([xlim[0], xlim[1], 0, 0.2])
for i in df:
    y=chi2.pdf(x, i)
    plt.plot(x,y, lw=1, color='blue')
    # plt.pause(0.5)

plt.title(r'$\chi^2$ Distribution')
plt.yticks([0, 0.1, 0.2])
plt.show()

練習： $t$ Distribution

繪製 $t$ 分配的 PDF 圖，其中自由度 $\nu$ 從 0.1:0.1:1 延續到 3:3:30, 如下圖。

注意事項：

t 分配為人熟知的特點是，隨自由度升高，漸近於標準常態。在此練習裡，先畫一條較粗、明顯的標準常態 PDF 圖，再觀察 t 分配隨著自由度改變的變化（可以加入動畫的表現）。
需要組合向量 [0.1 0.2 … 0.9 1] 與 [3 6 … 30] ，方便 “for loop” 進行.
scipy 的 t 分配 pdf 為 scipy.stats.t.pdf(x, $\nu$ ).

練習： $\beta$ Distribution

繪製 $\beta(a, b)$ 分配的 PDF 圖，其中參數 $a=9, 1 \leq b \leq 30$ , 如下圖。

注意事項：

scipy 套件的 $\beta(a, b)$ PDF 函數為scipy.stats.beta.pdf(x, a, b)。

習題： $\beta$ Distribution

$\beta(a, b)$ 分配有兩個參數影響分佈的形狀，試著嘗試各種可能的組合，看看能得到多少種「樣子」。譬如夏曆四張圖。

注意事項：

下圖採 subplots(2,2)，並分別對每個分割圖做迴圈繪圖。
下圖的採用 super title: plt.suptitle()

習題： $F$ Distribution

$F$ 分配同樣有兩個參數左右分佈的形狀。如前的練習，嘗試各種可能的組合，盡可能地挖掘出 $F(n_1, n_2)$ 分配可能的樣子、兩個參數分佈的影響是甚麼？

注意事項：

$F$ 呈現右偏的趨勢是否會因為參數的選擇而變成不偏或甚至左偏？值得寫成是好好探索。
繪多條函數線圖時，顏色的選擇也常扮演關鍵角色。下圖的顏色碼是 #3399FF，取自網站：RGB Color Codes Chart

範例 2: Binomial Distribution

繪製二項分配 $B(n, p)$ 的 PMF(probability mass function）圖，其中參數 $n = 20, p =0.7$ 。

注意事項：

scipy 套件的二項分配 pmf 是 scipy.stats.binom.pmf。
離散型分配的 PMF 圖，可以採 STEM 或 BAR 圖。下圖採用 STEM 的方式，讀者可以試著用 BAR 長條圖。
在套件 matplotlib 裡的 STEM 圖，可以採指令 stem 或在 plt.plot() 的 attribute 選擇加入 drawstype=”steps-pre”，如下列程式碼所示。
下圖套用了 538 的 style: plt.style.use(‘fivethirtyeight’) 。
請留意，stem 指令有些 attributes 供選擇，變化很多，不妨多試試。

import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt

n, p = 20, 0.7
x = np.arange(n + 1)
y = binom.pmf(x, n, p)
plt.stem(x, y, linefmt='g-', markerfmt='o', basefmt = 'C1--')
# plt.plot(x, y, drawstyle = 'steps-pre')
plt.title('$B(20, 0.7)$')
plt.show()

範例 3: Binomial Distribution

本範例示範三種不同的表現方式來繪製 PMF 圖。另外加上一張 CDF 圖。用了四種繪製離散分配的方式。

注意事項：

下圖運用了 4 x 1 的分割圖，除了採用 super title 外，還運用分享 X 軸的方式（share x ticks and ticklabels column-wise），避免重複多次相同的 X 軸。
bar 指令可以定義 width, color, edgeclor, …等（甚至可以細到針對每一個長條）
階梯圖（Step plot）通常用來繪製離散分配的 CDF 圖。

import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt

n, p = 20, 0.7
x = np.arange(n + 1)
y = binom.pmf(x, n, p)
fig, ax = plt.subplots(4,1, sharex = 'col', figsize = [6, 9])
ax[0].stem(x, y, linefmt='k-', markerfmt='ko', basefmt = 'C2--')
ax[1].bar(x, y, width = 0.9, color = 'r', edgecolor = 'y' )
ax[2].vlines(x, 0, y, lw = 5, alpha = 0.5)
Y = binom.cdf(x, n, p)
ax[3].plot(x, Y, drawstyle = 'steps-pre')
plt.suptitle('The PMF and CDF of Bino({}, {})'.format(n, p))
plt.show()