linear regression

兩個(gè)變量之間的相關(guān)性越高，用一個(gè)來(lái)預(yù)測(cè)另一個(gè)變量就越容易。例如，如果我知道我的牛排的價(jià)格與牛排的大小成正比，我就可以創(chuàng)建一個(gè)公式來(lái)幫助我預(yù)測(cè)我的牛排要付多少錢。

我們這樣做的方法是線性回歸。線性回歸給出了一個(gè)公式。如果我們把一個(gè)變量的值代入這個(gè)公式，我們就得到了另一個(gè)變量的值

y = mx + b

m 的值就是斜率（slope），b的值就是截距（intercept）
m =

image.png

我們使用cov函數(shù)計(jì)算協(xié)方差cov(x,y),.var()計(jì)算x的方差,其中cov函數(shù)返回一個(gè)矩陣，必須用0和1來(lái)索引

from numpy import cov
slope_density = cov(wine_quality["density"],wine_quality["quality"])[0,1] / wine_quality["density"].var()
print(slope_density)

output
-90.9423999421

然后計(jì)算b的值，即截距，

b =

image.png

使用定義斜率函數(shù)calc_slope來(lái)計(jì)算截距

from numpy import cov

# This function will take in two columns of data, and return the slope of the linear regression line.
def calc_slope(x, y):
    return cov(x, y)[0, 1] / x.var()

intercept_density = wine_quality["quality"].mean() - calc_slope(wine_quality["density"],wine_quality["quality"]) * wine_quality["density"].mean()

print(intercept_density)

output
96.2771445761

然后我們可以測(cè)試一下預(yù)測(cè)所有的x值對(duì)應(yīng)的y值，預(yù)測(cè)的好與壞完全取決于兩個(gè)變量的相關(guān)性大小

from numpy import cov

def calc_slope(x, y):
    return cov(x, y)[0, 1] / x.var()

# Calculate the intercept given the x column, y column, and the slope
def calc_intercept(x, y, slope):
    return y.mean() - (slope * x.mean())

def calc_y(x):              #通過(guò)x的值計(jì)算y值
    y = m * x + b
    return y

m = calc_slope(wine_quality["density"],wine_quality["quality"])
b = calc_intercept(wine_quality["density"],wine_quality["quality"],m)
predicted_quality = wine_quality["density"].apply(calc_y)
print(predicted_quality)

output
0       5.243802
1       5.880399
2       5.780362
3       5.734891
4       5.734891
5       5.780362
6       5.798551
7       5.243802

linregress線性回歸函數(shù)

使用linregress函數(shù)使線性回歸變得簡(jiǎn)單
linregress(x,y)
能夠返回5個(gè)參數(shù)，分別是slope（斜率）, intercept（截距）, r_value（相關(guān)性）, p_value, stderr_slope

residuals(殘差)

殘差描述了每一個(gè)真實(shí)值和預(yù)測(cè)值之間的差距，如果把殘差的平方sum起來(lái)就得到殘差平方和，yi是預(yù)測(cè)值，最終得到的RSS就是殘差和

sum of squared residuals

from scipy.stats import linregress
slope, intercept, r_value, p_value, stderr_slope = linregress(wine_quality["density"], wine_quality["quality"])
print(slope)
print(intercept)
y = slope * wine_quality["density"] + intercept
print(type(y)) #為了查看寫法是否正確看一下類型，確實(shí)是series型
residuals = (wine_quality["quality"] - y) ** 2
rss = sum(residuals)
print(rss)    


output
-90.9423999421
96.2771445761
<class 'pandas.core.series.Series'>
3478.68946969

standard error(標(biāo)準(zhǔn)錯(cuò)誤值)

standard error標(biāo)準(zhǔn)誤差跟標(biāo)準(zhǔn)差類似，可以讓我們快速判斷一個(gè)線性模型在預(yù)測(cè)中的好壞，其定義如下

standard error

n-2這是由于整個(gè)人群和樣本之間的差異造成的
下面的函數(shù)舉例來(lái)查看數(shù)據(jù)中差值在1個(gè)2個(gè)或3個(gè)標(biāo)準(zhǔn)錯(cuò)誤值范圍內(nèi)的數(shù)據(jù)比例各是多少

from scipy.stats import linregress
import numpy as np
slope, intercept, r_value, p_value, stderr_slope = linregress(wine_quality["density"], wine_quality["quality"])
predicted_y = np.asarray([slope * x + intercept for x in wine_quality["density"]])
residuals = (wine_quality["quality"] - predicted_y) ** 2
rss = sum(residuals)  #殘差和
stderr = (rss/ len(wine_quality["quality"] -2) ) ** (1/2) #這里1/2一定要加括號(hào)，標(biāo)準(zhǔn)錯(cuò)誤值
print(stderr)
def within_proportion(y,y_pre,stderr,num): #定義一個(gè)求不同個(gè)數(shù)個(gè)標(biāo)準(zhǔn)錯(cuò)誤值內(nèi)的比例有多少的函數(shù)
    diff = abs(y_pre - y)
    within = [item for item in diff if item <= (stderr * num)]
    proportion = len(within) / len(y)
    return proportion
 
within_one = within_proportion(wine_quality["quality"],predicted_y,stderr,1)
within_two = within_proportion(wine_quality["quality"],predicted_y,stderr,2)
within_three = within_proportion(wine_quality["quality"],predicted_y,stderr,3)
print(within_one,within_two,within_three)
 
output
0.842749378428
0.6845651286239282 0.9356880359330338 0.9936708860759493

Distributions and sampling

random.seed隨機(jī)種子
隨機(jī)種子為下一次隨機(jī)規(guī)定隨機(jī)策略，種子一致，隨機(jī)到的數(shù)列就一致，如下

random.seed(10)
print([random.randint(0,10) for _ in range(5)])
random.seed(10)
# Same sequence as above.
print([random.randint(0,10) for _ in range(5)])
random.seed(20)
new_sequence = [random.randint(0,10) for _ in range(5)]
print(new_sequence)

output
[9, 0, 6, 7, 9]
[9, 0, 6, 7, 9]
[10, 2, 4, 10, 10]

random.sample隨機(jī)抽樣
在若總數(shù)中抽取若干個(gè)數(shù)為樣本,同樣需要首先確定隨機(jī)種子

shopping = [300, 200, 100, 600, 20]
random.seed(1)
new_sample = random.sample(shopping,4)
print(new_sample)

output
[200, 300, 20, 600]

隨機(jī)抽樣的優(yōu)勢(shì)在于抽取的數(shù)據(jù)能夠以正態(tài)分布，不會(huì)出現(xiàn)集中某段數(shù)據(jù)的狀態(tài)，如下列代碼中求income表中隨機(jī)抽取100行為一個(gè)樣本，重復(fù)取1000次作為一個(gè)隨機(jī)樣本數(shù)列，畫出這個(gè)數(shù)列就能發(fā)現(xiàn)其實(shí)以正態(tài)分布分布。

import random
def select_random_sample(count):
    random_indices = random.sample(range(0, income.shape[0]), count)#這部分調(diào)用random.sample函數(shù)隨機(jī)取count個(gè)結(jié)果，作為一個(gè)樣本數(shù)列
    return income.iloc[random_indices]      #iloc[行]組成新df

random.seed(1)

random_sample = [select_random_sample(100)["median_income"].mean() for _ in range(10000)]#求一個(gè)樣本數(shù)列的mean以后，取1000次這樣的樣本數(shù)列平均值
plt.hist(random_sample,20)
plt.show()

output

random_sample

二項(xiàng)分布概率計(jì)算

math框架
計(jì)算數(shù)學(xué)公式使用，計(jì)算連接概率時(shí)需要計(jì)算組合的數(shù)量時(shí)，可以調(diào)用math函數(shù)中階乘計(jì)算方法math.factorial，眾所周知組合的公式如下：

combinations

計(jì)算10天中7天晴天的組合數(shù)，用代碼完成公式如下：

import math
def find_outcome_combinations(N, k):
    # Calculate the numerator of our formula.
    numerator = math.factorial(N)
    # Calculate the denominator.
    denominator = math.factorial(k) * math.factorial(N - k)
    # Divide them to get the final value.
    return numerator / denominator
sunny7 = find_outcome_combinations(10,7)

單獨(dú)組合的概率為

single combination

某一事件發(fā)生在N次中發(fā)生k次的概率為

multi combination

我們自定義求解某個(gè)時(shí)間發(fā)生若干次的概率分布函數(shù)基本如下,先設(shè)定一個(gè)list范圍30，統(tǒng)計(jì)在30天中多少天使用bike的人超過(guò)5000人的概率分布

import math
import matplotlib.pyplot as plt
outcome_counts = list(range(31)) #創(chuàng)建一個(gè)0-30的list
def probability(p,k,N):#計(jì)算單獨(dú)事件發(fā)生概率p，在N次中發(fā)生k次的概率
    q = 1-p
    single = p ** k * q**(N-k) 
    combination = math.factorial(N)/(math.factorial(k)*math.factorial(N-k))
    pro = single*combination
    return pro
p = 0.39
outcome_probs = [probability(p,item,30) for item in outcome_counts] #組成可能性分布list
print(outcome_probs)
plt.bar(outcome_counts, outcome_probs)
plt.show()

image.png

Scipy框架
python中使用Scipy框架中中的binom.pmf函數(shù)來(lái)分析二項(xiàng)分布概率，來(lái)簡(jiǎn)化上方我們自己來(lái)寫的函數(shù)，binom函數(shù)需要三個(gè)變量，x為橫軸list，n為計(jì)數(shù)總次數(shù)，p為單次發(fā)生的概率，具體函數(shù)如下,得到和上圖一致的分布圖，只要p是不變值，概率分布曲線看起來(lái)都是一樣的：

import scipy
from scipy import linspace
from scipy.stat import binom
outcome = linspace(0,30,31)
distribution = binom.pmf(outcome,30,0.39)
print(type(distribution))
plt.bar(outcome,distribution)
plt.show()

output
<class 'numpy.ndarray'>

image.png

概率分布期望值

expected value

概率分布標(biāo)準(zhǔn)差

standard deviation of a probability distribution

概率累計(jì)函數(shù)（cumulative density function）

如果一個(gè)硬幣投擲三次概率分布是

image.png

那么投擲三次的累計(jì)概率分布就是,將低于自己的次數(shù)的概率累計(jì)相加，到最后一種可能就變成1

image.png

我們使用binom.cdf函數(shù)來(lái)計(jì)算概率累計(jì)分布，同binom.pmf也加載X,N,P三個(gè)參數(shù)，如下計(jì)算30天內(nèi)騎車超過(guò)5000人的累計(jì)概率分布圖

from scipy import linspace
from scipy.stats import binom
outcome_counts = linspace(0,30,31)
comulative_dis = binom.cdf(outcome_counts,30,0.39)
plt.plot(outcome_counts,comulative_dis)
plt.show()

cdf

我們能夠看出最終在不到20天時(shí)候概率就累計(jì)為1了，說(shuō)明超過(guò)20天或者更多天數(shù)出現(xiàn)5000人騎車的可能性已經(jīng)完全沒(méi)有了.

如果要計(jì)算在某一個(gè)值累計(jì)的概率的值，可以在X位置只代入一個(gè)參數(shù)，比如求16天的概率累計(jì)值,1-所得值就得到了16天右側(cè)的概率累計(jì)值

left_16 = None
right_16 = None
left_16 = binom.cdf(16,30,.39)
print(left_16)
right_16 = 1 - left_16

output
0.962300376605

Significance Testing（顯著性鑒定）

顯著性鑒定用來(lái)分析兩組盲眼測(cè)試的結(jié)果是能夠接受實(shí)驗(yàn)有效的alternative hypothesis（則一假設(shè)），還是沒(méi)有效果的null hypothesis（空假設(shè)），blind experiment（盲眼實(shí)驗(yàn)）的目的能夠規(guī)避掉人的心理對(duì)于測(cè)試影響，即the control group（對(duì)照組）和the treatment group（實(shí)驗(yàn)組）不知道自己是實(shí)驗(yàn)組還是對(duì)照組。通過(guò)一個(gè)減肥藥是否對(duì)固定人群體重變化有效的實(shí)驗(yàn)來(lái)進(jìn)一步說(shuō)明顯著性鑒定。

首先要進(jìn)行research design，他能夠幫助我們?cè)诟钊胙芯繒r(shí)候克服研究中的缺陷。如果在實(shí)驗(yàn)中有一個(gè)有顯著意義的區(qū)別，我們叫做這叫做statistically significant（統(tǒng)計(jì)學(xué)意義），對(duì)于減肥藥實(shí)驗(yàn)我們分為兩組a和b來(lái)進(jìn)行對(duì)比，其中具有統(tǒng)計(jì)學(xué)意義的變量我們認(rèn)為是平均體重的下降量

import numpy as np
import matplotlib.pyplot as plt
mean_group_a = np.mean(weight_lost_a)
mean_group_b = np.mean(weight_lost_b)
print(mean_group_a,mean_group_b)
plt.hist(weight_lost_a,20)
plt.show()
plt.hist(weight_lost_b,20)
plt.show()

output
2.82 5.34

能看到組a下降平均量為2.82，組b為5.34，差值為2.52，為了檢測(cè)出我們這一組數(shù)據(jù)是一種很輕易就能隨機(jī)出現(xiàn)的差值，還是一個(gè)極具特殊性的很低幾率出現(xiàn)的值，我們通過(guò)采樣分布來(lái)觀察這組數(shù)據(jù)是哪一種，計(jì)算采樣分布我們使用的方法是將兩組數(shù)據(jù)統(tǒng)一混合后，分辨隨機(jī)再分為兩組，然后對(duì)兩組數(shù)據(jù)計(jì)算平均差值，取1000次樣本的差值分布數(shù)組，畫出分布圖

mean_difference = 2.52
print(all_values)
mean_differences = []
for _ in range(1000):
    group_a = []
    group_b = []
    for item in all_values:
        a = numpy.random.rand()
        if a >= 0.5:
            group_a.append(item)
        else:
            group_b.append(item)
    iteration_mean_difference = numpy.mean(group_b) - numpy.mean(group_a)
    mean_differences.append(iteration_mean_difference)
#print(mean_differences)
plt.hist(mean_differences,20)
plt.show()

sampling distribution

這時(shí)候就可以看出我們?cè)?000次樣本的情況下，均值差從未打到過(guò)2.52這么大，當(dāng)然你可以說(shuō)是采樣次數(shù)不足，在進(jìn)行采樣設(shè)計(jì)時(shí)我們需要權(quán)衡樣本數(shù)量和計(jì)算速度的數(shù)值。

之后我們把mean_differences中沒(méi)一個(gè)值出現(xiàn)的次數(shù)統(tǒng)計(jì)成一個(gè)字典，使用get（a，F(xiàn)alse）函數(shù)如果在字典key = a有值就返回value，如果沒(méi)有值就返回False，來(lái)做一個(gè)統(tǒng)計(jì)。

sampling_distribution = {}
for item in mean_differences:
    if sampling_distribution.get(item,False):
        var = sampling_distribution[item]
        var += 1
        sampling_distribution[item] = var
    else:
        sampling_distribution[item] = 1
print(sampling_distribution)

我們可以通過(guò)計(jì)算超過(guò)2.52平均值的次數(shù)出現(xiàn)的概率，等于超過(guò)2.52的次數(shù)總和/1000次，這個(gè)概率我們稱作p value，通常我們通過(guò)設(shè)定一個(gè)p value threshold來(lái)衡量事件出現(xiàn)是偏向于空假設(shè)還是則一假設(shè)。我們來(lái)計(jì)算一下我們的數(shù)據(jù)p值

frequencies = []
for key in sampling_distribution.keys():
    if key >=2.52:
        frequencies.append(sampling_distribution.get(key))
p_value = sum(frequencies)/1000
print(p_value)

output
0.0

如果p值小于閾值：

• 拒絕零假設(shè)，即兩組參與者減去的平均體重沒(méi)有差別，
• 接受則一假設(shè)，即吃減肥藥的人體重減輕了，
• 結(jié)論：減肥藥丸確實(shí)會(huì)影響減肥者的體重。（在這個(gè)例子中體現(xiàn)為很低概率出現(xiàn)平均值達(dá)到實(shí)驗(yàn)組數(shù)據(jù)體重平均值差值的情況，實(shí)驗(yàn)為非常特殊的數(shù)據(jù)）

如果p值大于臨界值，：

• 接受零假設(shè)，即兩組參與者減去的平均體重沒(méi)有差別，
• 拒絕則一假設(shè)即那些服用減肥藥丸的人體重減輕了，
• 結(jié)論：減肥藥丸似乎并不能有效地幫助人們減輕體重（在這個(gè)例子中體現(xiàn)為很高概率出現(xiàn)平均值達(dá)到實(shí)驗(yàn)組數(shù)據(jù)體重平均值差值的情況）。

我們通常p值的閾值設(shè)定為5%,即0.05，也就是說(shuō)，只有5%的概率是隨機(jī)的，大多數(shù)研究人員都認(rèn)為是可以接受的。
我們的p值最終為0,我們傾向于相信這是一組則一假設(shè)，很低概率再出現(xiàn)如此差值的事件，不是簡(jiǎn)單的偶然事件。

最終我們的實(shí)驗(yàn)設(shè)計(jì)和p閾值設(shè)定是對(duì)顯著性鑒定有著重大影響

• 研究設(shè)計(jì)是非常重要的，會(huì)影響你的結(jié)果。例如，如果A組的參與者意識(shí)到他們服用了安慰劑，他們可能會(huì)改變他們的行為并影響結(jié)果。
• 你設(shè)置的p值閾值也會(huì)影響你到達(dá)的結(jié)論。
  如果你設(shè)置的p值閾值過(guò)高，你可能會(huì)錯(cuò)誤地接受另一個(gè)假設(shè)，拒絕零假設(shè)。這就是所謂的type I error。
  如果你設(shè)置的p值閾值太低，你可能會(huì)不正確地拒絕另一種假設(shè)，而選擇接受零假設(shè)。這是type II error。

Chi-squared tests 卡方測(cè)試

卡方測(cè)試使我們能夠確定觀察一組分類值的統(tǒng)計(jì)意義，比如我們觀察到兩組分類數(shù)據(jù)男人和女人在全美收入調(diào)查表(US income and demographics)的的結(jié)果分布如下

male and female

我們知道有些東西看起來(lái)不太好，為什么男的這么多女的這么少，但是我們不知道如何量化觀察和期望的值有多么不同。我們也沒(méi)有辦法確定兩個(gè)組之間是否存在統(tǒng)計(jì)上的顯著差異，如果我們需要進(jìn)一步調(diào)查的話。
這就是卡方測(cè)試可以提供幫助的地方?？ǚ綔y(cè)試使我們能夠**量化觀察到的和期望的分類值之間的

卡方測(cè)試的公式舉例顯示為

image.png

平方誤差可以確保所有的差異不等于零（你不能有負(fù)的平方），給一個(gè)非零數(shù)，可以用來(lái)評(píng)估統(tǒng)計(jì)意義。
然后對(duì)每一個(gè)分類數(shù)據(jù)的卡方值求和，就能夠得到組數(shù)據(jù)中的卡方和

我們來(lái)計(jì)算1000次上面的數(shù)據(jù)表中的卡方分布

import numpy as np
from numpy.random import random
chi_squared_values = []

for _ in range(1000):    #循環(huán)1000次隨機(jī)男女類別的數(shù)據(jù)分配
    male_count = 0
    female_count = 0
    vector = random((32561,))
    #print(type(vector))
    #vector[vector < 0.5] = 0
    #vector[vector >=.5] = 1
    male_count = len(vector[vector < 0.5])   #小于0.5判定為男,統(tǒng)計(jì)隨機(jī)生成的男人數(shù)量
    female_count = len(vector[vector >= 0.5]) #大于0.5判定為女，統(tǒng)計(jì)隨機(jī)生成的女人數(shù)量
    male_diff = (male_count - 16280.5)**2 /16280.5  #計(jì)算一次隨機(jī)數(shù)的男卡方值
    female_diff = (female_count -16280.5)**2 /16280.5  #計(jì)算一次隨機(jī)數(shù)的女卡方值
    total = male_diff + female_diff#求一次的卡方和
    chi_squared_values.append(total)#變成卡方數(shù)組
plt.hist(chi_squared_values)#畫出卡方分布圖
plt.show()

chis_quare

可以從圖中看出卡方值大于8的數(shù)幾乎就聊聊無(wú)幾了，基本上都集中在0-2，這說(shuō)明分布值遠(yuǎn)離期望的隨機(jī)分布出現(xiàn)的數(shù)量極其少。

• 從下面卡方和的公式可以看出來(lái)，卡方值是隨著樣本數(shù)量的增加線性提升的

  
  
  image.png
  
  

  比如樣本數(shù)量提高了10倍，再左側(cè)分子部分值就擴(kuò)大了100倍，再除以各自分母后總值就是擴(kuò)大了10倍

• 但從另一方面說(shuō)，卡方值會(huì)在數(shù)據(jù)量很少向很多提升的時(shí)候，更趨近與正態(tài)分布，又會(huì)有一定的減小，比如扔色子

  
  
  flipping a coin 10 times
  
  
  在扔的總數(shù)少的時(shí)候可能出現(xiàn)與總數(shù)偏離比較嚴(yán)重的數(shù)據(jù)分布，但隨著樣本數(shù)量的提升，這種情況可能會(huì)很少見
  flipping a coin 1000 times
  
  如果不是硬幣出老千的話，正常值不會(huì)有如此多的極端分布，所以在一定范圍內(nèi)卡方值會(huì)有變小的趨勢(shì)，這兩種效應(yīng)相互抵消，而在每次迭代中采樣200個(gè)項(xiàng)目時(shí)，所構(gòu)建的一個(gè)卡方抽樣分布看起來(lái)與一個(gè)抽樣千項(xiàng)相同。

degree of freedom(自由度)
之前我們討論的男女問(wèn)題，只有兩個(gè)分類，但其實(shí)只有一個(gè)自由度，因?yàn)榭倲?shù)減去女人就是男人，所以在多自由度分類數(shù)據(jù)中，也可以同樣運(yùn)用卡方測(cè)試，下面是各種族收入表，我們計(jì)算一下他的卡方值

image.png

可以使用函數(shù)chisquare_value,pvalue，來(lái)直接計(jì)算兩個(gè)nparray的卡方值和pvalue，上一節(jié)講過(guò)pvalue，如果小于百分之5的數(shù)據(jù)可以被證實(shí)是有統(tǒng)計(jì)顯著性，大于百分之5更多的被認(rèn)為是隨機(jī)也能夠產(chǎn)生的結(jié)果。
chisquare_value, pvalue = chisquare(observed, expected)

詳細(xì)代碼如下

from scipy.stats import chisquare
import numpy as np
observed = np.array([27816,3124,1039,311,271]) #轉(zhuǎn)換為nparray，list是不可以的
expected = np.array([26146.5,3939.9,944.3,260.5,1269.8])
chisquare_value,race_pvalue = chisquare(observed,expected)
print(chisquare_value,race_pvalue)

output
1080.48593659 1.28484946749e-232

卡方分布的p值非常小，這里的pvalue計(jì)算可能選取了非常大的樣本數(shù)量才能得到這樣小的值，但小于百分之5依舊足以證明這組數(shù)據(jù)是有統(tǒng)計(jì)顯著性的