机器学习初步

linear regression

两个变量之间的相关性越高，用一个来预测另一个变量就越容易。例如，如果我知道我的牛排的价格与牛排的大小成正比，我就可以创建一个公式来帮助我预测我的牛排要付多少钱。

我们这样做的方法是线性回归。线性回归给出了一个公式。如果我们把一个变量的值代入这个公式，我们就得到了另一个变量的值

y = mx + b

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

image.png
我们使用cov函数计算协方差cov(x,y),.var()计算x的方差,其中cov函数返回一个矩阵，必须用0和1来索引

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

然后计算b的值，即截距，

b = image.png
使用定义斜率函数calc_slope来计算截距

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

然后我们可以测试一下预测所有的x值对应的y值，预测的好与坏完全取决于两个变量的相关性大小

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):              #通过x的值计算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线性回归函数

使用linregress函数使线性回归变得简单
linregress(x,y)
能够返回5个参数，分别是slope（斜率）, intercept（截距）, r_value（相关性）, p_value, stderr_slope

residuals(残差)

残差描述了每一个真实值和预测值之间的差距，如果把残差的平方sum起来就得到残差平方和，yi是预测值，最终得到的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)) #为了查看写法是否正确看一下类型，确实是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(标准错误值)

standard error标准误差跟标准差类似，可以让我们快速判断一个线性模型在预测中的好坏，其定义如下

standard error

n-2这是由于整个人群和样本之间的差异造成的
下面的函数举例来查看数据中差值在1个2个或3个标准错误值范围内的数据比例各是多少

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一定要加括号，标准错误值
print(stderr)
def within_proportion(y,y_pre,stderr,num): #定义一个求不同个数个标准错误值内的比例有多少的函数
    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随机种子
随机种子为下一次随机规定随机策略，种子一致，随机到的数列就一致，如下

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随机抽样
在若总数中抽取若干个数为样本,同样需要首先确定随机种子

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

output
[200, 300, 20, 600]

随机抽样的优势在于抽取的数据能够以正态分布，不会出现集中某段数据的状态，如下列代码中求income表中随机抽取100行为一个样本，重复取1000次作为一个随机样本数列，画出这个数列就能发现其实以正态分布分布。

import random
def select_random_sample(count):
    random_indices = random.sample(range(0, income.shape[0]), count)#这部分调用random.sample函数随机取count个结果，作为一个样本数列
    return income.iloc[random_indices]      #iloc[行]组成新df

random.seed(1)

random_sample = [select_random_sample(100)["median_income"].mean() for _ in range(10000)]#求一个样本数列的mean以后，取1000次这样的样本数列平均值
plt.hist(random_sample,20)
plt.show()

output

random_sample

二项分布概率计算

math框架
计算数学公式使用，计算连接概率时需要计算组合的数量时，可以调用math函数中阶乘计算方法math.factorial，众所周知组合的公式如下：

combinations

计算10天中7天晴天的组合数，用代码完成公式如下：

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)

单独组合的概率为

single combination

某一事件发生在N次中发生k次的概率为

multi combination
我们自定义求解某个时间发生若干次的概率分布函数基本如下,先设定一个list范围30，统计在30天中多少天使用bike的人超过5000人的概率分布

import math
import matplotlib.pyplot as plt
outcome_counts = list(range(31)) #创建一个0-30的list
def probability(p,k,N):#计算单独事件发生概率p，在N次中发生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函数来分析二项分布概率，来简化上方我们自己来写的函数，binom函数需要三个变量，x为横轴list，n为计数总次数，p为单次发生的概率，具体函数如下,得到和上图一致的分布图，只要p是不变值，概率分布曲线看起来都是一样的：

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

概率分布标准差

standard deviation of a probability distribution

概率累计函数（cumulative density function）

如果一个硬币投掷三次概率分布是

image.png

那么投掷三次的累计概率分布就是,将低于自己的次数的概率累计相加，到最后一种可能就变成1

image.png
我们使用binom.cdf函数来计算概率累计分布，同binom.pmf也加载X,N,P三个参数，如下计算30天内骑车超过5000人的累计概率分布图

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天时候概率就累计为1了，说明超过20天或者更多天数出现5000人骑车的可能性已经完全没有了.

如果要计算在某一个值累计的概率的值，可以在X位置只代入一个参数，比如求16天的概率累计值,1-所得值就得到了16天右侧的概率累计值

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（显著性鉴定）

显著性鉴定用来分析两组盲眼测试的结果是能够接受实验有效的alternative hypothesis（则一假设），还是没有效果的null hypothesis（空假设），blind experiment（盲眼实验）的目的能够规避掉人的心理对于测试影响，即the control group（对照组）和the treatment group（实验组）不知道自己是实验组还是对照组。通过一个减肥药是否对固定人群体重变化有效的实验来进一步说明显著性鉴定。

首先要进行research design，他能够帮助我们在更深入研究时候克服研究中的缺陷。如果在实验中有一个有显著意义的区别，我们叫做这叫做statistically significant（统计学意义），对于减肥药实验我们分为两组a和b来进行对比，其中具有统计学意义的变量我们认为是平均体重的下降量

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，为了检测出我们这一组数据是一种很轻易就能随机出现的差值，还是一个极具特殊性的很低几率出现的值，我们通过采样分布来观察这组数据是哪一种，计算采样分布我们使用的方法是将两组数据统一混合后，分辨随机再分为两组，然后对两组数据计算平均差值，取1000次样本的差值分布数组，画出分布图

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

这时候就可以看出我们在1000次样本的情况下，均值差从未打到过2.52这么大，当然你可以说是采样次数不足，在进行采样设计时我们需要权衡样本数量和计算速度的数值。

之后我们把mean_differences中没一个值出现的次数统计成一个字典，使用get（a，False）函数如果在字典key = a有值就返回value，如果没有值就返回False，来做一个统计。

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)

我们可以通过计算超过2.52平均值的次数出现的概率，等于超过2.52的次数总和/1000次，这个概率我们称作p value，通常我们通过设定一个p value threshold来衡量事件出现是偏向于空假设还是则一假设。我们来计算一下我们的数据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值小于阈值：

• 拒绝零假设，即两组参与者减去的平均体重没有差别，
• 接受则一假设，即吃减肥药的人体重减轻了，
• 结论：减肥药丸确实会影响减肥者的体重。（在这个例子中体现为很低概率出现平均值达到实验组数据体重平均值差值的情况，实验为非常特殊的数据）

如果p值大于临界值，：

• 接受零假设，即两组参与者减去的平均体重没有差别，
• 拒绝则一假设即那些服用减肥药丸的人体重减轻了，
• 结论：减肥药丸似乎并不能有效地帮助人们减轻体重（在这个例子中体现为很高概率出现平均值达到实验组数据体重平均值差值的情况）。

我们通常p值的阈值设定为5%,即0.05，也就是说，只有5%的概率是随机的，大多数研究人员都认为是可以接受的。
我们的p值最终为0,我们倾向于相信这是一组则一假设，很低概率再出现如此差值的事件，不是简单的偶然事件。

最终我们的实验设计和p阈值设定是对显著性鉴定有着重大影响

• 研究设计是非常重要的，会影响你的结果。例如，如果A组的参与者意识到他们服用了安慰剂，他们可能会改变他们的行为并影响结果。
• 你设置的p值阈值也会影响你到达的结论。
  如果你设置的p值阈值过高，你可能会错误地接受另一个假设，拒绝零假设。这就是所谓的type I error。
  如果你设置的p值阈值太低，你可能会不正确地拒绝另一种假设，而选择接受零假设。这是type II error。

Chi-squared tests 卡方测试

卡方测试使我们能够确定观察一组分类值的统计意义，比如我们观察到两组分类数据男人和女人在全美收入调查表(US income and demographics)的的结果分布如下

male and female
我们知道有些东西看起来不太好，为什么男的这么多女的这么少，但是我们不知道如何量化观察和期望的值有多么不同。我们也没有办法确定两个组之间是否存在统计上的显著差异，如果我们需要进一步调查的话。
这就是卡方测试可以提供帮助的地方。卡方测试使我们能够**量化观察到的和期望的分类值之间的

卡方测试的公式举例显示为

image.png

平方误差可以确保所有的差异不等于零（你不能有负的平方），给一个非零数，可以用来评估统计意义。
然后对每一个分类数据的卡方值求和，就能够得到组数据中的卡方和

我们来计算1000次上面的数据表中的卡方分布

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

for _ in range(1000):    #循环1000次随机男女类别的数据分配
    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判定为男,统计随机生成的男人数量
    female_count = len(vector[vector >= 0.5]) #大于0.5判定为女，统计随机生成的女人数量
    male_diff = (male_count - 16280.5)**2 /16280.5  #计算一次随机数的男卡方值
    female_diff = (female_count -16280.5)**2 /16280.5  #计算一次随机数的女卡方值
    total = male_diff + female_diff#求一次的卡方和
    chi_squared_values.append(total)#变成卡方数组
plt.hist(chi_squared_values)#画出卡方分布图
plt.show()

chis_quare

可以从图中看出卡方值大于8的数几乎就聊聊无几了，基本上都集中在0-2，这说明分布值远离期望的随机分布出现的数量极其少。

• 从下面卡方和的公式可以看出来，卡方值是随着样本数量的增加线性提升的

  
  image.png
  

  比如样本数量提高了10倍，再左侧分子部分值就扩大了100倍，再除以各自分母后总值就是扩大了10倍

• 但从另一方面说，卡方值会在数据量很少向很多提升的时候，更趋近与正态分布，又会有一定的减小，比如扔色子

  
  flipping a coin 10 times
  在扔的总数少的时候可能出现与总数偏离比较严重的数据分布，但随着样本数量的提升，这种情况可能会很少见 flipping a coin 1000 times
  如果不是硬币出老千的话，正常值不会有如此多的极端分布，所以在一定范围内卡方值会有变小的趋势，这两种效应相互抵消，而在每次迭代中采样200个项目时，所构建的一个卡方抽样分布看起来与一个抽样千项相同。

degree of freedom(自由度)
之前我们讨论的男女问题，只有两个分类，但其实只有一个自由度，因为总数减去女人就是男人，所以在多自由度分类数据中，也可以同样运用卡方测试，下面是各种族收入表，我们计算一下他的卡方值

image.png
可以使用函数chisquare_value,pvalue，来直接计算两个nparray的卡方值和pvalue，上一节讲过pvalue，如果小于百分之5的数据可以被证实是有统计显著性，大于百分之5更多的被认为是随机也能够产生的结果。
chisquare_value, pvalue = chisquare(observed, expected)

详细代码如下

from scipy.stats import chisquare
import numpy as np
observed = np.array([27816,3124,1039,311,271]) #转换为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计算可能选取了非常大的样本数量才能得到这样小的值，但小于百分之5依旧足以证明这组数据是有统计显著性的