Histgram直方圖適合於單個變數的value分佈圖形
seaborn在matplotlib基礎上做了更高層的抽象，方便對基礎的圖表繪製。也可以繼續使用matplotlib直接繪圖，但是呼叫seabon的set()方法就能獲得好看的樣式。
# Import plotting modules import matplotlib.pyplot as plt import seaborn as sns # Set default Seaborn style sns.set() # Plot histogram of versicolor petal lengths plt.hist(versicolor_petal_length,bins = 20) plt.xlabel('petal length (cm)') plt.ylabel('count') # Show histogram plt.show()
但是直方圖有以下缺點:
1. 受bins引數的影響，可能相同的資料圖形化後卻得到完全不同的解讀
2. 我們使用bin將資料做了切分，因此並不能展示全部資料
3. 直方圖只能用於單個變數的數值分佈，並不能很好地看出資料分佈和另一個類別變數之間的關係
為了克服這些缺點，我們可以使用Bee swarm(蜂群)圖表，也可以稱為swarm圖
import seabon as sns # Create bee swarm plot with Seaborn's default settings _ = sns.swarmplot(x='species',y='petal length (cm)',data=df) # Label the axes _ = plt.xlabel('species') _ = plt.ylabel('petal length (cm)') # Show the plot plt.show()
注意在swarm圖中水平軸的寬度沒有實質關係，因為數量很多的話會自然沿著類別型變數來做伸展，主要看資料密集度以及其y軸高度的分佈。
從上圖我們可以看出以下資訊:
不同花的類別對應的花瓣長度有明顯的不同，因此花瓣長度可以作為一個非常重要的變數。
swarm圖的缺點：
如果資料量過大的話，資料之間互相重疊，不容易看出單個x範圍內資料的分佈情況，這時我們可能要使用ECDF(empirical cumulative distribution function)
ECDF:
ECDF可以全景看到資料是如何分佈的，類似於統計中的概率累積分佈圖。主要用於給定一個一維陣列，我們們研究該資料的分佈情況
def ecdf(data):
"""Compute ECDF for a one-dimensional array of measurements."""
n = len(data)
# x-data for the ECDF: x
x = np.sort(data)
# y-data for the ECDF: y
y = np.arange(1,n+1)/n
return x, y
# Compute ECDF for versicolor data: x_vers, y_vers x_vers,y_vers = ecdf(versicolor_petal_length) # Generate plot plt.plot(x_vers,y_vers,marker='.',linestyle ='none') # Label the axes plt.xlabel('lenth') plt.ylabel('ECDF') # Display the plot plt.show()
為了進行對比，我們也可以將三種花的CDF進行堆疊，
# Compute ECDFs x_set,y_set = ecdf(setosa_petal_length) x_vers,y_vers = ecdf(versicolor_petal_length) x_virg,y_virg = ecdf(virginica_petal_length) # Plot all ECDFs on the same plot plt.plot(x_set,y_set,marker='.',linestyle='none') plt.plot(x_vers,y_vers,marker='.',linestyle='none') plt.plot(x_virg,y_virg,marker='.',linestyle='none') # Annotate the plot plt.legend(('setosa', 'versicolor', 'virginica'), loc='lower right') _ = plt.xlabel('petal length (cm)') _ = plt.ylabel('ECDF') # Display the plot plt.show()
產生如下的效果:
summary統計數值計算和圖表化
mean/median/percentile/IQR(25%-75%)
mean_length_vers = np.mean(versicolor_petal_length) meadian_length_vers = np.median(versicolor_petal_length) meadian_length_vers = np.percentile(versicolor_petal_length,[25,50,75])
其中median中位數就是50%的percentile
box plot用於圖形化顯示percentile資料的圖表。一般x軸為類別性變數，y軸為連續性數值變數
# Create box plot with Seaborn's default settings sns.boxplot(x='species',y='petal length (cm)',data=df) # Label the axes _ = plt.xlabel('species') _ = plt.ylabel('petal length (cm)') # Show the plot plt.show()
我們也可以把CDF圖和對應的分位數標記重疊畫出整合的一個圖來。
# Plot the ECDF _ = plt.plot(x_vers, y_vers, '.') _ = plt.xlabel('petal length (cm)') _ = plt.ylabel('ECDF') # Overlay percentiles as red diamonds. # ptiles_vers = [2.5,25,50,75,97.5] # percentiles = [ 2.5, 25. , 50. , 75. , 97.5] _ = plt.plot(ptiles_vers, percentiles/100, marker='D', color='red', linestyle='none') # Show the plot plt.show()
variance/mse/std
方差反映了資料的分散程度，方差越大，說明變化越大。標準差則是方差取平方根的結果。
# Array of differences to mean: differences means = np.mean(versicolor_petal_length)*np.ones(len(versicolor_petal_length)) differences = versicolor_petal_length - means # Square the differences: diff_sq diff_sq = differences**2 # Compute the mean square difference: variance_explicit variance_explicit = np.sum(diff_sq)/len(versicolor_petal_length) # Compute the variance using NumPy: variance_np variance_np = np.var(versicolor_petal_length)
variancestd_np = np.std(versicolor_petal_length)
# variance_np == variance_explicit print(variance_explicit,variance_np) # Print the results
covriance,協方差係數(Pearson r)及協方差矩陣
$$covariance = \frac{1}{n}\sum_{i=1}^{n} (x_i-\bar x)(y_i-\bar y)$$
$$\rho = Pearson \enspace correlation = \frac{covariance}{(std \enspace of \enspace x)(std \enspace of \enspace y)} \in [-1,1]$$
協方差可以反映兩個變數之間的線性關係，特別是協方差係數由於去除了量綱的影響，更加客觀。當協方差係數為正時為正相關，當等於0時，不相關。下面的程式碼我們先畫出散點圖，直觀看到兩個變數之間的線性關係，隨後通過np.cov來計算起協方差矩陣。
# Make a scatter plot plt.scatter(versicolor_petal_length,versicolor_petal_width) # Label the axes plt.xlabel('length') plt.ylabel('width') # Show the result plt.show() covariance_matrix = np.cov(versicolor_petal_length,versicolor_petal_width)
corr_mat = np.corrcoef(x,y) # 計算皮爾遜協方差係數 print(covariance_matrix)
統計推斷(statistical inferrence)
概率思想的核心就是給定我們一組資料，以概率語言表達：如果我們繼續取樣下一組資料，我們將能期待數字特徵是什麼，可以得到什麼推斷。因為可以肯定的是，當我們拿到下一組資料時，相關統計指標邊然不同。這些統計指標本身就具有不確定性和隨機性，但是卻具有概率意義的穩定性。
統計推斷的目標就是指出如果我們繼續收集更多的資料，則能得出什麼概率性的結論。也就是說，通過較少的資料能夠得出更加通用性的結論。
統計模擬(simulation)
對於概率性的事件，我們可以通過實際實驗來觀察結果，大量的實驗結果就呈現出一定的統計規律性。但是，很多時候重複性的實驗是不現實的，甚至是不可能的，這時，我們就可以通過numpy的random偽隨機引擎來不斷做模擬，以便得到我們感興趣的事件的統計特性。
統計模擬hacker stats probability的步驟
0. 抽象定義清楚你要研究問題的隨機變數；
1. 找到一個如何模擬產生實驗資料從而得到隨機變數值的計算方法；
2. 模擬非常非常多的次數，得到每次的結果；
3. 概率的計算就是使用對應模擬過程結果中屬於我們所感興趣結果的資料發生次數除以實驗總次數
只要我們能夠使用np.random來模擬一個實驗得到結果，那麼我們就可以得到相應的概率分佈PDF
首個使用numpy random模擬一個數字發生器，看看其直方圖中各個bin的分佈情況:
# Seed the random number generator np.random.seed(42) # Initialize random numbers: random_numbers random_numbers = np.empty(1000000) # Generate random numbers by looping over range(1000000) for i in range(1000000): random_numbers[i] = np.random.random(1)[0] # Plot a histogram _ = plt.hist(random_numbers) # Show the plot plt.show()
一個完整的銀行貸款壞賬概率模擬計算例子(基於np.random.random從頭開始做起):
# 模擬n重伯努利實驗函式 def perform_bernoulli_trials(n, p): """Perform n Bernoulli trials with success probability p and return number of successes.""" # Initialize number of successes: n_success n_success = 0 # Perform trials for i in range(n): # Choose random number between zero and one: random_number random_number = np.random.random() # If less than p, it's a success so add one to n_success if random_number < p: n_success +=1 return n_success # Seed random number generator np.random.seed(42) # Initialize the number of defaults: n_defaults n_defaults = np.empty(1000) # 模擬執行1000次100重0.05概率的伯努利實驗，檢視每百次 # 實驗中出現拖欠的次數分佈情況 for i in range(1000): n_defaults[i] = perform_bernoulli_trials(n=100,p=0.05) # 將模擬實驗結果拖欠隨機變數數的直方圖畫出來 _ = plt.hist(x=n_defaults, density=True) _ = plt.xlabel('number of defaults out of 100 loans') _ = plt.ylabel('probability') # Show the plot plt.show() # Compute ECDF: x, y x,y = ecdf(n_defaults) # 計算並且畫出每百次貸款拖欠數這個隨機變數的CDF概率分佈圖 # Plot the ECDF with labeled axes plt.plot(x,y,marker='.',linestyle='none') plt.xlabel('defaults') plt.ylabel('probability') # Show the plot plt.show() # 計算每100次貸款中有超過10次拖欠的次數 # Compute the number of 100-loan simulations with 10 or more defaults: n_lose_money n_defaults2 = n_defaults >=10 n_lose_money = 0 for i in range(len(n_defaults)): n_lose_money += n_defaults2[i] # Compute and print probability of losing money # 計算百次貸款拖欠10人次以上事件的概率 print('Probability of losing money =', n_lose_money / len(n_defaults))
改進型模擬，直接使用np.random.binomial取樣畫圖，得到相同的CDF模擬結果
# Take 10,000 samples out of the binomial distribution: n_defaults n_defaults = np.random.binomial(100,0.05,size=10000) # Compute CDF: x, y x, y = ecdf(n_defaults) # Plot the CDF with axis labels plt.plot(x,y,marker='.',linestyle='none') plt.xlabel('defaults') plt.ylabel('propabability') # Show the plot plt.show()
plotting the binomial PMF(使用直方圖模擬，目前還不會畫純粹的PMF)
# Compute bin edges: bins bins = np.arange(0, max(n_defaults)+1.5) - 0.5 # Generate histogram plt.hist(n_defaults,bins=bins,normed=True) # Label axes plt.xlabel('defaults number') plt.ylabel('proportion') # Show the plot plt.show()
泊松分佈(Poisson process/Poison distribution)
story: 在每特定間隔內發生事件數均值為$\lambda$次的泊松過程中，在給定事件內發生泊松過程r符合泊松分佈。比如：如果一個網站在1小時內平均會有6次訪問，那麼在一個小時內網站的訪問次數這個隨機變數就符合泊松分佈；
泊松分佈可以看做是Binomial分佈的極限：當n重伯努利實驗中事件概率極低，而試驗次數非常巨大時，也就是說針對小概率事件下的伯努利分佈。
由於泊松分佈只有一個引數，相對Binomial伯努利分佈簡單，因此常常可以使用泊松分佈來近似逼近伯努利分佈。
連續性隨機變數分佈
normal distribution
通過使用np.random.normal從相同mean值，但是不同標準差的正態分佈中取樣100000個資料，隨後將這些資料的直方圖畫出來，就可以近似作為正態分佈的密度函式PDF了.
# Draw 100000 samples from Normal distribution with stds of interest: samples_std1, samples_std3, samples_std10 samples_std1 = np.random.normal(20,1,100000) samples_std3 = np.random.normal(20,3,100000) samples_std10 = np.random.normal(20,10,100000) # Make histograms plt.hist(samples_std1,bins = 100,normed = True,histtype='step') plt.hist(samples_std3,bins = 100,normed = True,histtype='step') plt.hist(samples_std10,bins = 100,normed = True,histtype='step') # Make a legend, set limits and show plot _ = plt.legend(('std = 1', 'std = 3', 'std = 10')) plt.ylim(-0.01, 0.42) plt.show() # Generate CDFs x_std1,y_std1 = ecdf(samples_std1) x_std3,y_std3 = ecdf(samples_std3) x_std10,y_std10=ecdf(samples_std10) # Plot CDFs plt.plot(x_std1,y_std1,marker='.',linestyle = 'none') plt.plot(x_std3,y_std3,marker='.',linestyle = 'none') plt.plot(x_std10,y_std10,marker='.',linestyle = 'none') # Make a legend and show the plot _ = plt.legend(('std = 1', 'std = 3', 'std = 10'), loc='lower right') plt.show()
使用樣本均值和樣本方差作為總體均值和方差，探究樣本是否真正符合正態分佈
# Compute mean and standard deviation: mu, sigma mu = np.mean(belmont_no_outliers) sigma = np.std(belmont_no_outliers) # Sample out of a normal distribution with this mu and sigma: samples samples = np.random.normal(mu,sigma,size=10000) # Get the CDF of the samples and of the data x_theor,y_theor = ecdf(samples) x,y = ecdf(belmont_no_outliers) # Plot the CDFs and show the plot _ = plt.plot(x_theor, y_theor) _ = plt.plot(x, y, marker='.', linestyle='none') _ = plt.xlabel('Belmont winning time (sec.)') _ = plt.ylabel('CDF') plt.show()
使用上面的mu,sigma資料定義的正態分佈曲線，我們就假設對應賽馬時間符合該分佈，來計算小於144秒的概率能達到多大
# Take a million samples out of the Normal distribution: samples samples = np.random.normal(mu,sigma,size=1000000) # Compute the fraction that are faster than 144 seconds: prob prob = np.sum(samples<144)/1000000 # Print the result print('Probability of besting Secretariat:', prob) ##### 列印出0.000691
指數分佈 Exponential Distribution
story: 兩個Poison過程發生時的等待時間符合指數分佈。
只有一個引數：平均等待時長
LSR
# Plot the illiteracy rate versus fertility _ = plt.plot(illiteracy, fertility, marker='.', linestyle='none') plt.margins(0.02) _ = plt.xlabel('percent illiterate') _ = plt.ylabel('fertility') # Perform a linear regression using np.polyfit(): a, b a, b = np.polyfit(illiteracy,fertility,1) # Print the results to the screen print('slope =', a, 'children per woman / percent illiterate') print('intercept =', b, 'children per woman') # Make theoretical line to plot x = illiteracy y = a * x + b # Add regression line to your plot _ = plt.plot(x, y) # Draw the plot plt.show()
Bootstrapping
bootstrapping通常使用有放回地重取樣構造的多個資料樣本來執行統計推斷，比如對於均值這個統計量，原來只有一個樣本，因此均值就是從這個樣本通過定義計算出來的，但是我們如何對這個統計量來研究呢？可行的辦法就是通過bootstrapping重新取樣獲取多個樣本，這樣就能計算出多個均值來，我們據此就可以研究這個樣本均值隨機變數的統計特性，比如它的均值，方差，置信區間。。。
Bootstrap Replicate:有放回重抽樣
在前面介紹的統計方法研究資料時，往往對已知的一個資料樣本來計算相應的統計特徵，比如mean, median, mse等。但是我們想一想如果有可能的話，再做了一組實驗得到新的資料，這些統計值是否相同呢？顯然這些統計值本身也是隨機變數，因此新的資料樣本的mean,median肯定是不相同的。那麼我們又有什麼信心就告訴我們目前這個樣本算出來的均值本身就是可靠的呢？ 有沒有什麼辦法模擬出更多類似的樣本呢？一個可行的方法就被稱為有放回抽樣的Bootstrap Replicate方法.Resample的過程就像是在重複相同的實驗一樣不斷獲取新的資料。。。
bootstrap sample:指通過bootstrap有放回抽樣重組的新樣本陣列;
bootstrap replicate:指從resampled array（bootstrap sample）中計算出來的統計量值
通過bootstrap方法創造新資料，和原始dataset重疊比較
for _ in range(50): # Generate bootstrap sample: bs_sample bs_sample = np.random.choice(rainfall, size=rainfall.size) # Compute and plot ECDF from bootstrap sample x, y = ecdf(bs_sample) _ = plt.plot(x, y, marker='.', linestyle='none',color='gray',alpha=0.1) # Compute and plot ECDF from original data x, y = ecdf(rainfall) _ = plt.plot(x, y, marker='.') # Make margins and label axes plt.margins(0.02) _ = plt.xlabel('yearly rainfall (mm)') _ = plt.ylabel('ECDF') # Show the plot plt.show()
資料集統計量置信區間的數值方法
我們知道，對於有命名的著名的統計分佈，可以通過分析的方法計算出對應的置信區間。但是，對於計算機來說，往往更加通用化的做法是使用bootstrapping方法計算出來。。。
通過bootstrap方法創造新的資料集，並且計算對應的均值，方差。。。，得到均值方差的一組值，我們再使用np.percentile(bs_replicates,[2.5,97.5])就可以得到相應均值或者方差的95%置信區間。
# Draw bootstrap replicates of the mean no-hitter time (equal to tau): bs_replicates bs_replicates = draw_bs_reps(nohitter_times,np.mean,size=10000) # Compute the 95% confidence interval: conf_int conf_int = np.percentile(bs_replicates,[2.5,97.5]) # Print the confidence interval print('95% confidence interval =', conf_int, 'games') # Plot the histogram of the replicates _ = plt.hist(bs_replicates, bins=50, normed=True) _ = plt.xlabel(r'$\tau$ (games)') _ = plt.ylabel('PDF') # Show the plot plt.show()
注意：當我們使用bootstrap方法計算統計指標的置信區間時，並不會對隨機變數的分佈做任何假設，沒有任何的未知引數，僅僅是通過數值計算得出。但是當我們做類似線性迴歸模型的引數估計時，我們就有斜率和截距兩個引數，被稱為parametric estimate.這種場景下，最優引數也是由樣本資料來決定的，這一點和均值，方差等統計量隨機變數是一樣的，也就是說，如果我們拿到另一組資料，則我們就可以得到另外一組最優的斜率和截距引數。我們也可以通過bootstrap estimate的方法來計算得到斜率和截距這兩個隨機變數的置信區間，也就是說，我們think probabilistically.
Pairs Bootstrap
以LSR線性迴歸為例，我們要建模投票總數和奧巴馬得票數之間的線性關係，這時如果簡單使用每個州縣隨機的投票總數和奧巴馬得票數可能並不合理，因為這兩者之間本來就具有某種線性關係，如果完全隨機會破壞這種隨機性。
我們使用稱為Pairs Bootstrap的方法來有放回地重新抽樣重組資料，每次resample，得到一組資料，我們按照原始資料集相同個數來做resample獲得bootstrap sample,重新計算斜率和截距，得到一個bootstrap replicate。這個流程重複10000次，就能得到1萬組斜率和截距的估計值。這樣就能夠使用np.percentile有效地計算引數估計的置信區間confidence interval了。
def draw_bs_pairs_linreg(x, y, size=1): """Perform pairs bootstrap for linear regression.""" # Set up array of indices to sample from: inds inds = np.arange(len(x)) # Initialize replicates: bs_slope_reps, bs_intercept_reps bs_slope_reps = np.empty(size) bs_intercept_reps = np.empty(size) # Generate replicates for i in range(size): bs_inds = np.random.choice(inds, size=len(inds)) bs_x, bs_y = x[bs_inds], y[bs_inds] bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x,bs_y,1) return bs_slope_reps, bs_intercept_reps # Generate replicates of slope and intercept using pairs bootstrap bs_slope_reps, bs_intercept_reps = draw_bs_pairs_linreg(illiteracy,fertility,size=1000) # Compute and print 95% CI for slope print(np.percentile(bs_slope_reps,[2.5,97.5])) # Plot the histogram _ = plt.hist(bs_slope_reps, bins=50, normed=True) _ = plt.xlabel('slope') _ = plt.ylabel('PDF') plt.show()
從這個例子我們就能看到文盲率和生育率之間線性關係的slope引數估計的pdf概率分佈，當然就容易計算得出其95%的置信區間.
我們將通過bootstrap方法計算得到的前100個斜率和截距的值所代表的直線繪製出來，就有一個大致的概念，也就是說，如果重新做實驗，則新的直線引數估計形成以下直線集
# Generate array of x-values for bootstrap lines: x x = np.array([0,100]) # Plot the bootstrap lines for i in range(100): _ = plt.plot(x, bs_slope_reps[i]*x + bs_intercept_reps[i], linewidth=0.5, alpha=0.2, color='red') # Plot the data _ = plt.plot(illiteracy,fertility,marker='.',linestyle = 'none') # Label axes, set the margins, and show the plot _ = plt.xlabel('illiteracy') _ = plt.ylabel('fertility') plt.margins(0.02) plt.show()
Formulating and simulating a hypothesis
上面的學習中，我們假設資料符合線性模型，隨後對線性模型的引數進行估計。但是問題是：我們的假設是否合理？這就需要用到統計學中的假設檢驗。我們假設一個推斷是正確的(比如假設變數之間是線性關係，再比如，我們假設乘坐頭等艙和經濟艙的乘客在泰坦尼克沉船事故中存活的概率分佈是相同的)，那麼我們的觀測資料是否合理，多大程度上是合理的？？
permutation:將陣列裡面的資料順序打亂就稱為permutation.
permutation sample:打亂順序後的資料重新擷取形成新的sample就稱為permutation sample.