概率與統計分析學習筆記

注：該文是上了開智學堂數據科學基礎班的課后做的筆記，主講人是肖凱老師。

概率與統計分析

描述性分析

用一個數字描述一組數字的特征。用一個數字來歸納一組數字，這個數字稱為統計量或統計指標。

• 均值、中位數：描述一組數據的集中趨勢

• 方差、標准差、四分位距：描述一組數據的離散趨勢

• 相關系數：上面兩大類指標都是對一個變量或一組數據的特征描述，如果要描述兩個變量或兩組數據的相關性，可以使用相關系數

from scipy import stats
from scipy import optimize
import numpy as np
import random
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns # 用於可視化

sns.set(style="whitegrid")

x = np.array([3.5, 1.1, 3.2, 2.8, 6.7, 4.4, 0.9, 2.2]) # 新建一個向量
print np.mean(x) # 均值
print np.median(x) # 中位數
print x.min(), x.max() # 最小最大值
print x.var() # 方差
print x.std() # 標准差
print x.var(ddof=1) # 樣本方差
print x.std(ddof=1) # 樣本標准差

3.1
3.0
0.9 6.7
3.07
1.75214154679
3.50857142857
1.87311810321

random.seed(123456789) # 種子不同，產生的隨機數序列也不同

random.random() # 生成 0 到 1 之間的隨機數，random 是隨機數模塊

0.6414006161858726

random.randint(0, 10)  # 生成 0 到 10 之間的一個整數

5

np.random.seed(123456789)

np.random.rand() # 生成 0 到 1 之間的數

0.532833024789759

np.random.randn() # 服從 0，1 的標准正態分布

0.8768342101492541

np.random.rand(5) # 同時生成 5 個從 0 到 1 均勻分布的隨機數

array([ 0.71356403,  0.25699895,  0.75269361,  0.88387918,  0.15489908])

np.random.randn(2, 4) # 生成 2 行 4 列從標准正態分布產生的數

array([[ 3.13325952,  1.15727052,  1.37591514,  0.94302846],
       [ 0.8478706 ,  0.52969142, -0.56940469,  0.83180456]])

np.random.randint(10, size=10) # 生成從 0 到 9 的整數，生成 10 個

array([0, 3, 8, 3, 9, 0, 6, 9, 2, 7])

np.random.randint(low=10, high=20, size=(2, 10)) # 生成從 10 到 19 的整數，2 行 10 列

array([[12, 18, 18, 17, 14, 12, 14, 10, 16, 19],
       [15, 13, 15, 18, 11, 17, 17, 10, 13, 17]])

fig, axes = plt.subplots(1, 3, figsize=(12, 3))

axes[0].hist(np.random.rand(10000)) # 從 0 到 1 均勻分布
axes[0].set_title("rand")
axes[1].hist(np.random.randn(10000)) # 標准正態分布
axes[1].set_title("randn")
axes[2].hist(np.random.randint(low=1, high=10, size=10000), bins=9, align='left') # 從 1 到 9 的整數的均勻分布
axes[2].set_title("randint(low=1, high=10)")

fig.tight_layout()

除了從已知分布中抽取隨機數，還可以使用 choice 做一些有放回或無放回的抽樣。

np.random.choice(10, 5, replace=False) # 從 0 到 9 之間抽 5 個，無放回。即做抽樣。replace = False 表示無放回抽樣

array([9, 0, 5, 8, 1])

上面介紹的隨機數種子都是全局種子，還有一種局部種子，如下例所示。

prng = np.random.RandomState(123456789) # 定義局部種子
prng.rand(2, 4)

array([[ 0.53283302,  0.5341366 ,  0.50955304,  0.71356403],
       [ 0.25699895,  0.75269361,  0.88387918,  0.15489908]])

prng.chisquare(1, size=(2, 2)) # 卡方分布

array([[  1.00418922e+00,   1.26859720e+00],
       [  2.02731988e+00,   2.52605129e-05]])

prng.standard_t(1, size=(2, 3)) # t 分布

array([[ 0.59734384, -1.27669959,  0.09724793],
       [ 0.22451466,  0.39697518, -0.19469463]])

prng.f(5, 2, size=(2, 4)) # F 分布

array([[  0.77372119,   0.1213796 ,   1.64779052,   1.21399831],
       [  0.45471421,  17.64891848,   1.48620557,   2.55433261]])

prng.binomial(10, 0.5, size=10) # 二項分布

array([8, 3, 4, 2, 4, 5, 4, 4, 7, 5])

prng.poisson(5, size=10) # 泊松分布

array([7, 1, 3, 4, 6, 4, 9, 7, 3, 6])

思考如何生成 100 個服從同一概率分布的隨機數？

如果是要從標准分布中生成隨機數，如 random.random 可生成 0 到 1 之間的服從均勻分布的隨機數，random.randint(0, 10) 可生成 0 到 10 的均勻分布的整數，np.random.randn() 可生成服從標准正態分布的隨機數。

如果要從非標准分布中生成隨機數，在建模時，這種情況經常遇到，譬如研究者能提出一個新的噪聲過程或現有分布的組合。用計算方法來解決復雜采樣問題，一般要依賴現有的簡單分布采樣方法。

• 離散變量的逆轉換采樣

• 連續變量的逆轉換采樣

• 拒絕采樣

具體可參考Computational Statistics with Matlab。

概率分布

一個隨機變量可以通過某個分布函數來描述特征。

連續分布：

• 正態分布、均勻分布

• 概率分布函數 Probability distribution function

離散分布：

• 泊松分布、離散均勻分布

• 概率質量函數 Probability mass function

np.random.seed(123456789) # 定義個種子

X = stats.norm(1, 0.5) # 定義一個正態分布，期望是 1，標准差是 0.5
print type(X) # 可以看出 X 其實是個分布
print X.mean() # 期望
print X.median() # 中位數
print X.std() # 標准差
print X.var() # 方差
print [X.moment(n) for n in range(5)] # 計算各階中心矩，下面會詳細介紹
print X.stats() # 打印統計量
print X.pdf([0, 1, 2]) # 概率密度函數 pdf，給定 x 值，給出密度函數的 y 值
print X.cdf([0, 1, 2]) # 累積概率函數，其實是 pdf 積分，也就是密度曲線下的面積，整個的面積為 1
X.rvs(10) # 用於生成隨機數，跟前面 random 模塊功能類似
X.interval(0.95) # 整個函數比較有用，計算曲線下面積為 0.95 時對應的兩個 x 值

<class 'scipy.stats._distn_infrastructure.rv_frozen'>
1.0
1.0
0.5
0.25
[1.0, 1.0, 1.25, 1.75, 2.6875]
(array(1.0), array(0.25))
[ 0.10798193  0.79788456  0.10798193]
[ 0.02275013  0.5         0.97724987]





(0.020018007729972975, 1.979981992270027)

穿插個知識點：上面的 moment 函數是計算隨機變量的各階中心矩，這里簡單介紹下統計中矩的概念。

定義 設 X 為隨機變量，c 為常數，k 為正整數。則量\(E[(X-c)^k]\)稱為 X 關於 c 點的 k 階矩。

有兩種比較重要的情況：

當 \(c=0\) 時，\(\alpha_{k}=E(X^k)\)稱為 X 的 k 階原點矩；

當 \(c=E(X)\) 時，\(\mu_{k}=E[(X-EX)^k]\) 稱為 X 的 k 階中心矩。

一階原點矩就是期望，一階中心矩為 0，二階中心矩就是 X 的方差。在統計學上，高於四階的極少使用，三、四階有些應用，但也不多。

矩的應用之一是用三階中心距 \(\mu_{3}\) 去衡量分布是否有偏。設 X 的概率密度函數為 f(x)，若 f(x) 關於某點 a 對稱，即

\[f(a+x)=f(a-x) \]

則 a 必等於 \(E(X)\)，且 \(\mu_{3}=E[X-E(X)]^{3}=0\)。如果\(\mu_{3}>0\)，則稱分布為正偏或右偏。如果\(\mu_{3}<0\)，則稱分布為負偏或左偏。特別地，對正態分布而言有 \(\mu_{3}=0\)，故如 \(\mu_{3}\) 顯著異於 0，則是分布與正態有較大偏離的標志。由於 \(\mu_{3}\) 的因次是 X 的因次的三次方，為抵消這一點，以 X 的標准差的三次方，即 \(\mu_{2}^{3/2}\)去除 \(\mu_{3}\)，其商

\[\beta_{1}=\frac{\mu_{3}}{\mu_{2}^{3/2}} \]

稱為 X 或其分布的“偏度系數”。

def plot_rv_distribution(X, axes=None):
    """Plot the PDF, CDF, SF and PPF of a given random variable"""
    if axes is None:
        fig, axes = plt.subplots(1, 2, figsize=(8, 3))
    
    x_min_999, x_max_999 = X.interval(0.999)
    x999 = np.linspace(x_min_999, x_max_999, 1000)

    x_min_95, x_max_95 = X.interval(0.95)
    x95 = np.linspace(x_min_95, x_max_95, 1000)

    if hasattr(X.dist, 'pdf'): # 判斷有沒有 pdf，即是不是連續分布
        axes[0].plot(x999, X.pdf(x999), label="PDF")
        axes[0].fill_between(x95, X.pdf(x95), alpha=0.25) # alpha 是透明度，alpha=0 表示 100% 透明，alpha=100 表示完全不透明
    else: # 離散分布
        x999_int = np.unique(x999.astype(int))
        axes[0].bar(x999_int, X.pmf(x999_int), label="PMF") # pmf 和 pdf 是類似的
    axes[1].plot(x999, X.cdf(x999), label="CDF")
    
    for ax in axes:
        ax.legend()
    
    return axes

fig, axes = plt.subplots(3, 2, figsize=(12, 9)) # 畫布，3 行 3 列

X = stats.norm() # 標准正態分布
plot_rv_distribution(X, axes=axes[0, :])
axes[0, 0].set_ylabel("Normal dist.")
X = stats.f(2, 50) # F 分布
plot_rv_distribution(X, axes=axes[1, :])
axes[1, 0].set_ylabel("F dist.")
X = stats.poisson(5) # 泊松分布
plot_rv_distribution(X, axes=axes[2, :])
axes[2, 0].set_ylabel("Poisson dist.")

fig.tight_layout()

概率密度函數（Probability density function）和累積概率分布函數（Cumulative distribution function）分別是什么？它們之間有什么樣的關系？

定義 設 X 為一隨機變量，則函數

\[P(X\le x)=F(x) (-\infty<x<\infty) \]

稱為 X 的分布函數。這里並未限定 X 為離散型的，它對任何隨機變量都有定義。

定義 設連續性隨機變量 X 有概率分布函數 F(x)，則 F(x) 的導數 \(f(x)=F'(x)\) 稱為 X 的概率密度函數。

CDF 是 PDF 的積分，也就是 PDF 曲線下的面積。

離散概率分布與連續概率分布的累積分布曲線看起來有什么不同？

離散概率分布的累積分布曲線是階梯狀的，連續概率分布的累積分布函數是平滑的。如上圖右邊的三幅圖都是累積分布曲線，上兩幅是連續概率分布的累積分布函數，最下面一幅是離散概率分布的累積分布曲線。

假設檢驗

假設檢驗在統計學里占有非常重要的地位，步驟如下：

• 構造原假設和備擇假設

• 選擇某個檢驗統計量，例如 T 檢驗統計量

• 收集數據

• 根據數據計算統計量和相應的 P 值

• 若 P 值很小，拒絕原假設

np.random.seed(123456789)
mu, sigma = 1.0, 0.5 # 均值，標准差
X = stats.norm(mu-0.2, sigma) # 均值為 0.8，標准差為 0.5
n = 100 
X_samples = X.rvs(n) # 生成 100 個隨機數
plt.hist(X_samples);

上面的采樣數據真正的參數是，均值 0.8，標准差 0.5。我們假裝不知道這些數據是由均值 0.8、標准差 0.5 的正態分布定義的，假裝只知道這 100 個數字，那么我們的問題是，它們的分布參數的均值是不是 1 呢？

如果樣本真的是由 mu=1 的分布產生的，那么樣本的均值應該離 mu 不太遠。而樣本的均值是可以通過 mean 函數計算出來的。

X_samples.mean() # 求樣本的均值

0.85830510224950851

樣本均值 0.85 離 1 不太遠，有可能樣本就是從均值為 1 的分布出來的。如果樣本均值是 0 或者 10，那么離 1 就太遠了，我們肯定就懷疑這樣本不是從均值為 1 的分布產生的。

首先定義樣本均值和總體均值的距離。

z = (X_samples.mean() - mu)/(sigma/np.sqrt(n)) # 樣本均值和總體參數均值的距離
print z

-2.83389795501

z 表示樣本均值和總體均值的距離。z 越大，表示樣本均值距離總體均值越大，那么就越要懷疑原假設，原假設即樣本是由均值為 mu 的分布產生的。z 越小，表示距離越近，則可以同意原假設，或者說不反對原假設。

那 z 多大算大，多小算小呢？這就需要一個臨界值，我們是通過抽樣分布的特點來判斷臨界值的。

這里 z 是服從標准正態分布的，所以下面可以方便計算 z 的概率。

在具體計算距離時，有兩種情況，一種是知道總體標准差的，如上式；另一種情況是不知道總體標准差，要用樣本標准差來代替，如下式。

t = (X_samples.mean() - mu)/(X_samples.std(ddof=1)/np.sqrt(n))
print t

-2.96803385457

stats.norm().ppf(0.025)

-1.9599639845400545

2 * stats.norm().cdf(-abs(z))

0.0045984013290753566

2 * stats.t(df=(n-1)).cdf(-abs(t))

0.0037586479674227209

以上是整個計算過程細節，如果你不需要了解這些細節，那可以直接用 stats 模塊中的單變量 T 檢驗 ttest_lsamp 來做這個事情。

t, p = stats.ttest_1samp(X_samples, mu) # 傳入證據及參數，數據即證據，注意這里並沒有傳入標准差 sigma，對於 T 檢驗，標准差是未知的
print t
print p

-2.96803385457
0.00375864796742

p 值小於 0.05，可以拒絕原假設，即認為數據不是從總體均值為 mu 的分布產生的。

fig, ax = plt.subplots(figsize=(8, 3))

sns.distplot(X_samples, ax=ax) # 畫真實數據的圖形，藍色線，分布其實是以 0.8 為期望值
x = np.linspace(*X.interval(0.999), num=100)
ax.plot(x, stats.norm(loc=mu, scale=sigma).pdf(x)) # 理論上值 mu=1

fig.tight_layout()

以上分析是單樣本假設檢驗，還有多樣本假設檢驗，如兩樣本，我們要看兩樣本是不是來自同一分布。

mu1, mu2 = np.random.rand(2) # 隨機產生兩個不同的參數
print mu1,mu2
X1 = stats.norm(mu1, sigma) # 標准差一樣
X1_sample = X1.rvs(n) # 生成第 1 組數據
X2 = stats.norm(mu2, sigma)
X2_sample = X2.rvs(n) # 生成第 2 組數據
t, p = stats.ttest_ind(X1_sample, X2_sample) # 兩樣本 T 檢驗
print t
print p # p 值太小，可以拒絕原假設

0.329631563852 0.692241009083
-5.48067917083
1.2796633829e-07

sns.distplot(X1_sample);
sns.distplot(X2_sample);

兩樣本 T 檢驗可用於判斷兩個生產線上的產品是否有顯著區別等。

如何理解零假設跟備擇假設？

零假設是一個在做實驗之前原有的假設，備擇假設是指，一旦你決定否定原假設，則這假設可備你選擇。

假設檢驗就是，根據觀察或實驗結果去檢驗零假設是否成立。

接受零假設意味着，從所獲數據來看，並無足夠的根據認為零假設不對，而不是說，從所獲數據證明了零假設是對的，因此，問題多少仍處於未決的局面。反之，否定原假設則意味着，按所獲數據有充足的理由認為零假設不對，即有充足的理由認為對立假設成立，故在一定限度內，可以說問題有了一個明確的結論。

如何理解 p 值，p 值代表了什么事件的概率？

p 值就是數據對原假設的支持程度。

p 值大於 0.05 時，應該拒絕備擇假設嗎？

p 值表達了實驗結果對零假設的支持程度，指定一個界限，若 p 小於所指定的界限，則決定不接受原假設，不然就接受原假設。這個界限最常用的是 0.05，其次是 0.01，0.1 等等。

p 值大於 0.05 時，要不要拒絕備擇假設，還是要看具體情況，不是所有的情況都非得把界限定為 0.05，或更小的 0.01。

當界限是 0.05 時，若 p 值大於 0.05，只能說明，並無足夠的理由認為零假設不對，而不是說備擇假設就是錯的，所以其實問題仍是未決的局面，還有待做更多的研究。

p 的具體值也很重要，如果 p=0.051，則盡管大於 0.05，但已很接近界限，因此可以說，雖然試驗結果仍維持了原假設，但已很不利於它，這時拒絕備擇假設會猶豫不決；如果 p=0.85，則可以說試驗結果不僅維持了原假設，還很有利於它，這時候拒絕備擇假設就沒什么疑問了。

一類錯誤和二類錯誤

在檢驗一個假設時，可能搞對了（原假設成立而接受它，或者原假設不對而否定它），也可能犯以下兩類錯誤之一：

• 原假設成立，但被否定了，這稱為第一類錯誤。

• 原假設不對，但被接受了，這稱為第二類錯誤。

是否犯錯誤，取決於所用的檢驗方法及所獲得的樣本即試驗數據。而犯哪一類錯誤，則取決於原假設是否成立，若原假設成立，則不會犯第二類錯誤；若原假設不成立，則不會犯第一類錯誤。

這么說不難理解，但難以把這兩類錯誤聯系起來，下面舉個很不嚴謹的例子來體會下。

我們經常說，法律不會冤枉一個好人，但也不會放過一個壞人。這里的冤枉好人，就是第一類錯誤。放過壞人，即不是好人，但被誤認為是好人，就是第二類錯誤。

在對嫌疑人一無所知時，應假定人是好人，即零假設是，人是好人。

「不冤枉好人」和「不放過壞人」之間有沒有沖突呢？

有。比如，假設某年某地出了個大冤案，一個好人被冤枉住了二十年牢才被發現是冤案，震驚全國，於是當年全國很多地方在判案時就更加嚴格，沒有十足的把握，沒有實打實的證據，就認為嫌疑人是好人放走，這樣大大減少了冤枉好人的概率，但同時也增加了放走壞人的概率，好多壞人因為證據不足而釋放。再比如，某些年份犯罪橫行，於是當局在全國組織了嚴厲打擊犯罪的行為，只要有犯罪證據就抓起來，最終確實處理了一批壞人，但也有不少好人被冤枉。

可見，要想「不冤枉一個好人」，就很有可能放過不少壞人，要想「不放過一個壞人」，就很有可能冤枉不少好人。就看需要證據的多少，也就是個界限，需要多少證據才能斷定嫌疑人不是好人，在假設檢驗中也就是檢驗水平 \(\alpha\)。

這里好多人會有困惑，干嘛非得在兩種錯誤「冤枉好人」和「放過壞人」之間游移，難道就不能全部正確斷案？

很難。雖然嫌疑人有沒有干壞事是確定無疑的，但我們並不知道這個事實，我們只有看證據，看能搜集到多少證據，看證據的力度。在假設檢驗里，就是看數據，你說零假設是對的，但在零假設的情況下，這數據就基本不會出現，現在出現了，你說零假設是正確的，誰信？這里的數據是采樣數據，不是總體，所以才需要假設檢驗。如果有全部證據，那判案就沒什么難的了。

要想保證判案正確率，就要搜集更多證據，要花更多時間和精力。假設檢驗只是在現有數據的情況下，給出的檢驗結果。可見，統計方法雖是一種有力的科學方法，但無法搞無米之炊。

非參數推斷

對一組數據，我們希望用一個分布來描述它的特征，比如正態分布，可以用一個公式來描述數據，但公式要能夠跟數據擬合，才能做參數推斷。但現實情況下，數據不能跟公式很好地擬合，這就需要讓數據自己說話，這就是非參數推斷。

np.random.seed(0) # 先定義個種子
X = stats.chi2(df=5) # 卡方分布
X_samples = X.rvs(100) # 用卡方分布生成 100 個隨機數
print X_samples

[ 11.87740141   5.62133928  12.48424465   2.04604944   4.03651267
   5.65798023   6.98206807   4.70149541   5.77577501   5.39173922
  10.38780279   3.75698915   0.32533609   6.55639806  15.03828046
   1.30710924  10.5932091   10.25799631   2.20999289   0.72662779
   9.05799941   8.92433257   1.9214374    1.35311843   2.9994597
   3.16738619   1.10704806   3.73673066   2.19559054   5.57505371
   4.25089547   5.72057502   2.72133868   3.35077109   2.35324329
   0.97911165   1.08783953   5.84349338   6.85360346   4.72424779
   5.62899185   2.61378879   3.47933568   4.50079157   1.72867406
   7.56307388   5.85384885   1.20166469   8.81234066   3.8249378
   3.25148226   9.02029235   5.46938978   6.76401061   4.71779307
   5.62775727   1.56601507   7.85948981   3.2261648    2.48433029
  12.48425043   7.58536446   3.59126794   7.1509955    6.40418286
   7.65468192   5.53864837   1.83570203   5.27195502   9.52840427
   3.17436839  12.37551801   2.4389755    6.12259088   2.71804987
   6.64531598   5.60681827   1.84624872   4.84275986   6.48454567
   2.16327754   8.52369442   4.13552411  11.58741286   4.04991086
   2.65887737   1.75703528   3.1679936    3.02630659  12.8572994
   1.63427047   7.32466393   2.00537696   1.19102088   7.64888505
   5.34099642   1.94603823   6.66554219   2.37237427   2.60384033]

假設我們不知道這些數是通過卡方分布產生的，現在要推斷分布的形狀，只能采取一些經驗的方法，其中一種方法是核分布。

kde = stats.kde.gaussian_kde(X_samples) # 高斯核函數，kde 是核函數密度估計，核分布是種非參數估計方法
kde_low_bw = stats.kde.gaussian_kde(X_samples, bw_method=0.25) # bw_method 為窗寬參數，該值越小，就越受到數據本身的影響
x = np.linspace(0, 20, 100)

fig, axes = plt.subplots(1, 3, figsize=(12, 3))

axes[0].hist(X_samples, normed=True, alpha=0.5, bins=25)
axes[1].plot(x, kde(x), label="KDE")
axes[1].plot(x, kde_low_bw(x), label="KDE (low bw)")
axes[1].plot(x, X.pdf(x), label="True PDF")
axes[1].legend()
sns.distplot(X_samples, bins=25, ax=axes[2])

fig.tight_layout()

可以看到，理論分布 True PDF 是非常平滑的，因為是公式定義出來的，而核密度估計 KDE 是直接由數據擬合出來的，所以又稱為經驗分布。bw_method 是窗寬參數，該值越小，核密度估計就越受到數據的影響，上面的數據是有雙峰的，kde_low_bw 確實受數據影響比較厲害。

參數方法和非參數方法都是用來估計給定數據的分布，都是要從實際數據推斷背后分布的樣子，區別在於，參數方法背后是由數學公式定義的，而非參數方法則沒有。

用非參數方法估計出數據背后的經驗分布后，可以使用該經驗分布來進一步抽樣。

kde.resample(10) # 使用非參數估計的經驗分布來進一步抽樣

array([[ 1.75376869,  0.5812183 ,  8.19080268,  1.38539326,  7.56980335,
         1.16144715,  3.07747215,  5.69498716,  1.25685068,  9.55169736]])

這里經驗分布 KDE 沒有累積分布函數，我們可以手工計算。

def _kde_cdf(x):
    return kde.integrate_box_1d(-np.inf, x) # 得到經驗分布的累積分布圖形

kde_cdf = np.vectorize(_kde_cdf) # 向量化函數

fig, ax = plt.subplots(1, 1, figsize=(8, 3))

sns.distplot(X_samples, bins=25, ax=ax)
x = np.linspace(0, 20, 100)
ax.plot(x, kde_cdf(x))

fig.tight_layout()

思考何種情況下使用非參數推斷？

對一組數據，我們希望用一個分布來描述它的特征，比如正態分布，可以用一個公式來描述數據，但公式要能夠跟數據擬合，才能做參數推斷。但現實情況下，數據不能跟公式很好地擬合，這就需要讓數據自己說話，這就是非參數推斷。

你用過 Python 的內置函數 map 嗎？Numpy 里面有類似的方法嗎？

map 可以把一個函數作用到向量的每個元素，np.vectorize 也可實現這種功能，把函數向量化，函數就可以處理向量了。

補充閱讀

• scipy.stats 官方文檔
• Scipy lectures notes 第 15 章
• Numerical Python 第 13 章
• 統計學的世界 第 3、4 部分
• Think Stats 1-9 章

參考資料

• 統計學漫話
• 概率論與數理統計