汪群超 Chun-Chao Wang

Dept. of Statistics, National Taipei University, Taiwan

Lesson 6: 隨機亂數與相關圖形(Random numbers and the associated graphs)

Objective:

  1. 學習自幾個著名的離散型與連續型分配產生隨機亂數。
  2. 學習不同的亂數產生方式。
  3. 隨機亂數的視覺化表現,譬如直方圖(Histogram)、盒鬚圖(Boxplot)、機率圖(Probability Plot)與經驗 CDF 圖(Empirical CDF)。

Prerequisite:

  1. 懂得機率分配的意義與分配圖形的繪製(PDF, CDF)。

範例 1:直方圖 Histogram

從常態分配 N(\mu, \sigma^2) 抽取隨機亂數並繪製直方圖。

注意事項:

  1. 本範例呈現兩種生成亂數的方式: numpy method 及 scipy function。

  2. 了解 “seed number” 扮演的角色。

  3. 直方圖的做法有二: by counts or by frequency(如圖).

  4. 繪製直方圖要特別留意樣本數造成的影響。另一個影響因素是 bin numbers(組距數)。針對相同樣本數(特別是中等數量,譬如 100, 200),試試不同的 bin 數量所造成的視覺變化。

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

# random number generator with random seed
rng = np.random.default_rng() # (seed=...)
x = rng.normal(loc = 0, scale = 1, size = 1000)
# x = norm.rvs(loc = 0, scale = 1, size = 1000)

fig, (ax1, ax2) = plt.subplots(2, 1)
# typical histogram by counts
ax1.hist(x, bins = 10, alpha = 0.6, color = 'b', edgecolor = 'y', linewidth = 1)
# histogram by frequency 
ax2.hist(x, bins = 10, alpha = 0.6, color = 'b', edgecolor = 'y', linewidth = 1, density = True, rwidth = 0.9)
plt.suptitle('Histogram by counts and frequency')
plt.show()


範例 2:盒鬚圖 Boxplot

繪製兩組不同資料的盒鬚圖, 分別來自 N(0, 1) and N(1, 2)

注意事項:

  1. 示範使用與不使用 random seed 來產生亂數之不同。當使用固定的 random seed 時,是否會產生一樣的亂數?試著印出來看看。有些使用者習慣使用固定的 random seed,卻不知每次都產生相同的亂數,造成後續的計算得到相同的結果。

  2. 繪製盒鬚圖的技巧非常豐富,本範例僅提供一二,譬如, 使用了 dictionary “boxprops” and “flierprops”

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

rng = np.random.default_rng(0) # seed is fixed
n = 1000 # sampe size
x1 = norm.rvs(size = n, random_state = rng) # standard normal
x2 = norm.rvs(loc = 1, scale = 2, size = n)
# run a few times and check the random numbers generated 
# with and without fixing the random seed.
# print(x1[0:3])
# print(x2[0:3])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (9, 6))
ax1.boxplot(np.c_[x1, x2], notch = True, vert = True, labels = ['X1','X2'])
# plt.show()

# use boxplot properties
boxprops = dict(linestyle = '--', linewidth = 3, color = 'darkgoldenrod')
flierprops = dict(marker='o', markerfacecolor = 'green', 
    markersize = 8, linestyle = 'none') # define outliers
labels = ['N(0,1)','N(1,1)']
ax2.boxplot(np.c_[x1, x2], boxprops = boxprops,
    flierprops = flierprops, labels = labels)
plt.show()


範例 3:常態機率圖 Normal Probability plot

各別繪製兩組資料的常態機率圖,其中一組資料來自常態分配,另一組來自右偏的 \chi^2 分配。

注意事項:

  1. 本範例程式刻意採用 numpy 的方式產生隨機樣本。

  2. 機率圖(probability plot)可以設定為各種分配。當設定為常態分配時,可以用來測試該樣本是否來自常態。

  3. 當機率分配設定為常態以外的分配時,必須額外指定 “sparams” 作為該分配的參數。下列程式碼註解的部分設定了卡方分配(自由度為 2)。

import numpy as np
import scipy.stats as stats

x1 = np.random.normal(loc = 2, scale = 5, size=100) 
x2 = np.random.chisquare(df = 2, size = 100)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (12, 4))
stats.probplot(x1, dist = "norm", plot = ax1)
stats.probplot(x2, dist = "norm", plot = ax2)
# stats.probplot(x2, dist = "chi2", sparams = 2, plot = ax2)
ax1.set_title('Data from Normal')
ax2.set_title('Data from $\chi^2$')
plt.suptitle('probability Plot')
plt.show()

範例 4:經驗累積密度函數 Empirical CDF (ECDF)

繪製一組來自常態分配 N(0,1) 樣本的 ECDF 圖,並比較 ECDF 圖與真實 CDF 圖之差異。

注意事項:

  1. ECDF 圖來自樣本數為 n 的樣本,當 n 愈大時,與真實 CDF 愈接近。因此 ECDF 為真實 CDF 的模擬版本。

  2. scipy 模組提供了 “cumfreq” (cumulative frequency over sorted sample) 做為 ECDF。 如下列程式碼之 Approach 1。

  3. 模擬 CDF 最簡單的方式就是將每個樣本點的發生機率設為一樣,即 \frac{1}{n}. 如下列程式碼之 Approach 2。

  4. 將 ECDF 與真實 CDF 畫在一起, 並觀察樣本數的影響。

import numpy as np
from scipy.stats import norm
from scipy.stats import cumfreq  # for ECDF

n = 100
sample = norm.rvs(size = n)
# Approach 1
num_bins = n # cumulative frequency over each sample
res = cumfreq(sample, num_bins)
# print(res.lowerlimit, res.binsize, res.cumcount.size)
x = res.lowerlimit + np.linspace(0,
    res.binsize * res.cumcount.size, res.cumcount.size)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize = (6, 9))
# ax1.bar(x, res.cumcount/n, width=res.binsize)
ax1.plot(x, res.cumcount/n, drawstyle = 'steps-pre')
ax1.set_title('The Empirical CDF with {} steps'.format(num_bins), fontsize = 'small')
ax1.grid(True)

# Approach 2
x = np.sort(sample)
Y = np.arange(1, n+1) / n
ax2.plot(x, Y, drawstyle = 'steps-pre', label = 'Empirical CDF')
# draw the theoretical CDF of Z
Y_ = norm.cdf(x)
ax2.plot(x, Y_, color = 'r', linestyle = '--',
    linewidth = 5, alpha = 0.5, label = 'real CDF')
ax2.legend(), plt.grid(True)
ax2.set_title('The Empirical CDF with {} samples'.format(num_bins), fontsize = 'small')
plt.show()


練習:抽樣分配 Sampling Distribution Y = X^2

假設隨機變數 X 服從常態分配 N(0, 1),令 Y = X^2, 則隨機變數 Y 服從 \chi^2 分配且自由度為 1。本範例要藉由隨機抽樣來證實這個知名的定理。做法如下:

  1. 從常態分配 N(0, 1) 隨機抽取 n 個樣本,並計算每個樣本的平方。

  2. 繪製直方圖(貼上真實的 \chi^2(1) 的 PDF 圖)如下圖左上。

  3. 繪製 Boxplot, Probability plot 與 ECDF(貼上真實的 CDF 圖),如下圖,作為佐證,尤其是 ECDF 圖與真實 CDF 之相似度,隨樣本數增加而更像。

  4. 可以將 ECDF 的計算寫成一個函數,如 def(sample): …. return x, y

    注意事項:

    1. 繪製直方圖時,必須留意組距(bins)數的多寡,這將影響觀察的判斷。

    2. 同樣值得注意的是 Y 軸的調整,以便得到較好的觀察視野。當貼上 PDF 線圖時,直方圖必須採用相對頻率。


練習:抽樣分配 Sampling Distribution

Y = X_1 + X_2, where X_1 \sim \chi^2(2), X_2 \sim \chi^2(4)

利用隨機樣本並繪製直方圖(含真實 PDF 圖)與 ECDF 圖(含真實 CDF 圖)來證實 Y \sim \chi^2(6). 。

注意事項:

  1. 分別選擇樣本數 n=30, 100, 500, 1000,並觀察樣本數多寡的影響。


練習:抽樣分配之中央極限定理 Sampling Distribution for Central Limit Theorem (CLT)

中央極限定理之簡易描述:

  • Let X_1, X_2, \cdots, X_n \sim F(\mu, \sigma)

  • Let Y = (X_1 + X_2 + \cdots + X_n) / n

  • As n \longrightarrow \infty, \;\;Y \sim N(\mu, \sigma/\sqrt{n})

模擬實驗設定:

  • 令隨機樣本 F 來自 unif(0,1), t(3), t(10), \chi^2(2), \beta(9,2)

  • 令樣本數為 n=5, 10, 30, 50, 100, 300, 500, 1000

  • 令每次實驗共抽樣 N = 10000 次,即得到 N 個平均數 Y_i, i=1,2,\cdots, N

實驗結果觀察:(for each F and n)

  1. 繪製 N 個平均數 Y_i 的直方圖,並與定理期望的常態分配 N(\mu, \sigma/\sqrt{n}) 的 PDF 線對照,譬如下圖左的直方圖來自 F=\chi^2(2), n=100, 及 N=10000

  2. 計算 NY_i 的平均數 \hat{\mu} 與樣本數 n 的關係如下圖右上,並畫上一條理論的 \mu 值(譬如 F=\chi^2(2)\mu=2)。 同樣地,也計算 Y_i 的樣本標準差 \hat{\sigma} 與樣本數的折線圖,如下圖右下。

  3. 繪製下圖右的兩張圖時,必須謹慎選擇 Y 軸的範圍與 X 軸的 xticklabels,方能呈現好看的圖。這是從事統計計算必須具備的能力。不只是把圖畫出來而已。


範例 5:多組樣本的盒鬚圖 Boxplot for multiple data sets

盒鬚圖經常被用來呈現多組樣本的敘述統計量之間的差異。繪製如下圖的 N 個盒鬚圖,其中的資料來自獨立的標準常態分配。每組資料的樣本數皆為 n

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

n, N = 100, 20
X = norm.rvs(size=(n, N))
plt.boxplot(X)
plt.title('The boxplots of {} independent random samples from $N(0,1)$'.format(n))
plt.show()


習題:

  1. 選擇三個曾經學過的抽樣分配,透過隨機抽樣的模擬方式,與理論互相驗證。其中樣本數、模擬次數與繪製何種圖形,請自行決定。

  2. 給定四個數字 (2, 4, 9, 12)。從這四個數字中隨機抽取四個數字(取後放回)並計算其平均數。假設隨機變數 Y 代表這四個數字的平均數。請繪製隨機變數 Y 的 PMF。本題可以直接計算每個平均數的機率,但在此請使用隨機抽樣的方式,估計出這些機率值。其中抽樣的次數可以高至百萬以上。
    Hint: 取後放回的抽樣方式,可以使用 numpyrandom.choice。而從重複的數字中取出不重複的部分,可以使用 numpyunique。(本題出自 Statistical Inference by G. Casella & R.L. Berger)

Histogram of averages of samples with replacement from the four numbers {2,4,9,12}

專題:

  1. \{x_i, i=1,\cdots, n\} 代表來自標準常態 N(0,1)n 個隨機樣本。統計量 G_1 表示為

    G_1 = \sqrt{\frac{n}{6}} \hat{s}

    其中 \hat{s} 為偏態係數(skewness)的估計值。請利用蒙地卡羅模擬(Monte Carlo Simulation)驗證統計量 G_1 服從標準常態 N(0,1)。其中蒙地卡羅模擬的環境設定(scenarios)為:

    • 樣本數 n = 10, 20, 30, 50, 100, 300, 500

    • 針對每個樣本數 n,模擬次數皆為 N=10000

    • 繪製 n = 10n = 500 時,統計量 G_1 的直方圖與 ECDF 圖。並分別畫上對應的標準常態 PDF 與 CDF 圖。

  2. 同上,但令統計量為

    G_2 = \sqrt{\frac{n}{24}} (\hat{k} - 3)

    其中 \hat{k} 為峰態係數(Kurtosis)的估計值。同樣利用蒙地卡羅模擬,驗證統計量 G_2 服從標準常態 N(0,1)。蒙地卡羅模擬的環境設定同上。

  3. 同上,但統計量為

    G_3 = G_1^2 + G_2^2 = \sqrt{\frac{n}{6}} \left(\hat{s}^2 +\frac{(\hat{k} - 3)^2}{4}\right)

    同樣利用上述的蒙地卡羅模擬,驗證統計量 G_3 服從卡方分配 \chi^2(2)

商學院  7F16
ccw@gm.ntpu.edu.tw
(02)8674-1111 
ext 66777

部落格統計

  • 107,550 點擊次數
%d bloggers like this: