Python数据分析与挖掘
1、课程介绍
学习时间从2021.12.3-2021.12.15
这套笔记学习记录为《从零开始学Python数据分析与挖掘》（作者：刘顺祥）的学习。
重新熟悉Python并且巩固数据挖掘知识
1.1、资料
数据源和源代码：
源码数据资料github版（缺少第14章，无PPT）：https://github.com/SnakeLiu/Python-Data-Aanalysis-and-Miner
源码数据资料微信公众号版（完整版，带PPT）：https://www.wqyunpan.com/resourceDetail.html?id=123847&openId=oUgl9weCE59I_Q7hw-TQmM-z9kXM&qrcodeId=108103&sign=ZTVkOGM1MGRhZDJjLTE2Mzg1MjMzOTUzNzA=
百度云链接（自己存储备份的完整版，和微信公众号一样）：https://pan.baidu.com/s/1M0hYT_Fxs2UNKgJM69eqFg
提取码：tawj
视频课（仅供自己学习使用，侵权删）：
百度云链接：
链接：https://pan.baidu.com/s/11-I_f9cOFvaB6kO34ePEgg
提取码：nwiy
阿里云链接：
链接：https://www.aliyundrive.com/s/DLSvv2U3WZZ
1.2、工具选择与安装
python3、anaconda、jupyter notebok
之前安装过，安装过程很简单，可以直接下载anaconda进行安装，已经有python环境的的话也可直接安装jupyter notebok
# 安装
pip3 install jupyter
# 启动
jupyter notebok
在本机寻找合适目录D:\study\pylearn下，开启cmd，进入notebook，在此目录下进行学习。
我直接使用VsCode的jupyter插件进行开发
打开VsCode，Ctrl+Shift+P，打开所有命令窗口，然后输入>jupyter，选择创建新的空白jupyter笔记本
显示所有命令
新建jupyter
可以直接使用了
快捷键：
# 代码运行快捷键
Ctrl + Enter # 运行代码后会停留在当前窗口
Shift + Enter # 运行代码后会跳转到下一个窗口
# 代码注释
Ctrl + /
# 帮助查询（可以连续按好几次，可变换弹出的内容）
Shift + Tab
1.3、数据分析师的一天
1.3.1、数据分析挖掘SEMMA步骤分析
- Sample（搜集数据）：问卷，数据库，实验，仪器
- Explore（数据探索）：离散，连续，异常缺失，特征选择
- Modify （数据修正）：类型转换，数据一致性，异常值和缺失值，数据形态的转换log变化等
- Model （数据建模）：有监督的预测（回归，决策树，KNN），有监督的判别（Logistic聚类，贝叶斯，集成算法），无监督（Kmeans聚类，层次聚类，密度聚类），半监督（关联规则）
- Assess （模型评估）：RMSE，混淆矩阵，ROC曲线，KS曲线
1.3.2、数据分析与挖掘的区别
1.3.3、数据分析师的技能
数据搜集：SQL
数据清洗与探索：Python
数据建模：Python
结果：PPT
2、Python基础
2.1、基础数据结构
2.1.1、基本数据类型
基本类型
int 转整型，字符串的整数，浮点数
float 转浮点型
str 转字符型
bool 转布尔型 True False 非0非空即为真 None会返回假
strptime 转日期类型 datetime模块
字符串转int
# 字符型的3.88转换为int怎么转换
# 先将str转换为float,然后再将float转换为int
int(float('3.88'))
out:
3
日期测试
# 日期测试
import datetime
print(datetime.datetime.strptime('1999-7-1','%Y-%m-%d').date())
out:
1999-07-01
拓展运算符
% 两个数商的余
// 两个数商的整数部分
** 幂指数
2.1.2、字符串
字符串（属于序列）**单引号，双引号，三引号（三引号应用在多行字符串）
# 1.查询方法
str.index(sub[,start[,end]])
str.find(sub[,start[,end]])
sub:指定查询的目标子串
start：指定查询的起始位置
end：指定查询的结束位置
# 字符串是否以“2018年”开头
string7 = '2017年匆匆走过，迎来崭新的2018年'
print(string7.startswith('2018年'))
# 字符串是否以“2018年”年结尾
print(string7.endswith('2018年'))
# 2.压缩方法，就是删除指定，只能首尾
str.lstrip([chars]) 左压缩
str.rstrip([chars]) 右压缩
str.strip([chars]) 左右都压缩
chars:指定待压缩的首尾字符，不指定的话，默认删除空白
# 3.替换（可以按值替换，也可以使用切片达到按位置替换的效果）
str.replace(old,new)
# 测试替换手机号中间四位为****
tel='19862516287'
print(tel.replace(tel[3:7],'****'))
# 4.格式化插入
str.format(values)
# 格式化小数，保留两位有效数字
print('{:.2f}'.format(0.12138))
out：012
# 测试生成5个有规则的文本
for month in range(5):
print('现在是{0}月{0}日'.format(str(month+1)))
# 0代表第一次使用，后面接着使用
out:
现在是1月1日
现在是2月2日
现在是3月3日
现在是4月4日
现在是5月5日
# 5.分割方法
str.split(sep) # 按照sql分割
2.1.3、列表
列表（也属于序列），与字符串有相同之处
正向单索引
list1 = ['张三','男',33,'江苏','硕士','已婚',['身高178','体重72']]
# 取出第一个元素
print(list1[0])
# 取出第四个元素
print(list1[3])
# 取出最后一个元素
print(list1[6])
# 取出“体重72”这个值
print(list1[6][1])
out:
张三
江苏
['身高178', '体重72']
体重72
负向单索引
# 取出最后一个元素
print(list1[-1])
# 取出“身高178”这个值
print(list1[-1][0])
# 取出倒数第三个元素
print(list1[-3])
out:
['身高178', '体重72']
身高178
硕士
切片索引
list2 = ['江苏','安徽','浙江','上海','山东','山西','湖南','湖北']
# 取出“浙江”至“山西”四个元素
print(list2[2:6])
# 取出“安徽”、“上海”、“山西”三个元素
print(list2[1:6:2])
# 取出最后三个元素
print(list2[-3:-1])
out:
['浙江', '上海', '山东', '山西']
['安徽', '上海', '山西']
['山西', '湖南']
无限索引
# 取出头三个元素
print(list2[:3])
# 取出最后三个元素
print(list2[-3:])
# 取出所有元素
print(list2[::])
# 取出奇数位置的元素
print(list2[::2])
# 第一个:是指从列表第一个元素开始获取，第二个:是指到最后一个元素结束
print(list2[::3])
out:
['江苏', '安徽', '浙江']
['山西', '湖南', '湖北']
['江苏', '安徽', '浙江', '上海', '山东', '山西', '湖南', '湖北']
['江苏', '浙江', '山东', '湖南']
['江苏', '上海', '湖南']
列表元素的增加append,extend,insert
list3 = [1,10,100,1000,10000]
# 在列表末尾添加数字2
list3.append(2)
print(list3)
out:
[1, 10, 100, 1000, 10000, 2]
# 在列表末尾添加20,200,2000,20000四个值
list3.extend([20,200,2000,20000])
print(list3)
out:
[1, 10, 100, 1000, 10000, 2, 20, 200, 2000, 20000]
# 在数字10后面增加11这个数字
list3.insert(2,11)
print(list3)
# 在10000后面插入['a','b','c']
list3.insert(6,['a','b','c'])
print(list3)
out:
[1, 10, 11, 100, 1000, 10000, 2, 20, 200, 2000, 20000]
[1, 10, 11, 100, 1000, 10000, ['a', 'b', 'c'], 2, 20, 200, 2000, 20000]
列表元素的删除pop,remove,clear
# 删除list3中20000这个元素
list3.pop()
print(list3)
# 删除list3中11这个元素
list3.pop(2)
print(list3)
out:
[1, 10, 11, 100, 1000, 10000, ['a', 'b', 'c'], 2, 20, 200, 2000]
[1, 10, 100, 1000, 10000, ['a', 'b', 'c'], 2, 20, 200, 2000]
# 删除list3中的['a', 'b', 'c']
list3.remove(['a', 'b', 'c'])
print(list3)
out:
[1, 10, 100, 1000, 10000, 2, 20, 200, 2000]
# 删除list3中所有元素
list3.clear()
print(list3)
out:
[]
列表元素的修改
list4 = ['洗衣机','冰响','电视机','电脑','空调']
# 将“冰响”修改为“冰箱”
print(list4[1])
list4[1] = '冰箱'
print(list4)
out:
冰响
['洗衣机', '冰箱', '电视机', '电脑', '空调']
其他方法count,index,reverse,sort
list5 = [7,3,9,11,4,6,10,3,7,4,4,3,6,3]
# 计算列表中元素3的个数
print(list5.count(3))
# 找出元素6所在的位置 ，首次出现的位置
print(list5.index(6))
# 列表元素的颠倒
list5.reverse()
print(list5)
# 列表元素的降序
list5.sort(reverse=True)
print(list5)
out:
4
5
[3, 6, 3, 4, 4, 7, 3, 10, 6, 4, 11, 9, 3, 7]
[11, 10, 9, 7, 7, 6, 6, 4, 4, 4, 3, 3, 3, 3]
2.1.4、元组
元组是不可变对象
t = ('a','d','z','a','d','c','a')
# 计数
print(t.count('a'))
# 元素位置
print(t.index('c'))
out:
3
5
2.1.5、字典
dict1 = {'姓名':'张三','年龄':33,'性别':'男','子女':{'儿子':'张四','女儿':'张美'},'兴趣':['踢球','游泳','唱歌']}
# 打印字典
print(dict1)
# 取出年龄
print(dict1['年龄'])
# 取出子女中的儿子姓名
print(dict1['子女'])
print(dict1['子女']['儿子'])
# 取出兴趣中的游泳
print(dict1['兴趣'][1])
字典元素的添加setdefault,update,键索引
# 往字典dict1中增加户籍信息
dict1.setdefault('户籍','合肥')
print(dict1)
# 增加学历信息
dict1.update({'学历':'硕士'})
print(dict1)
# 增加身高信息
dict1['身高'] = 178
print(dict1)
字典元素的删除pop,popitem,clear
# 删除字典中的户籍信息
dict1.pop('户籍')
print(dict1)
# 删除字典中女儿的姓名
dict1['子女'].pop('女儿')
print(dict1)
# 删除字典中的任意一个元素
dict1.popitem()
print(dict1)
# 清空字典元素
dict1.clear()
print(dict1)
字典元素的修改
# 将学历改为本科
dict1.update({'学历':'本科'})
print(dict1)
# 将年龄改为35
dict1['年龄'] = 35
print(dict1)
# 将兴趣中的唱歌改为跳舞
dict1['兴趣'][2] = '跳舞'
print(dict1)
其他方法
dict2 = {'电影':['三傻大闹宝莱坞','大话西游之大圣娶亲','疯狂动物城'],
'导演':['拉吉库马尔·希拉尼','刘镇伟','拜伦·霍华德 '],
'评分':[9.1,9.2,9.2]}
# 取出键'评分'所对应的值
print(dict2.get('评分'))
# 取出字典中的所有键
print(dict2.keys())
# 取出字典中的所有值
print(dict2.values())
# 取出字典中的所有键值对
print(dict2.items())
out
[9.1, 9.2, 9.2]
dict_keys(['电影', '导演', '评分'])
dict_values([['三傻大闹宝莱坞', '大话西游之大圣娶亲', '疯狂动物城'], ['拉吉库马尔·希拉尼', '刘镇伟', '拜伦·霍华德 '], [9.1, 9.2, 9.2]])
dict_items([('电影', ['三傻大闹宝莱坞', '大话西游之大圣娶亲', '疯狂动物城']), ('导演', ['拉吉库马尔·希拉尼', '刘镇伟', '拜伦·霍华德 ']), ('评分', [9.1, 9.2, 9.2])])
2.2、控制流循环
# if elif else
# 输入18位身份证号码，识别性别
ID=input('输入18位身份证号码')
if int(ID)[-2]%2==0:
print('女')
else:
print('男')
# for
# 对列表中的偶数作三次方减10的处理【使用列表推导式】
list7 = [3,1,18,13,22,17,23,14,19,28,16]
result = [i ** 3 - 10 for i in list7 if i % 2 == 0]
print(result)
out:
[5822, 10638, 2734, 21942, 4086]
# while
# while更适合无具体迭代对象的循环
2.3、正则表达式
参考菜鸟教程：https://www.runoob.com/regexp/regexp-syntax.html
菜鸟教程的正则表达式在线测试以及常用正则表达式：https://c.runoob.com/front-end/854/
注意：.*?
**正则表达式：**描述或刻画字符串内在规律的表达式
使用场景：
- 无法使用切片
- 借助replace无法完成非固定值的替换
- 借助split无法按照多种值实现切割
import re
# 1.查询
findall(pattern, string, flags=0)
# pattern:正则表达式
# string:字符串
# flags:匹配模式，一般不用
# 2.替换
sub(pattern,repl,string,count=0,flags=0)
# repl:指定替换成的新值
# count:用于指定最多替换的次数，默认全部替换，一般不用
# 3.分割
split(pattern,string,maxsplit=0,flags=0)
# maxsplit:指定切割次数，一般不用
测试：
# 导入第三方包
import re
# 1.取出出字符中所有的天气状态
string8 = "{ymd:'2018-01-01',tianqi:'晴',aqiInfo:'轻度污染'},{ymd:'2018-01-02',tianqi:'阴~小雨',aqiInfo:'优'},{ymd:'2018-01-03',tianqi:'小雨~中雨',aqiInfo:'优'},{ymd:'2018-01-04',tianqi:'中雨~小雨',aqiInfo:'优'}"
print(re.findall("tianqi:'(.*?)'", string8))
out:
['晴', '阴~小雨', '小雨~中雨', '中雨~小雨']
# 2.取出所有含O字母的单词
string9 = 'Together, we discovered that a free market only thrives when there are rules to ensure competition and fair play, Our celebration of initiative and enterprise'
print(re.findall('\w*o\w*',string9, flags = re.I))
out:
['Together', 'discovered', 'only', 'to', 'competition', 'Our', 'celebration', 'of']
# 3.将标点符号、数字和字母删除
string10 = '据悉，这次发运的4台蒸汽冷凝罐属于国际热核聚变实验堆（ITER）项目的核二级压力设备，先后完成了压力试验、真空试验、氦气检漏试验、千斤顶试验、吊耳载荷试验、叠装试验等验收试验。'
print(re.sub('[，。、a-zA-Z0-9（）]','',string10))
out:
据悉这次发运的台蒸汽冷凝罐属于国际热核聚变实验堆项目的核二级压力设备先后完成了压力试验真空试验氦气检漏试验千斤顶试验吊耳载荷试验叠装试验等验收试验
# 4.将每一部分的内容分割开
string11 = '2室2厅 | 101.62平 | 低区/7层 | 朝南 \n 上海未来 - 浦东 - 金杨 - 2005年建'
split = re.split('[-\|\n]', string11)
print(split)
split_strip = [i.strip() for i in split]
print(split_strip)
out:
['2室2厅 ', ' 101.62平 ', ' 低区/7层 ', ' 朝南 ', ' 上海未来 ', ' 浦东 ', ' 金杨 ', ' 2005年建']
['2室2厅', '101.62平', '低区/7层', '朝南', '上海未来', '浦东', '金杨', '2005年建']
2.4、Python函数
2.4.1、匿名函数lambda
lambda parameters : fun
# parameters 形参
# fun 函数体
测试
# 统计列表中每个元素的频次
list6 = ['A','A','B','A','A','B','C','B','C','B','B','D','C']
# 构建空字典，用于频次统计数据的存储
dict3 = {}
# 循环计算
for i in set(list6):
dict3[i] = list6.count(i)
print(dict3)
# 取出字典中的键值对
key_value = list(dict3.items())
print(key_value)
# 列表排序
key_value.sort()
print(key_value)
# 按频次高低排序
key_value.sort(key = lambda x : x[1], reverse=True)
print(key_value)
out:
{'B': 5, 'D': 1, 'C': 3, 'A': 4}
[('B', 5), ('D', 1), ('C', 3), ('A', 4)]
[('A', 4), ('B', 5), ('C', 3), ('D', 1)]
[('B', 5), ('A', 4), ('C', 3), ('D', 1)]
2.4.2、自定义函数
# 1.猜数字
def game(min,max):
import random
number = random.randint(min,max) # 随机生成一个需要猜的数字
while True:
guess = float(input('请在%d到%d之间猜一个数字: ' %(min, max)))
if guess < number:
min = guess
print('不好意思，你猜的的数偏小了！请在%d到%d之间猜一个数！' %(min,max))
elif guess > number:
max = guess
print('不好意思，你猜的的数偏大了！请在%d到%d之间猜一个数！' %(min,max))
else:
print('恭喜你猜对了！')
print('游戏结束！')
break
# 2.计算1到n的平方和
def square_sum(n, p = 2):
result = sum([i ** p for i in range(1,n+1)])
return(n,p,result)
print('1到%d的%d次方和为%d！' %square_sum(200))
print('1到%d的%d次方和为%d！' %square_sum(200,3))
# 3.任意个数的数据求和
def adds(*args):
print(args)
s = sum(args)
return(s)
print('和为%d!' %adds(10,13,7,8,2))
print('和为%d!' %adds(7,10,23,44,65,12,17))
# 4.关键字参数
def info_collection(tel, birthday, **kwargs):
user_info = {} # 构造空字典，用于存储用户信息
user_info['tel'] = tel
user_info['birthday'] = birthday
user_info.update(kwargs)
# 用户信息返回
return(user_info)
# 调用函数
info_collection(13612345678,'1990-01-01',nickname='月亮',gender = '女',edu = '硕士',income = 15000,add = '上海市浦东新区',interest = ['游泳','唱歌','看电影'])
# 指定文件夹下数据的读取和合并
import pandas as pd
import os
Path = r'D:\study\pylearn\datas'
# 'r'是防止字符转义的， 如果字符串中出现'\n'的话 ，不加r的话，\n就会被转义成换行符,而加了'r'之后'\n'就能保留原有的样子
def readdatas(Path):
filenames = os.listdir(Path)
res = []
for file in filenames:
data = pd.read_csv(Path+'\\'+file)
res.append(data)
return pd.concat(res)
# 调用函数
readdatas(Path)
3、网络爬虫
三个案例，前两个（红牛，链家）都是静态网页，第三个（2345天气）为动态网页。
所需要的库
import requests # 发送请求
import re # 正则表达式
import bs4 # 解析
request.get # 基于URL，发送网络请求
re.findall # 正则表达式
bs4.BeautifulSoup # 对HTML源代码做解析，便于目标数据的拆解
3.1、爬取红牛分公司
爬取时间2021.12.5
网页：http://www.redbull.com.cn/about/branch
爬取网页的这些数据【公司，地区，邮编，电话】
# 网页：http://www.redbull.com.cn/about/branch
import requests
import re
url=r'http://www.redbull.com.cn/about/branch'
response = requests.get(url)
response
# 输出：<Response [200]>
# 字符串类型的html
response.text
# 输出：'<!DOCTYPE html><html lang="en"><head>....'
# 1.使用正则表达式获取值
# 分析源代码找到公司名
re.findall('<h2>红牛东莞分公司</h2>',response.text)
# 输出：['<h2>红牛东莞分公司</h2>']
re.findall('<h2>.*?</h2>',response.text)
# 输出：['<h2>红牛杭州分公司</h2>',
# '<h2>红牛广西分公司</h2>',
# '<h2>红牛广州分公司</h2>',
# '<h2>红牛深圳分公司</h2>',
# '<h2>红牛湖南分公司</h2>',
# ...]
company = re.findall('<h2>(.*?)</h2>',response.text)
company
# 输出：['红牛杭州分公司',
# '红牛广西分公司',
# '红牛广州分公司',
# '红牛深圳分公司',
# '红牛湖南分公司',
# '红牛福建分公司',
# ...]
re.findall("<p class=\'mapIco\'>.*?</p>",response.text)
#输出：["<p class='mapIco'>杭州市上城区庆春路29号远洋大厦11楼A座</p>",
# "<p class='mapIco'>南宁市金湖路59号地王国际商会中心50层D1、E1室</p>",
# "<p class='mapIco'>广东省广州市天河珠江新城华夏路10号富力中心写字楼1904房</p>",
# "<p class='mapIco'>广东省深圳市福田区福华三路88号财富大厦39楼BCD</p>",
# "<p class='mapIco'>湖南省长沙市天心区劳动西路289号嘉盛国际广场1626室</p>",
#...]
add = re.findall("<p class=\'mapIco\'>(.*?)</p>",response.text)
add
# 输出：['杭州市上城区庆春路29号远洋大厦11楼A座',
# '南宁市金湖路59号地王国际商会中心50层D1、E1室',
# '广东省广州市天河珠江新城华夏路10号富力中心写字楼1904房',
# '广东省深圳市福田区福华三路88号财富大厦39楼BCD',
# '湖南省长沙市天心区劳动西路289号嘉盛国际广场1626室',
# ...
# 2.下面使用bs4直接获取值
import bs4
soup = bs4.BeautifulSoup(response.text)
# for i in soup.findAll(name='p',attrs={'class':'mailIco'}):
# print(i.text)
# 列表表达式，与上面for循环是一样的
mall = [i.text for i in soup.findAll(name='p',attrs={'class':'mailIco'})]
mall
# 输出：['310009',
# '530021',
# '510623',
# '518048',
# ...]
tel = [i.text for i in soup.findAll(name='p',attrs={'class':'telIco'})]
tel
# 输出：['0571-87045279/7792',
# '0771-5592660/61/62',
# '020-38927681',
# '0755-23962001',
# ...]
import pandas as pd
pd.DataFrame({'company':company,'add':add,'mall':mall,'tel':tel})
最后的结果：
下图为原网站的部分代码截图：
3.2、爬取链家二手房
爬取时间2021.12.5
网页：https://sh.lianjia.com/ershoufang/pudong/pg1/
爬取网页上这些信息：小区，（户型，面积，朝向，装修，楼层，建筑年代，建筑样式），总售价
import requests
import re
import bs4
url=r'https://sh.lianjia.com/ershoufang/pudong/pg1/'
# 这里我直接请求，返回的也是200，但老师演示的是403，猜测可能是链家的反扒失效了？
response = requests.get(url)
response
# 输出：<Response [200]>
# 接下来按照请求403来解决反扒问题，模拟浏览器访问
F12，然后F5刷新，找到网络network，随便点击一个请求信息，点击表头headers，然后在下面请求标头Request Headers找到User-Agent，这就是浏览器的基本信息
# 模拟浏览器访问
Headers = {
'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/96.0.4664.45 Safari/537.36'
}
# response = requests.get(url) # 可能会出现403问题，反扒
response = requests.get(url,headers = Headers)
response
# 输出：<Response [200]>
soup = bs4.BeautifulSoup(response.text)
xaioqu = [i.text.strip() for i in soup.findAll(name='a',attrs={'data-el':'region'})]
xaioqu
# 输出：['东升家园(三期)',
# '齐爱佳苑',
# ...]
list = [i.text.split('|') for i in soup.findAll(name='div',attrs={'class':'houseInfo'})]
huxing = [i[0].strip() for i in list]
mianji = [i[1].strip() for i in list]
chaoxiang = [i[2].strip() for i in list]
zhuangxiu = [i[3].strip() for i in list]
louceng = [i[4].strip() for i in list]
jianzhuniandai = [i[5].strip() for i in list]
jianzhuyangshi = [i[-1].strip() for i in list]
# 输出：[['1室2厅 ', ' 79.32平米 ', ' 南 北 ', ' 毛坯 ', ' 低楼层(共15层) ', ' 2006年建 ', ' 板楼'],
# ['2室2厅 ', ' 82.6平米 ', ' 南 北 ', ' 精装 ', ' 中楼层(共6层) ', ' 2008年建 ', ' 板楼'],
# ...]
price = [i.text.strip() for i in soup.findAll(name='div',attrs={'class':'totalPrice totalPrice2'})]
price
# 输出：['218万',
# '449万',
# ...]
danjia = [i.text for i in soup.findAll(name='div',attrs={'class':'unitPrice'})]
danjia
# 输出：['27,484元/平',
# '54,359元/平',
# ...]
import pandas as pd
pd.DataFrame({'小区':xaioqu,'户型':huxing,'面积':mianji,'朝向':chaoxiang,'装修':zhuangxiu,'楼层':louceng,'建筑年代':jianzhuniandai,'建筑样式':jianzhuyangshi,'总价':price,'单价':danjia})
最后的dataframe
3.3、爬取2345天气
爬取时间2021.12.5
之前的红牛和链家都是静态网页，2345天气为动态网页
网页：http://tianqi.2345.com/wea_history/58362.htm
我们需要爬取【日期，最高温，最低温，天气，锋利风向，空气质量指数】
# 教程上的例子直接找到了原始js数据，但是这种请求在2021年被取消了
# http://tianqi.2345.com/t/wea_history/js/58362_{年}{月}.js
#例： http://tianqi.2345.com/t/wea_history/js/58362_201710.js
# 2021年以后的请求都是这种了
# http://tianqi.2345.com/Pc/GetHistory?areaInfo%5BareaId%5D=58362&areaInfo%5BareaType%5D=2&date%5Byear%5D={年}&date%5Bmonth%5D={月}
#例： http://tianqi.2345.com/Pc/GetHistory?areaInfo%5BareaId%5D=58362&areaInfo%5BareaType%5D=2&date%5Byear%5D=2021&date%5Bmonth%5D=12
# 先按照之前那的请求做
之前的请求是这样的js
直接使用正则表达式去切分
import requests
import re
import bs4
url =r'http://tianqi.2345.com/t/wea_history/js/58362_201710.js'
response = requests.get(url)
# 要加这一行，不然会乱码
response.encoding = 'gbk'
response.text
乱码：
修改后：
乱码问题解决后，直接正则表达式切分
date = re.findall("ymd:'([\d-]+)',",response.text)
print(date)
high = re.findall("bWendu:'([\d-]+)℃',",response.text)
print(high)
low = re.findall("yWendu:'([\d-]+)℃',",response.text)
print(low)
weather = re.findall("tianqi:'(.*?)',",response.text)
print(weather)
wind = re.findall("fengxiang:'(.*?)',",response.text)
print(wind)
api = re.findall("aqiInfo:'(.*?)',",response.text)
print(api)
# 最后转成dataframe
pd.DataFrame({'日期':date,'最高温':high,'最低温':low,'天气':weather,'风力风向':wind,'空气质量':api})
上面是按照之前的请求去做的，但是2021年2345天气官方废除了这种请求，所以新的数据获取不到。使用新的请求爬取
url =r'http://tianqi.2345.com/Pc/GetHistory?areaInfo%5BareaId%5D=58362&areaInfo%5BareaType%5D=2&date%5Byear%5D=2021&date%5Bmonth%5D=11'
response = requests.get(url)
response.text
# <td>2021-11-01
date = re.findall("<td>([\d-]+)",response.text)
print(date)
# <td style=\\"color:#ff5040;\\">22
high = re.findall('<td style=\\\\"color:#ff5040;\\\\">([\d]+)',response.text)
print(high)
# <td style=\\"color:#3097fd;\\" >8
low = re.findall('<td style=\\\\"color:#3097fd;\\\\" >([\d]+)',response.text)
print(low)
# <td>\\u591a\\u4e91~\\u6674<\\/td>
# 实在是没什么规律，只能从波浪线展开
list = re.findall('<td>([\\\\[a-zA-Z0-9_]+]*?~*[\\\\[a-zA-Z0-9_]+]*?)<\\\\/td>',response.text)
print(list)
# 终于搞定了，然后把这个uft8编码转为汉子就行了
list = [i.encode('utf-8').decode("unicode_escape") for i in list]
print(list)
weather = []
wind = []
for i in range(len(list)):
if i%2==0:
weather.append(list[i])
else:
wind.append(list[i])
print(weather)
print(wind)
# <td><span class=\\"history-aqi wea-aqi-2\\">67 \\u826f<\\/span><\\/td>
# <td><span class=\\"history-aqi wea-aqi-1\\">47 \\u4f18<\\/span><\\/td>
api = re.findall('<td><span class=\\\\"history-aqi wea-aqi-[0-9]\\\\">[0-9]+ (\\\\[a-zA-Z0-9_]+)<\\\\/span><\\\\/td>',response.text)
print(api)
api = [i.encode('utf-8').decode("unicode_escape") for i in api]
print(api)
print(len(date))
print(len(high))
print(len(low))
print(len(weather))
print(len(wind))
print(len(api))
# 最后转成dataframe
pd.DataFrame({'日期':date,'最高温':high,'最低温':low,'天气':weather,'风力风向':wind,'空气质量':api})
4、Numpy
数组的构造及优势
常用的数学函数与统计函数
4.1、构造numpy数组使用
生成默认数组
# 怎样构造numpy数组
import numpy as np
# 生成默认数组
a = np.arange(3) #一个参数 默认起点0，步长为1 输出：[0 1 2]
a = np.arange(3,9) #两个参数 默认步长为1 输出[3 4 5 6 7 8]
a = np.arange(0, 3, 0.1) #三个参数 起点为0，终点为3，步长为0.1 输出[ 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2. 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9]
a = np.arange(12).reshape(3,4) # [[ 0 1 2 3][ 4 5 6 7][ 8 9 10 11]]
已有数据转数组
# 已有数据
# 一维数组：
# 列表（元组同）
list = [3,10,8,7,34,11,28,72]
arr = np.array(list) # list转numpy.ndarray
print('一维数组：\n',arr)
# 一维数组元素的获取
arr[0] # 3
arr[[0]] # array([3])
arr[[1,5,7]] # array([10, 11, 72]) # 取多个
arr[arr<10] # array([3, 8, 7]) # 逻辑索引
# 二维数组：
# 嵌套元组（列表同）
tuple = ((8.5,6,4.1,2,0.7),(1.5,3,5.4,7.3,9),(3.2,3,3.8,3,3),(11.2,13.4,15.6,17.8,19))
arr = np.array(tuple) # tuple转numpy.ndarray
print('二维数组：\n',arr)
# 二维数组元素的获取
arr[1,2] # 第2行第3列元素
arr[2,:] # 第3行所有元素
arr[:,1] # 第2列所有元素
arr[1:4,1:5] # 第2至4行，2至5行
# ×第1行、最后一行和第2列、第4列构成的数组(被解释为是[0,1]和[-1,3],这不是我们想要的结果)
# arr[[0,-1],[1,3]] # [ 6. , 17.8]
# 第1行、最后一行和第2列、第4列构成的数组(这个是两行两列交叉的四个点)
arr2[np.ix_([0,-1],[1,3])] # [[ 6. , 2. ][13.4, 17.8]]
# 第1行、最后一行和第1列、第3列、第4列构成的数组(两行三列交叉的六个点)
arr[np.ix_([0,-1],[1,2,3])] # [[ 6. 4.1 2. ][13.4 15.6 17.8]]
查看数组结构
# 查看数组结构
type(arr) # numpy.ndarray
arr.ndim # 2 查看数组的维数
arr.shape # (4,5) 查看数组的行列数
arr.dtype # dtype('float64') 查看数组的数据类型
arr.size # 20 查看数组元素总个数
4.2、通过计算了解numpy效率
# 计算 身体质量指数 bmi=w/(h/100)**2
# 随机生成1千万个身高体重
import numpy as np
import random
h = []
w = []
for i in range(10000000):
h.append(random.randint(153,180))
w.append(random.uniform(51,88))
# 使用for遍历循环运算
%%time
bmi = []
for i in range(10000000):
bmi.append(w[i]/(h[i]/100)**2)
bmi[:5]
# 输出：Wall time: 4.36 s
# 使用numpy数组运算
%%time
H = np.array(h)
W = np.array(w)
BMI = W/(H/100)**2
BMI[:5]
# 输出：Wall time: 1.28 s
# 可以发现numpy的效率比for循环遍历高
4.3、numpy数学与统计函数
常见的数学函数
np.round(arr) # 各元素的四舍五入
np.sqrt(arr) # 各元素的算数平方根
np.square(arr) # 各元素的平方
np.exp(arr) # 计算以e为底的指数(比如np.exp(5)=e的5次方)
np.power(arr,α) # 计算各元素的指数(arr元素的五次方就是α=5)
np.log(arr) # 计算以e为底各元素的对数 # np.log2(arr) # 计算以2为底各元素的对数 # np.log10(arr) # 计算以10为底各元素的对数
arr7 = np.array([[1,2,10],[10,8,3],[7,6,5]])
arr8 = np.array([[2,2,2],[3,3,3],[4,4,4]])
# 取子集
# 从arr7中取出arr7大于arr8的所有元素
print('满足条件的二维数组元素获取：\n',arr7[arr7>arr8])
# 从arr9中取出大于10的元素
arr9 = np.array([3,10,23,7,16,9,17,22,4,8,15])
print('满足条件的一维数组元素获取：\n',arr9[arr9>10])
# 判断操作
# 将arr7中大于7的元素改成5，其余的不变
print('二维数组的条件操作：\n',np.where(arr7>7,5,arr7))
# 将arr9中大于10 的元素改为1，否则改为0
print('一维数组的条件操作：\n',np.where(arr9>10,1,0))
常见的统计函数
# axis两种取值，0或1，当axis=0,则纵向计算,当axis=1,则横向计算
np.min(arr,axis) # 按照轴的方向计算最小值
np.max(arr,axis) # 按照轴的方向计算最大值
np.mean(arr,axis) # 按照轴的方向计算平均值
np.average(arr,axis) # 同mean,average可以做加权平均值
np.median(arr,axis) # 按照轴的方向计算中位数
np.sum(arr,axis) # 按照轴的方向计算和
np.std(arr,axis) # 按照轴的方向计算标准差
np.var(arr,axis) # 按照轴的方向计算方差
mean和average的区别
np.mean()和np.average()都是计算均值。
不加权时，np.mean()和np.average()都一样。
np.average()可以计算加权平均。
#
a = np.array([1, 2, 3, 4, 5])
aw = np.array([0.1, 0.2, 0.3, 0.4, 0.5])
print('平均：', np.mean(a)) # 3.0
print('平均：', np.average(a)) # 3.0
print('加权平均：', np.average(a,weights = aw)) # 3.666
np.sum(arr,axis=1)和arr.sum(axis=1)是等价的
意思就是，np.fun 和arr.fun其实是等价的（包的函数=对象的方法）
广播运算
规则：
- 各输入数组的维度可以不相等，但必须确保从右到左的对应维度值相等
- 如果对应维度值不相等，就必须保证其中一个为1
- 各输入数组都向其shape最长的数组看齐，shape中不足的部分都通过在前面加1补齐
# 各输入数组维度一致，对应维度值相等
arr10 = np.arange(12).reshape(3,4)
arr11 = np.arange(101,113).reshape(3,4)
print('3×4的二维矩阵运算：\n',arr10 + arr11)
# 各输入数组维度不一致，对应维度值相等
arr12 = np.arange(60).reshape(5,4,3)
arr10 = np.arange(12).reshape(4,3)
print('维数不一致，但末尾的维度值一致：\n',arr12 + arr10)
# 各输入数组维度不一致，对应维度值不相等，但其中有一个为1
arr12 = np.arange(60).reshape(5,4,3)
arr13 = np.arange(4).reshape(4,1)
print('维数不一致，维度值也不一致，但维度值至少一个为1：\n',arr12 + arr13)
# 加1补齐
arr14 = np.array([5,15,25])
print('arr14的维度自动补齐为(1,3)：\n',arr10 + arr14)
4.4、numpy的随机数
之前使用的random一次只能生成一个随机数，我们需要大量随机数时，只能循环
而numpy提供了random方法
np.random.randint()生成的随机整数的取值区间是**前闭后开[1,10)区间；
random.randint()生成的随机整数的取值区间是前闭后闭[1,10]**区间。
import numpy as np
# 随机整数
np.random.randint(0,1,20) # 从[0,1)之间生成20个随机数
# 随机均匀分布
np.random.uniform(0,1,20)
# # 随机正态分布
np.random.normal(0,1,20)
测试案例：
# 理想赌场：输赢概率皆为0.5,输-8,赢+8,初始金额1000,赌1500次
# 随机生成1500个数
R = np.random.uniform(0,1,1500)
money = 1000
packages = [money]
for i in R:
if i < 0.5:
money-=8
else:
money+=8
packages.append(money)
# 作图
import matplotlib.pyplot as plt
plt.plot(range(1501),packages)
plt.show()
测试：
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats
# 生成各种正态分布随机数
np.random.seed(1234)
rn1 = np.random.normal(loc = 0, scale = 1, size = 1000)
rn2 = np.random.normal(loc = 0, scale = 2, size = 1000)
rn3 = np.random.normal(loc = 2, scale = 3, size = 1000)
rn4 = np.random.normal(loc = 5, scale = 3, size = 1000)
# 绘图
plt.style.use('ggplot')
sns.distplot(rn1, hist = False, kde = False, fit = stats.norm,
fit_kws = {'color':'black','label':'u=0,s=1','linestyle':'-'})
sns.distplot(rn2, hist = False, kde = False, fit = stats.norm,
fit_kws = {'color':'red','label':'u=0,s=2','linestyle':'--'})
sns.distplot(rn3, hist = False, kde = False, fit = stats.norm,
fit_kws = {'color':'blue','label':'u=2,s=3','linestyle':':'})
sns.distplot(rn4, hist = False, kde = False, fit = stats.norm,
fit_kws = {'color':'purple','label':'u=5,s=3','linestyle':'-.'})
# 呈现图例
plt.legend()
# 呈现图形
plt.show()
# 生成各种指数分布随机数
np.random.seed(1234)
re1 = np.random.exponential(scale = 0.5, size = 1000)
re2 = np.random.exponential(scale = 1, size = 1000)
re3 = np.random.exponential(scale = 1.5, size = 1000)
# 绘图
sns.distplot(re1, hist = False, kde = False, fit = stats.expon,
fit_kws = {'color':'black','label':'lambda=0.5','linestyle':'-'})
sns.distplot(re2, hist = False, kde = False, fit = stats.expon,
fit_kws = {'color':'red','label':'lambda=1','linestyle':'--'})
sns.distplot(re3, hist = False, kde = False, fit = stats.expon,
fit_kws = {'color':'blue','label':'lambda=1.5','linestyle':':'})
# 呈现图例
plt.legend()
# 呈现图形
plt.show()
4.5、线性代数的相关计算
使用numpy的linalg模块
参考菜鸟教程：https://www.runoob.com/numpy/numpy-linear-algebra.html
np.zeros() # 生成零矩阵
np.eye() # 生成单位矩阵
np.dot() # 计算两个数组的点积
np.diag() # 矩阵主对角线与一维数组间的转换
np.linalg.det() # 计算矩阵行列式
np.linalg.eigvals() # 计算方阵特征根
np.linalg.pinv() # 计算方阵的Moore-Penrose伪逆
np.linalg.lstsq # 计算Ax=b的最小二乘解
np.linalg.svd() # 计算奇异值分解
np.ones() # 生成所有元素为1的矩阵
np.transpose() # 矩阵转置
np.inner() # 计算两个数组的内积
np.trace() # 矩阵主对角线元素的和
np.linalg.eig() # 计算矩阵特征根与特征向量
np.linalg.inv() # 计算方阵的逆
np.linalg.solve() # 计算Ax=b的线性方程组的解
np.linalg.qr() # 计算QR分解
np.linalg.norm() # 计算向量或矩阵的范数
测试：
矩阵乘法：
# 一维数组的点积
vector_dot = np.dot(np.array([1,2,3]), np.array([4,5,6]))
print('一维数组的点积：\n',vector_dot)
# 二维数组的乘法
print('两个二维数组：')
print(arr10)
print(arr11)
arr2d = np.dot(arr10,arr11)
print('二维数组的乘法：\n',arr2d)
diag函数的使用
# diag的使用
arr15 = np.arange(16).reshape(4,-1)
print('4×4的矩阵：\n',arr15)
print('取出矩阵的主对角线元素：\n',np.diag(arr15))
print('由一维数组构造的方阵：\n',np.diag(np.array([5,15,25])))
特征根与特征向量
# 计算方阵的特征向量和特征根
arr16 = np.array([[1,2,5],[3,6,8],[4,7,9]])
print('计算3×3方阵的特征根和特征向量：\n',arr16)
print('求解结果为：\n',np.linalg.eig(arr16))
多元线性回归模型的解
# 计算偏回归系数
X = np.array([[1,1,4,3],[1,2,7,6],[1,2,6,6],[1,3,8,7],[1,2,5,8],[1,3,7,5],[1,6,10,12],[1,5,7,7],[1,6,3,4],[1,5,7,8]])
Y = np.array([3.2,3.8,3.7,4.3,4.4,5.2,6.7,4.8,4.2,5.1])
X_trans_X_inverse = np.linalg.inv(np.dot(np.transpose(X),X))
beta = np.dot(np.dot(X_trans_X_inverse,np.transpose(X)),Y)
print('偏回归系数为：\n',beta)
多元一次方程组的求解
# 多元线性方程组
A = np.array([[3,2,1],[2,3,1],[1,2,3]])
b = np.array([39,34,26])
X = np.linalg.solve(A,b)
print('三元一次方程组的解：\n',X)
范数的计算
# 范数的计算
arr17 = np.array([1,3,5,7,9,10,-12])
# 一范数
res1 = np.linalg.norm(arr17, ord = 1)
print('向量的一范数：\n',res1)
# 二范数
res2 = np.linalg.norm(arr17, ord = 2)
print('向量的二范数：\n',res2)
# 无穷范数
res3 = np.linalg.norm(arr17, ord = np.inf)
print('向量的无穷范数：\n',res3)
5、Pandas
外部数据的读取
数据概览
子集筛选
数据汇总
数据合并连接
5.1、序列Series与数据框Dataframe
序列Series（序列与一维数组很相似）
如果要做数学运算，使用numpy，（pandas在这方面比较匮乏）
如果对序列做统计运算，则都可以用，用pandas的话方法更为丰富（如偏度，峰度等这些numpy没有）
# 导入模块
import pandas as pd
import numpy as np
# 构造序列
gdp1 = pd.Series([2.8,3.01,8.99,8.59,5.18])
gdp2 = pd.Series({'北京':2.8,'上海':3.01,'广东':8.99,'江苏':8.59,'浙江':5.18})
gdp3 = pd.Series(np.array((2.8,3.01,8.99,8.59,5.18)))
print(type(gdp1)) # <class 'pandas.core.series.Series'>
print(gdp1)
print(gdp2)
print(gdp3)
# 取出gdp1中的第一、第四和第五个元素
print('行号风格的序列：\n',gdp1[[0,3,4]])
# 取出gdp2中的第一、第四和第五个元素
print('行名称风格的序列：\n',gdp2[[0,3,4]])
# 取出gdp2中上海、江苏和浙江的GDP值
print('行名称风格的序列：\n',gdp2[['上海','江苏','浙江']])
# 数学函数--取对数
print('通过numpy函数：\n',np.log(gdp1))
# 平均gdp
print('通过numpy函数：\n',np.mean(gdp1))
print('通过序列的方法：\n',gdp1.mean())
数据框DataFrame
# 构造数据框
df1 = pd.DataFrame([['张三',23,'男'],['李四',27,'女'],['王二',26,'女']])
df2 = pd.DataFrame({'姓名':['张三','李四','王二'],'年龄':[23,27,26],'性别':['男','女','女']})
df3 = pd.DataFrame(np.array([['张三',23,'男'],['李四',27,'女'],['王二',26,'女']]))
print('嵌套列表构造数据框：\n',df1)
print('字典构造数据框：\n',df2)
print('二维数组构造数据框：\n',df3)
5.2、外部数据读取
文本文件、电子表格、数据库
（1）文本文件的读取txt,csv
read_csv或者read_table
import pandas as pd
# 读取文本文件，读取csv文件
pd.read_csv(filepath_or_buffer,sep=',',header=None,names=None,usecols=None,
skiprows=None,skipfooter=None,converters=None,encoding=None)
filepath_or_buffer # 指定txt或csv文件所在的路径
sep # 指定源数据集中各字段之间的分隔符，默认为逗号","
header # 是否将第一行作为表头，默认将第一行作为字段名称(没有表头则header=None)
names # 添加表头
usecols # 读取指定字段
skiprows # 跳过首行数
skipfooter # 跳过尾行数
converters # 转换数据类型
encoding # 编码
comment # 指定注释符，如果行首碰到注释符，则跳过该行
parse_dates
# 1.如果参数值为列表,解析对应的日期列(嵌套列表的话,合并列为日期列)
# 2.如果参数值为True,尝试解析数据框的行索引
# 3.如果参数值为字典,解析对应的列，生成新的字段名
测试读取D:\study\pylearn\pandasdata\data_test01.txt数据
import pandas as pd
data01 = pd.read_csv(r'D:\study\pylearn\pandasdata\data_test01.txt',
sep=',',skiprows=2,skipfooter=3,comment='#',converters={'id':str},
parse_dates={'birthday':[1,2,3]},thousands='&',encoding='utf-8')
# 如果没有header(字段名)的话,可以这样
# header=None,names=['生日','编号','性别','职位','收入']
# 想要读取指定列，可以用
# usecols=['id','income']
print(type(data01)) # <class 'pandas.core.frame.DataFrame'>
data01
（2）读取电子表格Excel
read_excel
import pandas as pd
# 读取excel
pd.read_excel(io,sheet_name=0,header=0,skiprows=None,skipfooter=0,
index_col=None,names=None,na_values=None,thousands=None,convert_float=True)
io # 指定电子表格的具体路径
sheet_name # 指定需要读取的sheet(可以传递整数,也可以传递具体的sheet名称)
header # 是否需要将数据集的第一行用作表头,默认需要
skiprows # 跳过首行数
skipfooter # 跳过尾行数
index_col # 指定那些列用作数据框的行索引(标签)
测试读取D:\study\pylearn\pandasdata\data_test02.xlsx
pd.read_excel(r'D:\study\pylearn\pandasdata\data_test02.xlsx',header = None,
names = ['Prod_Id','Prod_Name','Prod_Color','Prod_Price'],converters={'ID':str})
（3）读取数据库（+操作）
参考菜鸟教程：https://www.runoob.com/python3/python3-mysql.html
pymysql.connect()
连接建表
import pandas as pd
import pymysql
db = pymysql.connect(host="localhost",user="root",password="123456",
database="pandas211206",port=3306,charset="utf8")
# 使用cursor()方法获取操作游标
cursor = db.cursor()
# 建表
sql1 = '''create table `user` (
`id` int,
`name` varchar(30),
`sex` varchar(2))'''
cursor.execute(sql1)
插入insert
# 插入
sql2 = """insert into `user` (`id`,`name`,`sex`)
values
(1,'张三','男'),
(2,'李四','女'),
(3,'王五','男'),
(4,'岩六','男')"""
try:
cursor.execute(sql2) # 执行sql语句
db.commit() # 提交到数据库执行
except:
db.rollback() # 如果发生错误则回滚
修改update
# 修改
sql3 = """update `user` set name = '张三12138' where `id` = 1"""
try:
cursor.execute(sql3) # 执行sql语句
db.commit() # 提交到数据库执行
except:
db.rollback() # 如果发生错误则回滚
删除delete
# 删除
sql4 = "delete from `user` where `id` = 3"
try:
cursor.execute(sql4) # 执行sql语句
db.commit() # 提交到数据库执行
except:
db.rollback() # 如果发生错误则回滚
查询select（pymysql原生）
# pymysql的原生查询(不建议使用，直接用下一个pd.read_sql)
# SQL 查询语句
sql = "select * from user"
try:
# 执行SQL语句
cursor.execute(sql)
# 获取所有记录列表
results = cursor.fetchall()
for row in results:
id = row[0]
name = row[1]
sex = row[2]
# 打印结果
print ("id=%s,name=%s,sex=%s" % \
(id, name, sex))
except:
print ("Error: 未找到数据")
查询select（pd.read_sql）
# 查找pd.read_sql
user = pd.read_sql('select * from user', db)
# 数据输出
user
关闭数据库
# 关闭数据库连接
db.close()
5.3、查看数据信息
# 查看user的一些内容信息
print(user.shape)
print(user.head)
print(user.columns)
print(user.dtypes)
print(user.describe(include = 'object'))
# 对列的筛选
print(user.name) # 不建议这样取数据
print(user['name']) # 建议这样取数据
# 对行的筛选
user[user['sex']=='男']
user.loc[(user.sex=='男') & (user['id']==1),['name','sex']]
5.4、数据的清洗（类型，冗余，异常，缺失）
# 1.数据类型的修改
pd.to_datetime
df.column.astype
# 2.冗余数据的识别与处理
df.duplicated
df.drop.duplicates
# 3.异常点
# z得分法，分位数法，距离法
# 4.缺失值
df.isnull
df.fillna
df.dropna
(1) 数据类型的修改
读取D:\study\pylearn\pandasdata\sec_cars.csv文件进行操作
import pandas as pd
sec_car = pd.read_csv(r'D:\study\pylearn\pandasdata\sec_cars.csv')
sec_car.head()
sec_car.dtypes
# 把日期Boarding_time的object转换为datatime类型
sec_car.Boarding_time = pd.to_datetime(sec_car.Boarding_time,format='%Y年%m月')
# 把New_price的object转换为float类型
sec_car.New_price = sec_car.New_price.str[:-1].astype(float) # 坑：这里需要加一个.str
sec_car.dtypes
# 数据的描述性统计
sec_cars.describe()
# 我们可以发现,这里面包括[非缺失值个数,平均值,标准差,最小值,下四分位数,中位数,上四分位数,最大值]
# 以售价Sec_price为例,平均价格25.65万,中位数为10.20万,上四分位数(75%)为23.80万,最高价为808万,很明显，平均值受到了极端高价的影响
然后我们想知道数据是否有偏以及是否属于‘尖峰厚尾’的特征，为了一次性统计数值型变量的偏度和峰度，如下：
# 数据的形状特征
# 挑出所有数值型变量
num_variables = sec_cars.columns[sec_cars.dtypes !='object'][1:]
# 自定义函数，计算偏度和峰度
def skew_kurt(x):
skewness = x.skew()
kurtsis = x.kurt()
# 返回偏度值和峰度值
return pd.Series([skewness,kurtsis], index = ['Skew','Kurt'])
# 运用apply方法
sec_cars[num_variables].apply(func = skew_kurt, axis = 0)
上面我们通过结果可以看到每个数值的偏度和峰值，这三个变量都属于右偏，因为偏度都大于0，三个变量都是尖峰，因为峰度都大于0
字符与日期类型处理
# 数据读入
df = pd.read_excel(r'D:\study\pylearn\pandasdata\data_test03.xlsx')
# 各变量数据类型
print(df.dtypes)
# 将birthday变量转换为日期型
df.birthday = pd.to_datetime(df.birthday, format = '%Y/%m/%d')
# 将手机号转换为字符串
df.tel = df.tel.astype('str')
# 新增年龄和工龄两列
df['age'] = pd.datetime.today().year - df.birthday.dt.year
df['workage'] = pd.datetime.today().year - df.start_work.dt.year
# 将手机号中间四位隐藏起来
df.tel = df.tel.apply(func = lambda x : x.replace(x[3:7], '****'))
# 取出邮箱的域名
df['email_domain'] = df.email.apply(func = lambda x : x.split('@')[1])
# 取出用户的专业信息
df['profession'] = df.other.str.findall('专业：(.*?)，')
# 去除birthday、start_work和other变量
df.drop(['birthday','start_work','other'], axis = 1, inplace = True)
df.head()
# 常用日期处理方法
dates = pd.to_datetime(pd.Series(['1989-8-18 13:14:55','1995-2-16']), format = '%Y-%m-%d %H:%M:%S')
print('返回日期值：\n',dates.dt.date)
print('返回季度：\n',dates.dt.quarter)
print('返回几点钟：\n',dates.dt.hour)
print('返回年中的天：\n',dates.dt.dayofyear)
print('返回年中的周：\n',dates.dt.weekofyear)
print('返回星期几的名称：\n',dates.dt.weekday_name)
print('返回月份的天数：\n',dates.dt.days_in_month)
(2) 冗余数据的识别与处理
data04.duplicated() # 整行有无重复
# data04.duplicated(subset='appname') # 指定名有无重复
data04.drop_duplicates() # 删除整行重复的 #默认不修改源数据
# data04.drop_duplicates(subset='appname') # 删除指定行(指定名重复的)
data04.drop_duplicates(inplace=True) # 修改为替换源数据
data04
(3) 异常值的识别与处理：
z得分法(标准差)，分位数法（箱线图），距离法
1.z得分法（要求数据服从正态分布）
以μ±nσ划线，落在外面认为异常
原理：正态分布中，μ±σ 68%；μ±2σ 95%；μ±3σ 99%+：
2.分位数法
如果数据不服从正态时使用，类似Z得分
# 数据读入
sunspots = pd.read_table(r'D:\study\pylearn\pandasdata\sunspots.csv', sep = ',')
# 异常值检测之标准差法
xbar = sunspots.counts.mean()
xstd = sunspots.counts.std()
print('标准差法异常值上限检测：\n',any(sunspots.counts > xbar + 2 * xstd))
print('标准差法异常值下限检测：\n',any(sunspots.counts < xbar - 2 * xstd))
# 异常值检测之箱线图法
Q1 = sunspots.counts.quantile(q = 0.25)
Q3 = sunspots.counts.quantile(q = 0.75)
IQR = Q3 - Q1
print('箱线图法异常值上限检测：\n',any(sunspots.counts > Q3 + 1.5 * IQR))
print('箱线图法异常值下限检测：\n',any(sunspots.counts < Q1 - 1.5 * IQR))
# 导入绘图模块
import matplotlib.pyplot as plt
# 设置绘图风格
plt.style.use('ggplot')
# 绘制直方图
sunspots.counts.plot(kind = 'hist', bins = 30, density = True, stacked=True)
# 绘制核密度图
sunspots.counts.plot(kind = 'kde')
# 图形展现
plt.show()
# 替换法处理异常值
print('异常值替换前的数据统计特征：\n',sunspots.counts.describe())
# 箱线图中的异常值判别上限
UL = Q3 + 1.5 * IQR
print('判别异常值的上限临界值：\n',UL)
# 从数据中找出低于判别上限的最大值
replace_value = sunspots.counts[sunspots.counts < UL].max()
print('用以替换异常值的数据：\n',replace_value)
# 替换超过判别上限异常值
sunspots.counts[sunspots.counts > UL] = replace_value
print('异常值替换后的数据统计特征：\n',sunspots.counts.describe())
3.距离法
手把手教你如何利用K均值聚类实现异常值的识别！
https://mp.weixin.qq.com/s/aWTDJtafY9XHZdHdOUaqXw
KNN除了可以做分类和预测，还知道它可以识别异常值吗？
https://mp.weixin.qq.com/s/728HfX6VFi0tN6MBkFrTsA
(4) 缺失值
打开文件D:\study\pylearn\pandasdata\data_test05.xlsx
import pandas as pd
data05 = pd.read_excel(r'D:\study\pylearn\pandasdata\data_test05.xlsx')
data05.head()
# 查看缺失值
data05.isnull() # 查看每一个元素缺失
data05.isnull().any(axis=0) #axis=0查看每一列有无缺失 axis=1查看每一行有无缺失
data05.isnull().sum(axis=0) # sum计算缺失数量
data05.isnull().sum(axis=0)/data05.shape[0] # 缺失比例
# 缺失值的处理
# 1.直接删除
data05.dropna() # 直接删除 # 不改变源数据
# data05.dropna(inplace=True) # 改变源数据
# 2.填充
# data05.fillna(value = 0) # 直接填充0
# 根据实际填充 # 坑：mode()[0], 一定要写[0]，因为众数可能不止一个，所以返回的是列表，要取第一个
data05.fillna(value = {'gender':data05.gender.mode()[0],'age':data05.age.mean(),'income':data05.income.median()},inplace=True)
data05
5.5、数据的汇总
透视表功能
pd.pivot_table(data,values=None,index=None,columns=None,aggfunc='mean',
fill_value=None,margins=False,dropna=True,margins_name='All')
data # 透视表的源数据集
values # '数值'
index # '行标签'
columns # '列标签'
aggfunc # 数值的统计函数,默认为统计均值,也可以指定numpy的其他统计函数
fill_value # 指定一个标量,用于填充缺失值
margins # bool类型参数,是否显示行或列的总计值,默认False
dropna # bool类型,是都删除整列为缺失的字段,默认为True
margin_name # 指定行或列的总计名称,默认为All
import pandas as pd
data = pd.read_csv(r'D:\study\pylearn\pandasdata\diamonds.csv')
data.head()
# 统计不同color(列字段)的平均价格(列字段)
pd.pivot_table(data,index='color',values='price',aggfunc='mean')
# 统计不同color(列字段)的所有不同clarity(列字段)的总数
pd.pivot_table(data,index='color',columns='clarity',values='price',aggfunc='size')
import numpy as np
# 分组出color和cut,对其进行color数目,carat最小值,price平均值,table最大值的统计
# 通过groupby方法指定分组变量
grouped = data.groupby(by=['color','cut'])
# 对分组变量进行统计汇总
result = grouped.aggregate({'color':np.size,'carat':np.min,
'price':np.mean,'table':np.max})
result
# 可以调整变量名的顺序
result = pd.DataFrame(result,columns=['color','carat','price','table'])
result
# 变量名的修改
result.rename(columns={'color':'counts','carat':'min_weight','price':'avg_price',
'table':'max_table'},inplace=True)
result
5.6、数据的合并与连接
合并
pd.concat(objs,axis=0,join='outer',join_axex=None,
ignore_index=False,keys=None)
objs # 需要合并的对象
axis # 合并的轴 默认是0为行,1为列
join # 数据合并的方式,默认为outer,表示合并所有的数据
join_axex # 合并数据后,指定保留的数据轴
ignore_index # bool类型,表示是否忽略原数据的索引,默认False,如果为True则忽略原索引生成新索引
keys # 为合并后的数据添加新索引,用于区分各个数据部分
# 构造数据集df1和df2
df1 = pd.DataFrame({'name':['张三','李四','王二'], 'age':[21,25,22],
'gender':['男','女','男']})
df2 = pd.DataFrame({'name':['丁一','赵五'], 'age':[23,22], 'gender':['女','女']})
# 数据集的纵向合并
pd.concat([df1,df2] , keys = ['df1','df2']).reset_index().drop(labels='level_1',axis=1).rename(columns={'level_0':'Class'})
# 如果df2数据集中的“姓名变量为Name”
df2 = pd.DataFrame({'Name':['丁一','赵五'], 'age':[23,22], 'gender':['女','女']})
# 数据集的纵向合并 # concat行合并列名需相同(顺序不限)
pd.concat([df1,df2]) # name和Name不同，所以出现下列情况
连接
内连接：用于返回两张表中的共同部分，可理解为集合的交集（缺点：两张表的数据都会有丢失）
左连接：保全左表的所有数据，将右表数据分配到左表中
pd.merge(left, right, how='inner', on=None, left_on=None, right_on=None,
left_index=False, right_index=False, sort=False, suffixes=('_x', '_y'))
left # 主表
right # 辅表
how # 指定连接方式，默认为inner内连，还有其他选项，如左连left、右连right和外连outer
on # 指定连接两张表的共同字段
left_on # 指定主表中需要连接的共同字段
right_on # 指定辅表中需要连接的共同字段
left_index # bool类型，是否将主表中的行索引用作表连接的共同字段，默认为False
right_index # bool类型，是否将辅表中的行索引用作表连接的共同字段，默认为False
sort # bool类型，是否对连接后的数据按照共同字段排序，默认为False
suffixes # 如果数据连接的结果中存在重叠的变量名，则使用各自的前缀进行区分
# 构造数据集
df3 = pd.DataFrame({'id':[1,2,3,4,5],'name':['张三','李四','王二','丁一','赵五'],
'age':[27,24,25,23,25],'gender':['男','男','男','女','女']})
df4 = pd.DataFrame({'Id':[1,2,2,4,4,4,5], 'score':[83,81,87,75,86,74,88],
'kemu':['科目1','科目1','科目2','科目1','科目2','科目3','科目1']})
df5 = pd.DataFrame({'id':[1,3,5],'name':['张三','王二','赵五'],
'income':[13500,18000,15000]})
print(df3)
print(df4)
print(df5)
# 三表的数据连接
# 首先df3和df4连接
merge1 = pd.merge(left = df3, right = df4, how = 'left', left_on='id', right_on='Id')
merge1
# 再将连接结果与df5连接
merge2 = pd.merge(left = merge1, right = df5, how = 'left')
merge2
6、Matplotlib
离散、数值、关系、多图形
6.1、离散
（1）饼图
# 饼图的绘制
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] # 设置微软雅黑字体,防止乱码
# 构造数据
edu = [0.2515,0.3724,0.3336,0.0368,0.0057]
labels = ['中专','大专','本科','硕士','其他']
# 绘制饼图
plt.pie(x = edu, # 绘图数据
labels=labels, # 添加教育水平标签
autopct='%.1f%%' # 设置百分比的格式，这里保留一位小数\
)
# 添加图标题
plt.title('失信用户的教育水平分布')
# 显示图形
plt.show()
# 添加修饰的饼图
explode = [0,0.1,0,0,0] # 生成数据，用于突出显示大专学历人群
colors=['#9999ff','#ff9999','#7777aa','#2442aa','#dd5555'] # 自定义颜色
# 中文乱码和坐标轴负号的处理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 将横、纵坐标轴标准化处理，确保饼图是一个正圆，否则为椭圆
plt.axes(aspect='equal')
# 绘制饼图
plt.pie(x = edu, # 绘图数据
explode=explode, # 突出显示大专人群
labels=labels, # 添加教育水平标签
colors=colors, # 设置饼图的自定义填充色
autopct='%.1f%%', # 设置百分比的格式，这里保留一位小数
pctdistance=0.8, # 设置百分比标签与圆心的距离
labeldistance = 1.1, # 设置教育水平标签与圆心的距离
startangle = 180, # 设置饼图的初始角度
radius = 1.2, # 设置饼图的半径
counterclock = False, # 是否逆时针，这里设置为顺时针方向
wedgeprops = {'linewidth': 1.5, 'edgecolor':'green'},# 设置饼图内外边界的属性值
textprops = {'fontsize':10, 'color':'black'}, # 设置文本标签的属性值
)
# 添加图标题
plt.title('失信用户的受教育水平分布')
# 显示图形
plt.show()
pandas也可以绘图
# 导入第三方模块
import pandas as pd
# 构建序列
data1 = pd.Series({'中专':0.2515,'大专':0.3724,'本科':0.3336,'硕士':0.0368,'其他':0.0057})
# 将序列的名称设置为空字符，否则绘制的饼图左边会出现None这样的字眼
data1.name = ''
# 控制饼图为正圆
plt.axes(aspect = 'equal')
# plot方法对序列进行绘图
data1.plot(kind = 'pie', # 选择图形类型
autopct='%.1f%%', # 饼图中添加数值标签
radius = 1, # 设置饼图的半径
startangle = 180, # 设置饼图的初始角度
counterclock = False, # 将饼图的顺序设置为顺时针方向
title = '失信用户的受教育水平分布', # 为饼图添加标题
wedgeprops = {'linewidth': 1.5, 'edgecolor':'green'}, # 设置饼图内外边界的属性值
textprops = {'fontsize':10, 'color':'black'} # 设置文本标签的属性值
)
# 显示图形
plt.show()
（2）条形图
# 条形图的绘制--垂直条形图
# 读入数据
GDP = pd.read_excel(r'D:\study\pylearn\matplotlibdata\Province GDP 2017.xlsx')
# 设置绘图风格（不妨使用R语言中的ggplot2风格）
plt.style.use('ggplot')
# 绘制条形图
plt.bar(x = range(GDP.shape[0]), # 指定条形图x轴的刻度值
height = GDP.GDP, # 指定条形图y轴的数值
tick_label = GDP.Province, # 指定条形图x轴的刻度标签
color = 'steelblue', # 指定条形图的填充色
)
# 添加y轴的标签
plt.ylabel('GDP(万亿)')
# 添加条形图的标题
plt.title('2017年度6个省份GDP分布')
# 为每个条形图添加数值标签
for x,y in enumerate(GDP.GDP):
plt.text(x,y+0.1,'%s' %round(y,1),ha='center')
# 显示图形
plt.show()
# 条形图的绘制--水平条形图
# 对读入的数据作升序排序
GDP.sort_values(by = 'GDP', inplace = True)
# 绘制条形图
plt.barh(y = range(GDP.shape[0]), # 指定条形图y轴的刻度值
width = GDP.GDP, # 指定条形图x轴的数值
tick_label = GDP.Province, # 指定条形图y轴的刻度标签
color = 'steelblue', # 指定条形图的填充色
)
# 添加x轴的标签
plt.xlabel('GDP(万亿)')
# 添加条形图的标题
plt.title('2017年度6个省份GDP分布')
# 为每个条形图添加数值标签
for y,x in enumerate(GDP.GDP):
plt.text(x+0.1,y,'%s' %round(x,1),va='center')
# 显示图形
plt.show()
# 条形图的绘制--堆叠条形图
# 读入数据
Industry_GDP = pd.read_excel(r'D:\study\pylearn\matplotlibdata\Industry_GDP.xlsx')
# 取出四个不同的季度标签，用作堆叠条形图x轴的刻度标签
Quarters = Industry_GDP.Quarter.unique()
# 取出第一产业的四季度值
Industry1 = Industry_GDP.GPD[Industry_GDP.Industry_Type == '第一产业']
# 重新设置行索引
Industry1.index = range(len(Quarters))
# 取出第二产业的四季度值
Industry2 = Industry_GDP.GPD[Industry_GDP.Industry_Type == '第二产业']
# 重新设置行索引
Industry2.index = range(len(Quarters))
# 取出第三产业的四季度值
Industry3 = Industry_GDP.GPD[Industry_GDP.Industry_Type == '第三产业']
# 绘制堆叠条形图
# 各季度下第一产业的条形图
plt.bar(x = range(len(Quarters)), height=Industry1, color = 'steelblue', label = '第一产业', tick_label = Quarters)
# 各季度下第二产业的条形图
plt.bar(x = range(len(Quarters)), height=Industry2, bottom = Industry1, color = 'green', label = '第二产业')
# 各季度下第三产业的条形图
plt.bar(x = range(len(Quarters)), height=Industry3, bottom = Industry1 + Industry2, color = 'red', label = '第三产业')
# 添加y轴标签
plt.ylabel('生成总值（亿）')
# 添加图形标题
plt.title('2017年各季度三产业总值')
# 显示各产业的图例
plt.legend()
# 显示图形
plt.show()
pandas也可以绘制条形图
# Pandas模块之垂直或水平条形图
# 绘图（此时的数据集在前文已经按各省GDP做过升序处理）
GDP.GDP.plot(kind = 'bar', width = 0.8, rot = 0, color = 'steelblue', title = '2017年度6个省份GDP分布')
# 添加y轴标签
plt.ylabel('GDP（万亿）')
# 添加x轴刻度标签
plt.xticks(range(len(GDP.Province)), #指定刻度标签的位置
GDP.Province # 指出具体的刻度标签值
)
# 为每个条形图添加数值标签
for x,y in enumerate(GDP.GDP):
plt.text(x-0.1,y+0.2,'%s' %round(y,1),va='center')
# 显示图形
plt.show()
seaborn模块之垂直或水平条形图
# seaborn模块之垂直或水平条形图
# 导入第三方模块
import seaborn as sns
sns.barplot(y = 'Province', # 指定条形图x轴的数据
x = 'GDP', # 指定条形图y轴的数据
data = GDP, # 指定需要绘图的数据集
color = 'steelblue', # 指定条形图的填充色
orient = 'horizontal' # 将条形图水平显示
)
# 重新设置x轴和y轴的标签
plt.xlabel('GDP（万亿）')
plt.ylabel('')
# 添加图形的标题
plt.title('2017年度6个省份GDP分布')
# 为每个条形图添加数值标签
for y,x in enumerate(GDP.GDP):
plt.text(x,y,'%s' %round(x,1),va='center')
# 显示图形
plt.show()
# 读入数据
Titanic = pd.read_csv(r'D:\study\pylearn\matplotlibdata\titanic_train.csv')
# 绘制水平交错条形图
sns.barplot(x = 'Pclass', # 指定x轴数据
y = 'Age', # 指定y轴数据
hue = 'Sex', # 指定分组数据
data = Titanic, # 指定绘图数据集
palette = 'RdBu', # 指定男女性别的不同颜色
errcolor = 'blue', # 指定误差棒的颜色
errwidth=2, # 指定误差棒的线宽
saturation = 1, # 指定颜色的透明度，这里设置为无透明度
capsize = 0.05 # 指定误差棒两端线条的宽度
)
# 添加图形标题
plt.title('各船舱等级中男女乘客的年龄差异')
# 显示图形
plt.show()
6.2、连续
（1）直方图
# matplotlib模块绘制直方图
# 检查年龄是否有缺失
any(Titanic.Age.isnull())
# 不妨删除含有缺失年龄的观察
Titanic.dropna(subset=['Age'], inplace=True)
# 绘制直方图
plt.hist(x = Titanic.Age, # 指定绘图数据
bins = 20, # 指定直方图中条块的个数
color = 'steelblue', # 指定直方图的填充色
edgecolor = 'black' # 指定直方图的边框色
)
# 添加x轴和y轴标签
plt.xlabel('年龄')
plt.ylabel('频数')
# 添加标题
plt.title('乘客年龄分布')
# 显示图形
plt.show()
# Pandas模块绘制直方图和核密度图
# 绘制直方图
Titanic.Age.plot(kind = 'hist', bins = 20, color = 'steelblue', edgecolor = 'black', density = True, label = '直方图')
# 绘制核密度图
Titanic.Age.plot(kind = 'kde', color = 'red', label = '核密度图')
# 添加x轴和y轴标签
plt.xlabel('年龄')
plt.ylabel('核密度值')
# 添加标题
plt.title('乘客年龄分布')
# 显示图例
plt.legend()
# 显示图形
plt.show()
# seaborn模块绘制分组的直方图和核密度图
# 取出男性年龄
Age_Male = Titanic.Age[Titanic.Sex == 'male']
# 取出女性年龄
Age_Female = Titanic.Age[Titanic.Sex == 'female']
# 绘制男女乘客年龄的直方图
sns.distplot(Age_Male, bins = 20, kde = False, hist_kws = {'color':'steelblue'}, label = '男性')
# 绘制女性年龄的直方图
sns.distplot(Age_Female, bins = 20, kde = False, hist_kws = {'color':'purple'}, label = '女性')
plt.title('男女乘客的年龄直方图')
# 显示图例
plt.legend()
# 显示图形
plt.show()
# 绘制男女乘客年龄的核密度图
sns.distplot(Age_Male, hist = False, kde_kws = {'color':'red', 'linestyle':'-'},
norm_hist = True, label = '男性')
# 绘制女性年龄的核密度图
sns.distplot(Age_Female, hist = False, kde_kws = {'color':'black', 'linestyle':'--'},
norm_hist = True, label = '女性')
plt.title('男女乘客的年龄核密度图')
# 显示图例
plt.legend()
# 显示图形
plt.show()
（2）箱线图
# 读取数据
Sec_Buildings = pd.read_excel(r'D:\study\pylearn\matplotlibdata\sec_buildings.xlsx')
# 绘制箱线图
plt.boxplot(x = Sec_Buildings.price_unit, # 指定绘图数据
patch_artist=True, # 要求用自定义颜色填充盒形图，默认白色填充
showmeans=True, # 以点的形式显示均值
boxprops = {'color':'black','facecolor':'steelblue'}, # 设置箱体属性，如边框色和填充色
# 设置异常点属性，如点的形状、填充色和点的大小
flierprops = {'marker':'o','markerfacecolor':'red', 'markersize':3},
# 设置均值点的属性，如点的形状、填充色和点的大小
meanprops = {'marker':'D','markerfacecolor':'indianred', 'markersize':4},
# 设置中位数线的属性，如线的类型和颜色
medianprops = {'linestyle':'--','color':'orange'},
labels = [''] # 删除x轴的刻度标签，否则图形显示刻度标签为1
)
# 添加图形标题
plt.title('二手房单价分布的箱线图')
# 显示图形
plt.show()
import numpy as np
# 二手房在各行政区域的平均单价
group_region = Sec_Buildings.groupby('region')
avg_price = group_region.aggregate({'price_unit':np.mean}).sort_values('price_unit', ascending = False)
# 通过循环，将不同行政区域的二手房存储到列表中
region_price = []
for region in avg_price.index:
region_price.append(Sec_Buildings.price_unit[Sec_Buildings.region == region])
# 绘制分组箱线图
plt.boxplot(x = region_price,
patch_artist=True,
labels = avg_price.index, # 添加x轴的刻度标签
showmeans=True,
boxprops = {'color':'black', 'facecolor':'steelblue'},
flierprops = {'marker':'o','markerfacecolor':'red', 'markersize':3},
meanprops = {'marker':'D','markerfacecolor':'indianred', 'markersize':4},
medianprops = {'linestyle':'--','color':'orange'}
)
# 添加y轴标签
plt.ylabel('单价（元）')
# 添加标题
plt.title('不同行政区域的二手房单价对比')
# 显示图形
plt.show()
seaborn模块
# 绘制分组箱线图
sns.boxplot(x = 'region', y = 'price_unit', data = Sec_Buildings,
order = avg_price.index, showmeans=True,color = 'steelblue',
flierprops = {'marker':'o','markerfacecolor':'red', 'markersize':3},
meanprops = {'marker':'D','markerfacecolor':'indianred', 'markersize':4},
medianprops = {'linestyle':'--','color':'orange'}
)
# 更改x轴和y轴标签
plt.xlabel('')
plt.ylabel('单价（元）')
# 添加标题
plt.title('不同行政区域的二手房单价对比')
# 显示图形
plt.show()
（3）小提琴图
# 读取数据
tips = pd.read_csv(r'D:\study\pylearn\matplotlibdata\tips.csv')
# 绘制分组小提琴图
sns.violinplot(x = "total_bill", # 指定x轴的数据
y = "day", # 指定y轴的数据
hue = "sex", # 指定分组变量
data = tips, # 指定绘图的数据集
order = ['Thur','Fri','Sat','Sun'], # 指定x轴刻度标签的顺序
scale = 'count', # 以男女客户数调节小提琴图左右的宽度
split = True, # 将小提琴图从中间割裂开，形成不同的密度曲线；
palette = 'RdBu' # 指定不同性别对应的颜色（因为hue参数为设置为性别变量）
)
# 添加图形标题
plt.title('每天不同性别客户的消费额情况')
# 设置图例
plt.legend(loc = 'upper center', ncol = 2)
# 显示图形
plt.show()
（4）折线图
%matplotlib
# 数据读取
wechat = pd.read_excel(r'D:\study\pylearn\matplotlibdata\wechat.xlsx')
# 绘制单条折线图
plt.plot(wechat.Date, # x轴数据
wechat.Counts, # y轴数据
linestyle = '-', # 折线类型
linewidth = 2, # 折线宽度
color = 'steelblue', # 折线颜色
marker = 'o', # 折线图中添加圆点
markersize = 6, # 点的大小
markeredgecolor='black', # 点的边框色
markerfacecolor='brown') # 点的填充色
# 添加y轴标签
plt.ylabel('人数')
# 添加图形标题
plt.title('每天微信文章阅读人数趋势')
# 显示图形
plt.show()
%matplotlib
# 绘制两条折线图
# 导入模块，用于日期刻度的修改
import matplotlib as mpl
# 绘制阅读人数折线图
plt.plot(wechat.Date, # x轴数据
wechat.Counts, # y轴数据
linestyle = '-', # 折线类型，实心线
color = 'steelblue', # 折线颜色
label = '阅读人数'
)
# 绘制阅读人次折线图
plt.plot(wechat.Date, # x轴数据
wechat.Times, # y轴数据
linestyle = '--', # 折线类型，虚线
color = 'indianred', # 折线颜色
label = '阅读人次'
)
# 获取图的坐标信息
ax = plt.gca()
# 设置日期的显示格式
date_format = mpl.dates.DateFormatter("%m-%d")
ax.xaxis.set_major_formatter(date_format)
# 设置x轴显示多少个日期刻度
# xlocator = mpl.ticker.LinearLocator(10)
# 设置x轴每个刻度的间隔天数
xlocator = mpl.ticker.MultipleLocator(7)
ax.xaxis.set_major_locator(xlocator)
# 为了避免x轴刻度标签的紧凑，将刻度标签旋转45度
plt.xticks(rotation=45)
# 添加y轴标签
plt.ylabel('人数')
# 添加图形标题
plt.title('每天微信文章阅读人数与人次趋势')
# 添加图例
plt.legend()
# 显示图形
plt.show()
pandas模块绘制折线图
# 读取天气数据
weather = pd.read_excel(r'D:\study\pylearn\matplotlibdata\weather.xlsx')
# 统计每月的平均最高气温
data = weather.pivot_table(index = 'month', columns='year', values='high')
# 绘制折线图
data.plot(kind = 'line',
style = ['-','--',':'] # 设置折线图的线条类型
)
# 修改x轴和y轴标签
plt.xlabel('月份')
plt.ylabel('气温')
# 添加图形标题
plt.title('每月平均最高气温波动趋势')
# 显示图形
plt.show()
6.3、关系型数据
（1）散点图
# 读入数据
iris = pd.read_csv(r'D:\study\pylearn\matplotlibdata\iris.csv')
# 绘制散点图
plt.scatter(x = iris.Petal_Width, # 指定散点图的x轴数据
y = iris.Petal_Length, # 指定散点图的y轴数据
color = 'steelblue' # 指定散点图中点的颜色
)
# 添加x轴和y轴标签
plt.xlabel('花瓣宽度')
plt.ylabel('花瓣长度')
# 添加标题
plt.title('鸢尾花的花瓣宽度与长度关系')
# 显示图形
pandas
# Pandas模块绘制散点图
# 绘制散点图
iris.plot(x = 'Petal_Width', y = 'Petal_Length', kind = 'scatter', title = '鸢尾花的花瓣宽度与长度关系')
# 修改x轴和y轴标签
plt.xlabel('花瓣宽度')
plt.ylabel('花瓣长度')
# 显示图形
plt.show()
seaborn
# seaborn模块绘制分组散点图
sns.lmplot(x = 'Petal_Width', # 指定x轴变量
y = 'Petal_Length', # 指定y轴变量
hue = 'Species', # 指定分组变量
data = iris, # 指定绘图数据集
legend_out = False, # 将图例呈现在图框内
truncate=True # 根据实际的数据范围，对拟合线作截断操作
)
# 修改x轴和y轴标签
plt.xlabel('花瓣宽度')
plt.ylabel('花瓣长度')
# 添加标题
plt.title('鸢尾花的花瓣宽度与长度关系')
# 显示图形
plt.show()
（2）气泡图
# 读取数据
Prod_Category = pd.read_excel(r'D:\study\pylearn\matplotlibdata\SuperMarket.xlsx')
# 将利润率标准化到[0,1]之间（因为利润率中有负数），然后加上微小的数值0.001
range_diff = Prod_Category.Profit_Ratio.max()-Prod_Category.Profit_Ratio.min()
Prod_Category['std_ratio'] = (Prod_Category.Profit_Ratio-Prod_Category.Profit_Ratio.min())/range_diff + 0.001
# 绘制办公用品的气泡图
plt.scatter(x = Prod_Category.Sales[Prod_Category.Category == '办公用品'],
y = Prod_Category.Profit[Prod_Category.Category == '办公用品'],
s = Prod_Category.std_ratio[Prod_Category.Category == '办公用品']*1000,
color = 'steelblue', label = '办公用品', alpha = 0.6
)
# 绘制技术产品的气泡图
plt.scatter(x = Prod_Category.Sales[Prod_Category.Category == '技术产品'],
y = Prod_Category.Profit[Prod_Category.Category == '技术产品'],
s = Prod_Category.std_ratio[Prod_Category.Category == '技术产品']*1000,
color = 'indianred' , label = '技术产品', alpha = 0.6
)
# 绘制家具产品的气泡图
plt.scatter(x = Prod_Category.Sales[Prod_Category.Category == '家具产品'],
y = Prod_Category.Profit[Prod_Category.Category == '家具产品'],
s = Prod_Category.std_ratio[Prod_Category.Category == '家具产品']*1000,
color = 'black' , label = '家具产品', alpha = 0.6
)
# 添加x轴和y轴标签
plt.xlabel('销售额')
plt.ylabel('利润')
# 添加标题
plt.title('销售额、利润及利润率的气泡图')
# 添加图例
plt.legend()
# 显示图形
plt.show()
（3）热力图
seaborn
# 读取数据
Sales = pd.read_excel(r"D:\study\pylearn\matplotlibdata\Sales.xlsx")
# 根据交易日期，衍生出年份和月份字段
Sales['year'] = Sales.Date.dt.year
Sales['month'] = Sales.Date.dt.month
# 统计每年各月份的销售总额
Summary = Sales.pivot_table(index = 'month', columns = 'year', values = 'Sales', aggfunc = np.sum)
# 绘制热力图
sns.heatmap(data = Summary, # 指定绘图数据
cmap = 'PuBuGn', # 指定填充色
linewidths = .1, # 设置每个单元格边框的宽度
annot = True, # 显示数值
fmt = '.1e' # 以科学计算法显示数据
)
#添加标题
plt.title('每年各月份销售总额热力图')
# 显示图形
plt.show()
6.4、多个图的合并
# 读取数据
Prod_Trade = pd.read_excel(r"D:\study\pylearn\matplotlibdata\Prod_Trade.xlsx")
# 衍生出交易年份和月份字段
Prod_Trade['year'] = Prod_Trade.Date.dt.year
Prod_Trade['month'] = Prod_Trade.Date.dt.month
# 设置大图框的长和高
plt.figure(figsize = (12,6))
# 设置第一个子图的布局
ax1 = plt.subplot2grid(shape = (2,3), loc = (0,0))
# 统计2012年各订单等级的数量
Class_Counts = Prod_Trade.Order_Class[Prod_Trade.year == 2012].value_counts()
Class_Percent = Class_Counts/Class_Counts.sum()
# 将饼图设置为圆形（否则有点像椭圆）
ax1.set_aspect(aspect = 'equal')
# 绘制订单等级饼图
ax1.pie(x = Class_Percent.values, labels = Class_Percent.index, autopct = '%.1f%%')
# 添加标题
ax1.set_title('各等级订单比例')
# 设置第二个子图的布局
ax2 = plt.subplot2grid(shape = (2,3), loc = (0,1))
# 统计2012年每月销售额
Month_Sales = Prod_Trade[Prod_Trade.year == 2012].groupby(by = 'month').aggregate({'Sales':np.sum})
# 绘制销售额趋势图
Month_Sales.plot(title = '2012年各月销售趋势', ax = ax2, legend = False)
# 删除x轴标签
ax2.set_xlabel('')
# 设置第三个子图的布局
ax3 = plt.subplot2grid(shape = (2,3), loc = (0,2), rowspan = 2)
# 绘制各运输方式的成本箱线图
sns.boxplot(x = 'Transport', y = 'Trans_Cost', data = Prod_Trade, ax = ax3)
# 添加标题
ax3.set_title('各运输方式成本分布')
# 删除x轴标签
ax3.set_xlabel('')
# 修改y轴标签
ax3.set_ylabel('运输成本')
# 设置第四个子图的布局
ax4 = plt.subplot2grid(shape = (2,3), loc = (1,0), colspan = 2)
# 2012年客单价分布直方图
sns.distplot(Prod_Trade.Sales[Prod_Trade.year == 2012], bins = 40, norm_hist = True, ax = ax4, hist_kws = {'color':'steelblue'}, kde_kws=({'linestyle':'--', 'color':'red'}))
# 添加标题
ax4.set_title('2012年客单价分布图')
# 修改x轴标签
ax4.set_xlabel('销售额')
# 调整子图之间的水平间距和高度间距
plt.subplots_adjust(hspace=0.6, wspace=0.3)
# 图形显示
plt.show()
7、线性回归
7.1、一元线性回归
原理
原理就是一元线性方程
实战
# 工作年限与收入之间的散点图
# 导入第三方模块
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# 导入数据集
income = pd.read_csv(r"D:\study\pylearn\线性回归data\Salary_Data.csv")
# 绘制散点图
sns.lmplot(x = 'YearsExperience', y = 'Salary', data = income, ci = None)
# 显示图形
plt.show()
# 简单线性回归模型的参数求解(手撸源算法,一般不用)
# 样本量
n = income.shape[0]
# 计算自变量、因变量、自变量平方、自变量与因变量乘积的和
sum_x = income.YearsExperience.sum()
sum_y = income.Salary.sum()
sum_x2 = income.YearsExperience.pow(2).sum()
xy = income.YearsExperience * income.Salary
sum_xy = xy.sum()
# 根据公式计算回归模型的参数
b = (sum_xy-sum_x*sum_y/n)/(sum_x2-sum_x**2/n)
a = income.Salary.mean()-b*income.YearsExperience.mean()
# 打印出计算结果
print('回归参数a的值：',a)
print('回归参数b的值：',b)
输出:
回归参数a的值： 25792.200198668666
回归参数b的值： 9449.962321455081
# 简单线性回归模型的参数求解(调用statsmodels)
# 导入第三方模块
import statsmodels.api as sm
# 利用收入数据集，构建回归模型
fit = sm.formula.ols('Salary ~ YearsExperience', data = income).fit()
# 返回模型的参数值
fit.params
输出:
Intercept 25792.200199
YearsExperience 9449.962321
dtype: float64
# y=a+bx: Intercept 对应截距a,YearsExperience 对应自变量b
# 所以最后的线性回归模型可以表示为
# y = Intercept + YearsExperience(x)
# y = 25792.20 + 9449.96x
income.columns
# 查看相关系数
income.Salary.corr(income.YearsExperience)
7.2、多元线性回归
原理
原理就是多元线性方程
实战
# 多元线性回归模型的构建和预测
# 导入模块
from sklearn import model_selection
# 导入数据
Profit = pd.read_excel(r"D:\study\pylearn\线性回归data\Predict to Profit.xlsx")
# 将数据集拆分为训练集和测试集
train, test = model_selection.train_test_split(Profit, test_size = 0.2, random_state=1234)
# 根据train数据集建模
model = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + C(State)', data = train).fit()
print('模型的偏回归系数分别为：\n', model.params)
# 删除test数据集中的Profit变量，用剩下的自变量进行预测
test_X = test.drop(labels = 'Profit', axis = 1)
pred = model.predict(exog = test_X)
print('对比预测值和实际值的差异：\n',pd.DataFrame({'Prediction':pred,'Real':test.Profit}))
# 这个结果的模型公式为
# Profit = 58581.52 +0.80RD_Spend - 0.06Administration + 0.01Marketing_Spend + 927.39Florida - 513.47New York
# 发现结果有5个回归参数,但是我们本来只有4个,因为我们把State作为了类别,总共三类(Florida,New York,California),默认砍掉了California(为什么会砍掉呢,因为State本来是作为一个整体的,切分出来的三个类具有高度的一致性,会对结果产生影响)
输出:
模型的偏回归系数分别为：
Intercept 58581.516503
C(State)[T.Florida] 927.394424
C(State)[T.New York] -513.468310
RD_Spend 0.803487
Administration -0.057792
Marketing_Spend 0.013779
dtype: float64
对比预测值和实际值的差异：
Prediction Real
8 150621.345802 152211.77
48 55513.218079 35673.41
14 150369.022458 132602.65
42 74057.015562 71498.49
29 103413.378282 101004.64
44 67844.850378 65200.33
4 173454.059691 166187.94
31 99580.888894 97483.56
13 128147.138396 134307.35
18 130693.433835 124266.90
# 人工自定义构造哑变量(自定义砍掉New York类别)
# 生成由State变量衍生的哑变量
dummies = pd.get_dummies(Profit.State)
# 将哑变量与原始数据集水平合并
Profit_New = pd.concat([Profit,dummies], axis = 1)
# 删除State变量和California变量（因为State变量已被分解为哑变量，New York变量需要作为参照组）
Profit_New.drop(labels = ['State','New York'], axis = 1, inplace = True)
# 拆分数据集Profit_New
train, test = model_selection.train_test_split(Profit_New, test_size = 0.2, random_state=1234)
# 建模
model2 = sm.formula.ols('Profit ~ RD_Spend + Administration + Marketing_Spend + Florida + California', data = train).fit()
print('模型的偏回归系数分别为：\n', model2.params)
# 这样结果就只剩下Florida,California,没有了New York
# 这个结果的模型公式为
# Profit = 58068.05 + 0.80RD_Spend - 0.06Administration + 0.01Marketing_Spend + 1440.86 Florida + 513.47California
输出:
模型的偏回归系数分别为：
Intercept 58068.048193
RD_Spend 0.803487
Administration -0.057792
Marketing_Spend 0.013779
Florida 1440.862734
California 513.468310
dtype: float64
# 查看某一变量Profit与其他变量的相关系数
Profit.drop('State',axis=1).corrwith(Profit['Profit'])
# 查看两两变量间的相关系数
Profit.drop('State',axis=1).corr()
7.3、线性回归模型的假设检验
（1）模型的显著性检验-F检验
-
提出问题的原假设和备择假设
-
在原假设的条件下，构造统计量F
-
根据样本信息，计算统计量的值
-
对比统计量的值和理论F分布的值，当统计量值超过理论值时，拒绝原假设，否则接受原假设
1.提出假设
H0为原假设，该假设认为模型的所有偏回归系数全为0，即认为没有一个自变量可以构成因变量的线性组合，H1为备择假设，正好是假设的对立面，即至少有一个可以构成因变量的线性组合，就F检验而言，研究者往往是更加希望通过数据推翻原假设H0，而接受备择假设H1的结论
2.构造统计量
ESS是残差平方和，误差平方和（因变量的实际值与预测值之间的离差平方和，随模型变化）
RSS是回归离差平方和（因变量的预测值与实际均值的离差平方和，随模型变化）
TSS是总的总的离差平方和（因变量与其均值之间的离差平方和，其值不随模型而改变）
TSS=ESS+RSS（由于TSS值不变，所以ESS与RSS存在负向关系，一增一减）
F的公式中，p和n-p-1分别为RSS与ESS的自由度，模型拟合的越好，ESS越小，RSS越大，则F越大
3.计算统计量
# 导入第三方模块
import numpy as np
# 这里用的train和model2均为上面的
# 1.手撸源算法(一般不用)
# 计算建模数据中，因变量的均值
ybar = train.Profit.mean()
# 统计变量个数和观测个数
p = model2.df_model
n = train.shape[0]
# 计算回归离差平方和
RSS = np.sum((model2.fittedvalues-ybar) ** 2)
# 计算误差平方和
ESS = np.sum(model2.resid ** 2)
# 计算F统计量的值
F = (RSS/p)/(ESS/(n-p-1))
print('F统计量的值：',F)
# 2.直接调用fvalue返回F
# 返回模型中的F值
model2.fvalue
# 手算与调用结果相同
输出:
F统计量的值： 174.637217168447
174.63721715703537
4.对比结果下结论
# 导入模块
from scipy.stats import f
# 计算F分布的理论值
F_Theroy = f.ppf(q=0.95, dfn = p,dfd = n-p-1)
print('F分布的理论值为：',F_Theroy)
输出:
F分布的理论值为： 2.502635007415366
结果显示，F分布的理论值为： 2.50，计算出来的F统计量的值： 174.64，所以应当拒绝原假设，即认为多元线性回归模型是显著的，也就说回归模型的偏回归系数都不全为0
（2）回归系数的显著性检验-t检验
模型通过F检验，只能说明关于因变量的线性组合是合理的，并不能说明每个自变量对因变量都有意义，所以要对回归系数做显著性t检验
1.提出假设
t检验的出发点就是验证每一个自变量是否都对因变量有意义。H0原假设是假定第j变量的偏回归系数为0，即认为该变量不是因变量的有效因素，而备择H1则相反，认为第j变量是影响因变量的有效因素。我们选择H0，驳回H0，从而得到H1
2.构造统计量
3.计算统计量
# 模型的概览信息
model2.summary()
# 模型包含三部分
# 第一部分：模型的信息(模型的判别系数R²,F统计值等)
# 第二部分：偏回归系数的信息(回归系数的估计值Coef,t统计量,回归系数的置信区间等)
# 第三部分：模型误差项的信息
4.对比结果下结论
只有当p值小于0.05时，才能拒绝原假设。
我们发现，只有Intercept和RD_Spend小于0.05，所以说明其余变量均没有通过系数的显著性检验，即在模型中，其余变量不是影响利润的主要因素
7.4、回归模型的诊断
（1）正态性检验
直方图法
# 正态性检验
# 直方图法
# 导入第三方模块
import scipy.stats as stats
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘制直方图
sns.distplot(a = Profit_New.Profit, bins = 10, fit = stats.norm, norm_hist = True,
hist_kws = {'color':'steelblue', 'edgecolor':'black'},
kde_kws = {'color':'black', 'linestyle':'--', 'label':'核密度曲线'},
fit_kws = {'color':'red', 'linestyle':':', 'label':'正态密度曲线'})
# 显示图例
plt.legend()
# 显示图形
plt.show()
# 从图上我们可以看到，核密度曲线与正态密度曲线的趋势比较吻合，故直观上可以认为利润变量服从正态分布
PP图与QQ图
PP图的思想就是比对正态分布的累计概率值和实际分布的累计概率值，而QQ图则比对正态分布的分位数和实际分布的分位数。判断变量是否近似服从正态分布的标准是：如果散点都比较均匀的散落在直线上，就说明变量近似服从正态分布，否则就认为数据不服从正态分布。
# 残差的正态性检验（PP图和QQ图法）
pp_qq_plot = sm.ProbPlot(Profit_New.Profit)
# 绘制PP图
pp_qq_plot.ppplot(line = '45')
plt.title('P-P图')
# 绘制QQ图
pp_qq_plot.qqplot(line = 'q')
plt.title('Q-Q图')
# 显示图形
plt.show()
# 结果显示，不管是PP图还是QQ图，绘制的点均在直线的附近，没有较大的偏移，故认为利润的变量服从正态分布
Shapiro检验和K-S检验
原假设被认为变量服从正态分布
数据低于5000，使用Shapiro，否则使用K-S
# 导入模块
import scipy.stats as stats
stats.shapiro(Profit_New.Profit)
# 下面输出第一个元素为shapiro检验的统计量值，第二个元素是对应的概率值p，由于p大于置信水平0.05，故接受利润变量服从正态分布的原假设
输出:
ShapiroResult(statistic=0.979339599609375, pvalue=0.5378903746604919)
# 生成正态分布和均匀分布随机数
rnorm = np.random.normal(loc = 5, scale=2, size = 10000)
runif = np.random.uniform(low = 1, high = 100, size = 10000)
# 正态性检验
KS_Test1 = stats.kstest(rvs = rnorm, args = (rnorm.mean(), rnorm.std()), cdf = 'norm')
KS_Test2 = stats.kstest(rvs = runif, args = (runif.mean(), runif.std()), cdf = 'norm')
print(KS_Test1)
print(KS_Test2)
# 下面输出说明，正态分布随机数的检验p值大于置信水平0.05，接受原假设，均匀分布随机数的检验p远小于0.05，拒绝原假设。
# 如果因变量的值不满足正态分布，则需要对因变量做某些数学转换，log,开方，开方倒数，倒数，平方，平方倒数等
输出:
KstestResult(statistic=0.006962358648418321, pvalue=0.7147379693091706)
KstestResult(statistic=0.05727391108862778, pvalue=5.914621987044006e-29)
（2）多重共线性检验
# 导入statsmodels模块中的函数
from statsmodels.stats.outliers_influence import variance_inflation_factor
# 自变量X(包含RD_Spend、Marketing_Spend和常数列1)
X = sm.add_constant(Profit_New.loc[:,['RD_Spend','Marketing_Spend']])
# 构造空的数据框，用于存储VIF值
vif = pd.DataFrame()
vif["features"] = X.columns
vif["VIF Factor"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
# 返回VIF值
vif
# 结果显示，两个自变量对应的方差膨胀因子均低于10，说明构建模型的数据并不存在多重共线性。
# 如果发现变量之间存在多重共线性的话，可以考虑删除变量或者重选模型（岭回归，LASSO）
（3）线性相关性检验
确保自变量与因变量之间存在线性关系
# 计算数据集Profit_New中每个自变量与因变量利润之间的相关系数
Profit_New.drop('Profit', axis = 1).corrwith(Profit_New.Profit)
# 通过输出结果可知，只有RD_Spend和Administration有较高的相关系数（0.98和0.74），其他都很低
输出:
RD_Spend 0.978437
Administration 0.205841
Marketing_Spend 0.739307
California -0.083258
Florida 0.088008
dtype: float64
# 散点图矩阵
# 导入模块
import matplotlib.pyplot as plt
import seaborn
# 绘制散点图矩阵（由于California与Florida都是哑变量，故不加入散点图中）
seaborn.pairplot(Profit_New.loc[:,['RD_Spend','Administration','Marketing_Spend','Profit']])
# 显示图形
plt.show()
# 通过图形可知，研发成本与利润之间的散点图几乎为一条上升的直线（左下角）
综合考虑相关系数，散点图矩阵，t检验的结果，最终确定只保留model2中的RD_Spend和Marketing_Spend两个自变量，下面重现对模型进行修订
# 模型修正
model3 = sm.formula.ols('Profit ~ RD_Spend + Marketing_Spend', data = train).fit()
# 模型回归系数的估计值
model3.params
# 根据输出结果可知，多元线性回归的预测方程为：
# Profit = 51902.11 + 0.79RD_Spend + 0.02Marketing_Spend
输出:
Intercept 51902.112471
RD_Spend 0.785116
Marketing_Spend 0.019402
dtype: float64
（4）异常值检验
帽子矩阵，DFFITS准则，Cook距离
# 异常值检验
outliers = model3.get_influence()
# 高杠杆值点（帽子矩阵）
leverage = outliers.hat_matrix_diag
# dffits值
dffits = outliers.dffits[0]
# 学生化残差
resid_stu = outliers.resid_studentized_external
# cook距离
cook = outliers.cooks_distance[0]
# 合并各种异常值检验的统计量值
contat1 = pd.concat([pd.Series(leverage, name = 'leverage'),pd.Series(dffits, name = 'dffits'),
pd.Series(resid_stu,name = 'resid_stu'),pd.Series(cook, name = 'cook')],axis = 1)
# 重设train数据的行索引
train.index = range(train.shape[0])
# 将上面的统计量与train数据集合并
profit_outliers = pd.concat([train,contat1], axis = 1)
profit_outliers.head()
# 计算异常值数量的比例
outliers_ratio = sum(np.where((np.abs(profit_outliers.resid_stu)>2),1,0))/profit_outliers.shape[0]
outliers_ratio
输出:
0.02564102564102564
# 挑选出非异常的观测点
none_outliers = profit_outliers.loc[np.abs(profit_outliers.resid_stu)<=2,]
# 应用无异常值的数据集重新建模
model4 = sm.formula.ols('Profit ~ RD_Spend + Marketing_Spend', data = none_outliers).fit()
model4.params
# 输出结果显示，经过异常点的排除，重新构模的偏回归系数发生了变化，改为
# Profit = 51827.42 + 0.80RD_Spend + 0.02Marketing_Spend
输出：
Intercept 51827.416821
RD_Spend 0.797038
Marketing_Spend 0.017740
dtype: float64
（5）独立性检验
残差的独立性检验，是对因变量y的独立性检验，因为在线性回归模型的等式左右两边都只有y和残差项属于随机变量，如果再加上正态分布，就构成了残差项独立同分布于正态分布的假设。
# Durbin-Watson统计量
# 模型概览
model4.summary()
（6）方差齐性检验
方差齐性是要求模型残差项的方差不随自变量的变动而呈现某种趋势，否则，残差的趋势就可以被自变量刻画。
图形法
# 方差齐性检验
# 设置第一张子图的位置
ax1 = plt.subplot2grid(shape = (2,1), loc = (0,0))
# 绘制散点图
ax1.scatter(none_outliers.RD_Spend, (model4.resid-model4.resid.mean())/model4.resid.std())
# 添加水平参考线
ax1.hlines(y = 0 ,xmin = none_outliers.RD_Spend.min(),xmax = none_outliers.RD_Spend.max(), color = 'red', linestyles = '--')
# 添加x轴和y轴标签
ax1.set_xlabel('RD_Spend')
ax1.set_ylabel('Std_Residual')
# 设置第二张子图的位置
ax2 = plt.subplot2grid(shape = (2,1), loc = (1,0))
# 绘制散点图
ax2.scatter(none_outliers.Marketing_Spend, (model4.resid-model4.resid.mean())/model4.resid.std())
# 添加水平参考线
ax2.hlines(y = 0 ,xmin = none_outliers.Marketing_Spend.min(),xmax = none_outliers.Marketing_Spend.max(), color = 'red', linestyles = '--')
# 添加x轴和y轴标签
ax2.set_xlabel('Marketing_Spend')
ax2.set_ylabel('Std_Residual')
# 调整子图之间的水平间距和高度间距
plt.subplots_adjust(hspace=0.6, wspace=0.3)
# 显示图形
plt.show()
# 如图所示，标准化残差并没有随自变量而呈现喇叭状，所有的散点几乎均匀的分布在参考线y=0附近，所以，可以说明模型的残差项满足方差齐性的前提假设
BP检验
# BP检验
sm.stats.diagnostic.het_breuschpagan(model4.resid, exog_het = model4.model.exog)
# 输出包含4个值，第一个为LM统计量，第二个为统计量对应的概率p，该值大于0.05，说明接受残差方差为常数的原假设，第三个为F统计量，用于检验残差平方项与自变量之间是否独立，如果独立则表明残差方差齐性，第四个值为F统计量的概率值p,同样大于0.05，则进一步说明残差项满足方差齐性的假设
输出：
(1.4675103668308807, 0.48010272699005274, 0.702975123716269, 0.501965974096277)
# 模型预测
# model4对测试集的预测
pred4 = model4.predict(exog = test.loc[:,['RD_Spend','Marketing_Spend']])
# 绘制预测值与实际值的散点图
plt.scatter(x = test.Profit, y = pred4)
# 添加斜率为1，截距项为0的参考线
plt.plot([test.Profit.min(),test.Profit.max()],[test.Profit.min(),test.Profit.max()],
color = 'red', linestyle = '--')
# 添加轴标签
plt.xlabel('实际值')
plt.ylabel('预测值')
# 显示图形
plt.show()
# 预测值与实际值之间的距离差异，如果两者接近，那么得到的散点图一定是在对角线附近波动。如图说明模型的预测结果不错
7.5、最终对产品的利润预测
最终结果为：
Profit = 51827.42 + 0.80RD_Spend + 0.02Marketing_Spend
8、岭回归与LASSO回归
岭回归与LASSO回归：为解决线性回归的弊端
8.1、岭回归
（1）线性回归的缺点
根据线性回归模型的参数估计公式β=(X′X)−1X′y可知，得到β的前提是矩阵X′X可逆，但在实际应用中，可能会出现自变量个数多于样本量或者自变量间存在多重共线性的情况，即X^′X的行列式为0。此时将无法根据公式计算回归系数的估计值β。
（2）岭回归原理
为解决多元线性回归模型中可能存在的不可逆问题，统计学家提出了岭回归模型。该模型解决问题的思路就是在线性回归模型的目标函数之上添加l2正则项（也称为惩罚项）。
（3）岭回归实战-糖尿病
可视化方法确定Lambda
下面利用糖尿病数据集绘制不同Lambda值对应回归系数的折线图
# 导入第三方模块
import pandas as pd
import numpy as np
from sklearn import model_selection
from sklearn.linear_model import Ridge,RidgeCV
import matplotlib.pyplot as plt
# 读取糖尿病数据集
diabetes = pd.read_excel(r"D:\study\pylearn\岭回归和LASSOdata\diabetes.xlsx")
# 构造自变量（剔除患者性别、年龄和因变量）
predictors = diabetes.columns[2:-1]
# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = model_selection.train_test_split(diabetes[predictors], diabetes['Y'],
test_size = 0.2, random_state = 1234 )
# 构造不同的Lambda值
Lambdas = np.logspace(-5, 2, 200)
# 构造空列表，用于存储模型的偏回归系数
ridge_cofficients = []
# 循环迭代不同的Lambda值
for Lambda in Lambdas:
ridge = Ridge(alpha = Lambda, normalize=True)
ridge.fit(X_train, y_train)
ridge_cofficients.append(ridge.coef_)
# 绘制Lambda与回归系数的关系
# 中文乱码和坐标轴负号的处理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
plt.plot(Lambdas, ridge_cofficients)
# 对x轴作对数变换
plt.xscale('log')
# 设置折线图x轴和y轴标签
plt.xlabel('Lambda')
plt.ylabel('Cofficients')
# 图形显示
plt.show()
# 如图得到，发现当Lambda在0.01附近时，绝大多数变量的回归系数趋于稳定，顾认为Lambda值可以选择在0.01附近
交叉验证法确定Lambda
# 岭回归模型的交叉验证
# 设置交叉验证的参数，对于每一个Lambda值，都执行10重交叉验证
ridge_cv = RidgeCV(alphas = Lambdas, normalize=True, scoring='neg_mean_squared_error', cv = 10)
# 模型拟合
ridge_cv.fit(X_train, y_train)
# 返回最佳的lambda值
ridge_best_Lambda = ridge_cv.alpha_
ridge_best_Lambda
# 根据输出结果可知，利用10重交叉验证方法得到的最佳Lambda值为0.0146，与可视化得到的在0.01附近的结果高度一致
# 该值的评判标准，对每一个Lambda计算MSE，然后选最小的，最为Lambda
输出：
0.014649713983072863
# 导入第三方包中的函数
from sklearn.metrics import mean_squared_error
# 基于最佳的Lambda值建模,这里的ridge_best_Lambda就是0.146
ridge = Ridge(alpha = ridge_best_Lambda, normalize=True)
ridge.fit(X_train, y_train)
# 返回岭回归系数
pd.Series(index = ['Intercept'] + X_train.columns.tolist(),data = [ridge.intercept_] + ridge.coef_.tolist())
#根据输出得到糖尿病的预测方程：
# Y = -322.00 + 6.21BMI + 0.93BP - 0.48S1 + 0.20S2 + 0.02S3 + 4.18S4 + 51.42S5 + 0.38S6
输出：
Intercept -321.996227
BMI 6.206337
BP 0.927093
S1 -0.479414
S2 0.202203
S3 0.016912
S4 4.183247
S5 51.424829
S6 0.384342
dtype: float64
# 预测
ridge_predict = ridge.predict(X_test)
# 预测效果验证
RMSE = np.sqrt(mean_squared_error(y_test,ridge_predict))
RMSE
# 结果RMSE = 53.12，RMSE越小，结果越好
输出：
53.11911788753519
8.2、LASSO
（1）原理
岭回归模型解决线性回归模型中矩阵X^′X不可逆的办法是添加l2正则的惩罚项，但缺陷在于始终保留建模时的所有变量，无法降低模型的复杂度。对于此，Lasso回归采用了l1正则的惩罚项。
（2）LASSO实战-糖尿病
# 导入第三方模块中的函数
from sklearn.linear_model import Lasso,LassoCV
# 构造空列表，用于存储模型的偏回归系数
lasso_cofficients = []
for Lambda in Lambdas:
lasso = Lasso(alpha = Lambda, normalize=True, max_iter=10000)
lasso.fit(X_train, y_train)
lasso_cofficients.append(lasso.coef_)
# 绘制Lambda与回归系数的关系
plt.plot(Lambdas, lasso_cofficients)
# 对x轴作对数变换
plt.xscale('log')
# 设置折线图x轴和y轴标签
plt.xlabel('Lambda')
plt.ylabel('Cofficients')
# 显示图形
plt.show()
# 由图可知，当Lambda落在0.05附近时，绝大多数变量的回归系数趋于稳定，所以在这确定合理的Lambda范围，然后通过交叉验证法确定准确的Lambda值
# LASSO回归模型的交叉验证
lasso_cv = LassoCV(alphas = Lambdas, normalize=True, cv = 10, max_iter=10000)
lasso_cv.fit(X_train, y_train)
# 输出最佳的lambda值
lasso_best_alpha = lasso_cv.alpha_
lasso_best_alpha
# 由输出结果知，合理的Lambda为0.0629，然后基于这个Lambda重新构建LASSO模型
输出：
0.06294988990221888
# 基于最佳的lambda值建模
lasso = Lasso(alpha = lasso_best_alpha, normalize=True, max_iter=10000)
lasso.fit(X_train, y_train)
# 返回LASSO回归的系数
pd.Series(index = ['Intercept'] + X_train.columns.tolist(),data = [lasso.intercept_] + lasso.coef_.tolist())
# 输出结果可知，LASSO回归模型方程为：
# Y = -278.56 + 6.19BMI + 0.86BP - 0.13S1 + 0.49S3 + 44.49S5 + 0.32S6
输出：
Intercept -278.560358
BMI 6.188602
BP 0.860826
S1 -0.127627
S2 -0.000000
S3 -0.488408
S4 0.000000
S5 44.487738
S6 0.324076
dtype: float64
# 预测
lasso_predict = lasso.predict(X_test)
# 预测效果验证
RMSE = np.sqrt(mean_squared_error(y_test,lasso_predict))
RMSE
# 输出结果为RMSE=53.06
# 相比岭回归模型的RMSE下降了0.8，在多数情况下，LASSO比岭回归更可靠
输出：
53.061437258225745
# 导入第三方模块
import statsmodels.api as sms
# 为自变量X添加常数列1，用于拟合截距项
X_train2 = sms.add_constant(X_train)
X_test2 = sms.add_constant(X_test)
# 构建多元线性回归模型
linear = sms.OLS(y_train, X_train2).fit()
# 返回线性回归模型的系数
linear.params
# 多元线性回归模型的方程：
# Y = -406.70 + 6.22BMI + 0.94BP - 1.26S1 + 0.90S2 + 0.96S3 + 6.69S4 + 71.61S5 + 0.38S6
输出：
const -406.699716
BMI 6.217649
BP 0.948245
S1 -1.264772
S2 0.901368
S3 0.962373
S4 6.694215
S5 71.614661
S6 0.376004
dtype: float64
# 模型的预测
linear_predict = linear.predict(X_test2)
# 预测效果验证
RMSE = np.sqrt(mean_squared_error(y_test,linear_predict))
RMSE
# RESE=53.43，可以看出，线性回归在三个模型中时最差的，RMSE最大。(当然，这里没有更好的构建多元线性回归模型)
输出：
53.426239397229885
9、Logistics
广义线性回归，解决0-1的分类问题
9.1、分类模型评估方法
混淆矩阵
A：表示正确预测负例的样本个数，用TN表示。
B：表示预测为负例但实际为正例的个数，用FN表示。
C：表示预测为正例但实际为负例的个数，用FP表示。
D：表示正确预测正例的样本个数，用TP表示。
A+B：表示预测负例的样本个数，用PN表示。
C+D：表示预测正例的样本个数，用PP表示。
A+C：表示实际负例的样本个数，用AN表示。
B+D：表示实际正例的样本个数，用AP表示。
准确率：表示正确预测的正负例样本数与所有样本数量的比值，即(A+D)/(A+B+C+D)。
正例覆盖率：表示正确预测的正例数在实际正例数中的比例，即D/(B+D) 。
负例覆盖率：表示正确预测的负例数在实际负例数中的比例，即A/(A+C)。
正例命中率：表示正确预测的正例数在预测正例数中的比例，即D/(C+D)，
ROC曲线
KS曲线
# python 没有自带KS曲线的实现，我们自己实现 一下
# 自定义绘制ks曲线的函数
def plot_ks(y_test, y_score, positive_flag):
# 对y_test,y_score重新设置索引
y_test.index = np.arange(len(y_test))
#y_score.index = np.arange(len(y_score))
# 构建目标数据集
target_data = pd.DataFrame({'y_test':y_test, 'y_score':y_score})
# 按y_score降序排列
target_data.sort_values(by = 'y_score', ascending = False, inplace = True)
# 自定义分位点
cuts = np.arange(0.1,1,0.1)
# 计算各分位点对应的Score值
index = len(target_data.y_score)*cuts
scores = target_data.y_score.iloc[index.astype('int')]
# 根据不同的Score值，计算Sensitivity和Specificity
Sensitivity = []
Specificity = []
for score in scores:
# 正例覆盖样本数量与实际正例样本量
positive_recall = target_data.loc[(target_data.y_test == positive_flag) & (target_data.y_score>score),:].shape[0]
positive = sum(target_data.y_test == positive_flag)
# 负例覆盖样本数量与实际负例样本量
negative_recall = target_data.loc[(target_data.y_test != positive_flag) & (target_data.y_score<=score),:].shape[0]
negative = sum(target_data.y_test != positive_flag)
Sensitivity.append(positive_recall/positive)
Specificity.append(negative_recall/negative)
# 构建绘图数据
plot_data = pd.DataFrame({'cuts':cuts,'y1':1-np.array(Specificity),'y2':np.array(Sensitivity),
'ks':np.array(Sensitivity)-(1-np.array(Specificity))})
# 寻找Sensitivity和1-Specificity之差的最大值索引
max_ks_index = np.argmax(plot_data.ks)
plt.plot([0]+cuts.tolist()+[1], [0]+plot_data.y1.tolist()+[1], label = '1-Specificity')
plt.plot([0]+cuts.tolist()+[1], [0]+plot_data.y2.tolist()+[1], label = 'Sensitivity')
# 添加参考线
plt.vlines(plot_data.cuts[max_ks_index], ymin = plot_data.y1[max_ks_index],
ymax = plot_data.y2[max_ks_index], linestyles = '--')
# 添加文本信息
plt.text(x = plot_data.cuts[max_ks_index]+0.01,
y = plot_data.y1[max_ks_index]+plot_data.ks[max_ks_index]/2,
s = 'KS= %.2f' %plot_data.ks[max_ks_index])
# 显示图例
plt.legend()
# 显示图形
plt.show()
导入数据测试一下：
# 导入第三方包
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# 导入虚拟数据
virtual_data = pd.read_excel(r"D:\study\pylearn\logisticsdata\virtual_data.xlsx")
# 应用自定义函数绘制k-s曲线
plot_ks(y_test = virtual_data.Class, y_score = virtual_data.Score,positive_flag = 'P')
# 如图所示，两条折线分别代表各分位点下的正例覆盖率和1-负例覆盖率，通过两条曲线很难对模型的好坏做评估，一般会选择KS值作为衡量指标。
# KS的计算公式为KS = Sensitivity - (1 - Specificity) = Sensitivity + Specificity + 1
# 对于KS值而言，越大越好，通常KS大于0.4，模型就基本可以接受
9.2、logistics回归实战-预测运动状态
# 导入第三方模块
import pandas as pd
import numpy as np
from sklearn import linear_model
from sklearn import model_selection
# 读取数据
sports = pd.read_csv(r"D:\study\pylearn\logisticsdata\Run or Walk.csv")
# 提取出所有自变量名称
predictors = sports.columns[4:]
# 构建自变量矩阵
X = sports.loc[:,predictors]
# 提取y变量值
y = sports.activity
# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size = 0.25, random_state = 1234)
# 利用训练集建模
sklearn_logistic = linear_model.LogisticRegression()
sklearn_logistic.fit(X_train, y_train)
# 返回模型的各个参数
print(sklearn_logistic.intercept_, sklearn_logistic.coef_)
# 根据输出得到Logistics回归模型：
# Xβ = 4.36 + 0.49acceleration(下标x) + 6.88acceleration(下标y) -2.45acceleration(下标z) - 0.01gyro_x -0.16gyro_y + 0.13gyro_z
输出：
[4.36637441] [[ 0.48695898 6.87517973 -2.44872468 -0.01385936 -0.16085022 0.13389695]]
# 模型预测
sklearn_predict = sklearn_logistic.predict(X_test)
# 预测结果统计
pd.Series(sklearn_predict).value_counts()
# 输出结果，显示步行状态有12119个，跑步状态有10028个，但看这两个数据无法判断准确性，所以需要下面的评估
输出：
0 12119
1 10028
dtype: int64
# 导入第三方模块
from sklearn import metrics
# 混淆矩阵
cm = metrics.confusion_matrix(y_test, sklearn_predict, labels = [0,1])
cm
# 输出可知，混淆矩阵，行表示的是实际的运动状态，列表示预测的运动状态，进而计算接下来的
输出：
array([[9969, 1122],
[2150, 8906]], dtype=int64)
Accuracy = metrics._scorer.accuracy_score(y_test, sklearn_predict)
Sensitivity = metrics._scorer.recall_score(y_test, sklearn_predict)
Specificity = metrics._scorer.recall_score(y_test, sklearn_predict, pos_label=0)
print('模型准确率为%.2f%%:' %(Accuracy*100))
print('正例覆盖率为%.2f%%' %(Sensitivity*100))
print('负例覆盖率为%.2f%%' %(Specificity*100))
# 计算得到了模型准确率为85.23%:，正例80.55，负例89.88，所以预测准确率非常好，当然，饿哦们也可以用混淆矩阵的可视化处理
输出：
模型准确率为85.23%:
正例覆盖率为80.55%
负例覆盖率为89.88%
# 混淆矩阵的可视化
# 导入第三方模块
import seaborn as sns
import matplotlib.pyplot as plt
# 绘制热力图
sns.heatmap(cm, annot = True, fmt = '.2e',cmap = 'GnBu')
# 图形显示
plt.show()
# 非常好看到主对角线的区域颜色明显，代表样本数多，说明准确预测的正例和负例样本数都很大
# 使用可视化对模型进行评估
# 1.首先使用ROC曲线
# y得分为模型预测正例的概率
y_score = sklearn_logistic.predict_proba(X_test)[:,1]
# 计算不同阈值下，fpr和tpr的组合值，其中fpr表示1-Specificity，tpr表示Sensitivity
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# 由图知，曲线下的面积高达0.93，远超评估标准0.8
# 2.调用自定义函数，绘制K-S曲线
plot_ks(y_test = y_test, y_score = y_score, positive_flag = 1)
# 由图可知，在40%分位处，KS = 0.71，通常KS=0.4就可以，达到0.71远超0.4
接下来使用另一种库statsmodels实现logistics
# -----------------------第一步 建模 ----------------------- #
# 导入第三方模块
import statsmodels.api as sm
# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size = 0.25, random_state = 1234)
# 为训练集和测试集的X矩阵添加常数列1
X_train2 = sm.add_constant(X_train)
X_test2 = sm.add_constant(X_test)
# 拟合Logistic模型
sm_logistic = sm.Logit(y_train, X_train2).fit()
# 返回模型的参数
sm_logistic.params
输出：
const 4.388537
acceleration_x 0.489617
acceleration_y 6.906590
acceleration_z -2.459638
gyro_x -0.014715
gyro_y -0.161164
gyro_z 0.134655
dtype: float64
# -----------------------第二步 预测构建混淆矩阵 ----------------------- #
# 模型在测试集上的预测
sm_y_probability = sm_logistic.predict(X_test2)
# 根据概率值，将观测进行分类，以0.5作为阈值
sm_pred_y = np.where(sm_y_probability >= 0.5, 1, 0)
# 混淆矩阵
cm = metrics.confusion_matrix(y_test, sm_pred_y, labels = [0,1])
cm
输出：
array([[9967, 1124],
[2149, 8907]], dtype=int64)
# -----------------------第三步 绘制ROC曲线 ----------------------- #
# 计算真正率和假正率
fpr,tpr,threshold = metrics.roc_curve(y_test, sm_y_probability)
# 计算auc的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# -----------------------第四步 绘制K-S曲线 ----------------------- #
# 调用自定义函数，绘制K-S曲线
sm_y_probability.index = np.arange(len(sm_y_probability))
plot_ks(y_test = y_test, y_score = sm_y_probability, positive_flag = 1)
10、决策树与随机森林
决策树：
决策树属于经典的十大数据挖掘算法之一，是一种类似于流程图的树结构，其规则就是IF...THEN...的思想，可以用于数值型因变量的预测和离散型因变量的分类。
该算法简单直观、通俗易懂，不需要研究者掌握任何领域知识或复杂的数学推理，而且算法的结果输出具有很强的解释性。
随机森林：
利用Bootstrap抽样法，从原始数据集中生成k个数据集，并且每个数据集都含有N个观测和P个自变量。
针对每一个数据集，构造一棵CART决策树，在构建子树的过程中，并没有将所有自变量用作节点字段的选择，而是随机选择p个字段。
让每一棵决策树尽可能地充分生长，使得树中的每个节点尽可能“纯净”，即随机森林中的每一棵子树都不需要剪枝。
针对k棵CART树的随机森林，对分类问题利用投票法，将最高得票的类别用于最终的判断结果；对回归问题利用均值法，将其用作预测样本的最终结果。
10.1、分类-乘客幸存预测
（1）预处理数据
# 导入第三方模块
import pandas as pd
# 读入数据
Titanic = pd.read_csv(r"D:\study\pylearn\决策树随机森林data\Titanic.csv")
Titanic.head()
# 删除无意义的变量，并检查剩余自字是否含有缺失值
Titanic.drop(['PassengerId','Name','Ticket','Cabin'], axis = 1, inplace = True)
Titanic.isnull().sum(axis = 0)
# 通过输出发现，Age有177个缺失值，Embarked有2个缺失项
输出：
Survived 0
Pclass 0
Sex 0
Age 177
SibSp 0
Parch 0
Fare 0
Embarked 2
dtype: int64
# 对Age使用均值填充，Embarked使用众数填充，由于Age缺失非常多，故不直接使用均值，采用分组均值填充
# 对Sex分组，用各组乘客的平均年龄填充各组中的缺失年龄
fillna_Titanic = []
for i in Titanic.Sex.unique():
update = Titanic.loc[Titanic.Sex == i,].fillna(value = {'Age': Titanic.Age[Titanic.Sex == i].mean()}, inplace = False)
fillna_Titanic.append(update)
Titanic = pd.concat(fillna_Titanic)
# 使用Embarked变量的众数填充缺失值
Titanic.fillna(value = {'Embarked':Titanic.Embarked.mode()[0]}, inplace=True)
Titanic.head()
# 对Sex和Embarked进行哑变量处理
# 将数值型的Pclass转换为类别型，否则无法对其哑变量处理
Titanic.Pclass = Titanic.Pclass.astype('category')
# 哑变量处理
dummy = pd.get_dummies(Titanic[['Sex','Embarked','Pclass']])
# 水平合并Titanic数据集和哑变量的数据集
Titanic = pd.concat([Titanic,dummy], axis = 1)
# 删除原始的Sex、Embarked和Pclass变量
Titanic.drop(['Sex','Embarked','Pclass'], inplace=True, axis = 1)
Titanic.head()
（2）构建决策树模型
# 导入第三方包
from sklearn import model_selection
# 取出所有自变量名称
predictors = Titanic.columns[1:]
# 将数据集拆分为训练集和测试集，且测试集的比例为25%
X_train, X_test, y_train, y_test = model_selection.train_test_split(Titanic[predictors], Titanic.Survived,
test_size = 0.25, random_state = 1234)
# 为防止决策树出现过拟合，需要对决策树进行预剪枝
# 导入第三方模块
from sklearn.model_selection import GridSearchCV
from sklearn import tree
# 预设各参数的不同选项值
max_depth = [2,3,4,5,6]
min_samples_split = [2,4,6,8]
min_samples_leaf = [2,4,8,10,12]
# 将各参数值以字典形式组织起来
parameters = {'max_depth':max_depth, 'min_samples_split':min_samples_split, 'min_samples_leaf':min_samples_leaf}
# 网格搜索法，测试不同的参数值，从而选择最佳的参数组合
grid_dtcateg = GridSearchCV(estimator = tree.DecisionTreeClassifier(), param_grid = parameters, cv=10)
# 模型拟合
grid_dtcateg.fit(X_train, y_train)
# 返回最佳组合的参数值
grid_dtcateg.best_params_
# 由输出结果可知，经过10重交叉验证的网格搜索，得到各参数的最佳组合值为6,8,6
# 根据经验，数据量较小，可以设置在10以内，反之20左右
输出：
{'max_depth': 6, 'min_samples_leaf': 8, 'min_samples_split': 6}
# 导入第三方模块
from sklearn import metrics
# 构建分类决策树（6,8,6）
CART_Class = tree.DecisionTreeClassifier(max_depth=6, min_samples_leaf = 8, min_samples_split=6)
# 模型拟合
decision_tree = CART_Class.fit(X_train, y_train)
# 模型在测试集上的预测
pred = CART_Class.predict(X_test)
# 模型的准确率
print('模型在测试集的预测准确率：\n',metrics.accuracy_score(y_test, pred))
# 根据输出，得到准确率为84.8%，还是很不错的
输出：
模型在测试集的预测准确率：
0.8475336322869955
# 为进一步验证预测效果，绘制ROC曲线
# 导入第三方包
import matplotlib.pyplot as plt
y_score = CART_Class.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# 根据ROC曲线的面积AUC为0.88，超过0.8，可以认为模型的拟合效果比较理想
# 决策树实际上是一个if then逻辑，为进一步展示决策树背后的逻辑，将决策树可视化展示
# 需要在电脑中安装Graphviz和配置环境变量（安装时可选）
# https://graphviz.gitlab.io/_pages/Download/Download_windows.html
# 导入第三方模块
from sklearn.tree import export_graphviz
from IPython.display import Image
import pydotplus
# from sklearn.externals.six import StringIO
from six import StringIO
# 绘制决策树
dot_data = StringIO()
export_graphviz(
decision_tree,
out_file=dot_data,
feature_names=predictors,
class_names=['Unsurvived','Survived'],
# filled=True,
rounded=True,
special_characters=True
)
# 决策树展现
graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
Image(graph.create_png())
局部展示（最上面两个）：
（3）构建随机森林模型
# 导入第三方包
from sklearn import ensemble
# 构建随机森林
RF_class = ensemble.RandomForestClassifier(n_estimators=200, random_state=1234)
# 随机森林的拟合
RF_class.fit(X_train, y_train)
# 模型在测试集上的预测
RFclass_pred = RF_class.predict(X_test)
# 模型的准确率
print('模型在测试集的预测准确率：\n',metrics.accuracy_score(y_test, RFclass_pred))
# 输出结果显示，利用随机森林进行分类，确实提高了准确率
输出：
模型在测试集的预测准确率：
0.852017937219731
# 绘制ROC曲线
# 计算绘图数据
y_score = RF_class.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
roc_auc = metrics.auc(fpr,tpr)
# 绘图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
plt.plot(fpr, tpr, color='black', lw = 1)
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
plt.show()
# AUC值为0.88，同样比单棵决策树的AUC高（这次结果相同）
# 最后利用理想的随机森林算法挑选出影响乘客是否幸存的重要因素
# 变量的重要性程度值
importance = RF_class.feature_importances_
# 构建含序列用于绘图
Impt_Series = pd.Series(importance, index = X_train.columns)
# 对序列排序绘图
Impt_Series.sort_values(ascending = True).plot(kind='barh')
plt.show()
# 由图知，最重要的三个因素为Age年龄，Fare票价,Sex_female是否为女性
10.2、预测-预测肾功能健康
（1）数据预处理
# 读入数据，数据已做过预处理
NHANES = pd.read_excel(r"D:\study\pylearn\决策树随机森林data\NHANES.xlsx")
print(NHANES.shape)
NHANES.head()
# 数据是关于患者肾小球滤过率
# 数据中的CKD_epi_eGFR变量为因变量，它是连续性变量
（2）构建决策树模型
# 取出自变量名称
predictors = NHANES.columns[:-1]
# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = model_selection.train_test_split(NHANES[predictors], NHANES.CKD_epi_eGFR,
test_size = 0.25, random_state = 1234)
# 预设各参数的不同选项值
max_depth = [18,19,20,21,22]
min_samples_split = [2,4,6,8]
min_samples_leaf = [2,4,8]
parameters = {'max_depth':max_depth, 'min_samples_split':min_samples_split, 'min_samples_leaf':min_samples_leaf}
# 网格搜索法，测试不同的参数值
grid_dtreg = GridSearchCV(estimator = tree.DecisionTreeRegressor(), param_grid = parameters, cv=10)
# 模型拟合
grid_dtreg.fit(X_train, y_train)
# 返回最佳组合的参数值
grid_dtreg.best_params_
# 由输出可知，参数组合为18,2,6
输出：
{'max_depth': 18, 'min_samples_leaf': 2, 'min_samples_split': 6}
# 构建用于回归的决策树
CART_Reg = tree.DecisionTreeRegressor(max_depth = 18, min_samples_leaf = 2, min_samples_split = 6)
# 回归树拟合
CART_Reg.fit(X_train, y_train)
# 模型在测试集上的预测
pred = CART_Reg.predict(X_test)
# 计算衡量模型好坏的MSE值
metrics.mean_squared_error(y_test, pred)
# MSE结果为1.93
输出：
1.9268545203952838
（3）构建决策森林
# 构建用于回归的随机森林
RF = ensemble.RandomForestRegressor(n_estimators=200, random_state=1234)
# 随机森林拟合
RF.fit(X_train, y_train)
# 模型在测试集上的预测
RF_pred = RF.predict(X_test)
# 计算模型的MSE值
metrics.mean_squared_error(y_test, RF_pred)
# MSE=0.89，比单棵树小一半多
输出：
0.895881530212451
# 基于随机森林的重要变量绘图
# 构建变量重要性的序列
importance = pd.Series(RF.feature_importances_, index = X_train.columns)
# 排序并绘图
importance.sort_values().plot(kind='barh')
plt.show()
# 根据图结果来看，影响结果预测准确率的三个主要因素是age_month年龄，cal_create某尿液细胞指标,CKD_stage慢性肾脏病
11、KNN
思想
寻找最近的k个邻居，基于最近样本做预测
步骤
1.确定k值（k过小，过拟合，k过大，欠拟合）
2.根据某种度量样本间相似度的指标（如欧氏距离（x1y1-x2y2的距离））将每一个未知样本的最近k个已知样本搜寻出来，形成一个簇，
3.对搜寻出来的已知样本进行投票，将各簇下类别最多的分类用作未知样本点的预测
k值的选择
权重法：当k比较大时，担心欠拟合，设置k近邻样本的投票权重（近的权重高，远的权重低）：
多重交叉验证法：k取不同的值，然后每种执行m重交叉验证，最后选出平均误差最小的k值。
距离度量
欧氏距离，x1y1-x2y2直线距离，勾股定理
曼哈顿距离，水平垂直距离，不斜着
余弦相似度，夹角查越小，越相似
实战-分类-学生成绩的高低
# 导入第三方包
import pandas as pd
# 导入数据
Knowledge = pd.read_excel(r"D:\study\pylearn\knndata\Knowledge.xlsx")
# 返回前5行数据
Knowledge.head()
# 行代表每一个被观察的学生，前五列分别为学习时长STG，重复次数SCG，相关科目的学习时长STR，相关科目的考试成绩LPR，目标科目的考试成绩LPR，最后一列为知识掌握程度UNS（四种程度Very Low，Low，Middle，High）
# 构造训练集和测试集
# 导入第三方模块
from sklearn import model_selection
# 将数据集拆分为训练集和测试集
predictors = Knowledge.columns[:-1]
X_train, X_test, y_train, y_test = model_selection.train_test_split(Knowledge[predictors], Knowledge.UNS,
test_size = 0.25, random_state = 1234)
# 导入第三方模块
import numpy as np
from sklearn import neighbors
import matplotlib.pyplot as plt
# 设置待测试的不同k值
K = np.arange(1,np.ceil(np.log2(Knowledge.shape[0])))
# 构建空的列表，用于存储平均准确率
accuracy = []
for k in K:
# 使用10重交叉验证的方法，比对每一个k值下KNN模型的预测准确率
cv_result = model_selection.cross_val_score(neighbors.KNeighborsClassifier(n_neighbors = int(k), weights = 'distance'),
X_train, y_train, cv = 10, scoring='accuracy')
accuracy.append(cv_result.mean())
# 从k个平均准确率中挑选出最大值所对应的下标
arg_max = np.array(accuracy).argmax()
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘制不同K值与平均预测准确率之间的折线图
plt.plot(K, accuracy)
# 添加点图
plt.scatter(K, accuracy)
# 添加文字说明
plt.text(K[arg_max], accuracy[arg_max], '最佳k值为%s' %int(K[arg_max]))
# 显示图形
plt.show()
# 重新构建模型，并将最佳的近邻个数设置为5
knn_class = neighbors.KNeighborsClassifier(n_neighbors = 5, weights = 'distance')
# 模型拟合
knn_class.fit(X_train, y_train)
# 模型在测试数据集上的预测
predict = knn_class.predict(X_test)
# 构建混淆矩阵
cm = pd.crosstab(predict,y_test)
cm
# 混淆矩阵，从对角线看，绝大多数样本都被正确分类，我们使用热力图来更好的观察
# 绘制热力图
# 导入第三方模块
import seaborn as sns
# 将混淆矩阵构造成数据框，并加上字段名和行名称，用于行或列的含义说明
cm = pd.DataFrame(cm)
# 绘制热力图
sns.heatmap(cm, annot = True,cmap = 'GnBu')
# 添加x轴和y轴的标签
plt.xlabel(' Real Lable')
plt.ylabel(' Predict Lable')
# 图形显示
plt.show()
# 每一行表示真实的样本类别，每一列表示预测的样本类别，颜色越深，代表对应数值越高，可以看到，对角线颜色都比较深，说明大部分预测分类正确。
from sklearn import metrics
# 模型整体的预测准确率
metrics._scorer.accuracy_score(y_test, predict)
# 预测的准确率为91.09%
# 准确率的计算公式其实就是，混淆矩阵上的数字和与所有数字和的商
输出：
0.9108910891089109
# 上面计算了整体的准确率，无法对每个类的精度进行计算，接下来计算各个类的预测效果
# 分类模型的评估报告
print(metrics.classification_report(y_test, predict))
# 评估报告，前四行代表因变量y中的各个类别值，最后一行为各指标的综合水平
# 第一列precision表示模型的预测精度（预测正确的个数/该类的预测所有个数）
# 第二列recall表示预测覆盖率（预测正确的个数/该类的实际所有个数）
# 第三列f1-score表示precision和recall的加权
# 第四列support表示实际的样本个数
实战-预测-预测发电量
# 读入数据
ccpp = pd.read_excel(r"D:\study\pylearn\knndata\CCPP.xlsx")
# 返回数据集的行数与列数
print(ccpp.shape)
ccpp.head()
# 前四个为自变量AT高炉温度，V炉内压力，AP高炉的相对湿度，RH高炉的排气量，PE为连续的因变量
# 导入第三方包
from sklearn.preprocessing import minmax_scale
# 对所有自变量数据作标准化处理
predictors = ccpp.columns[:-1]
X = minmax_scale(ccpp[predictors])
# 将数据集拆分为训练集和测试集
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, ccpp.PE,
test_size = 0.25, random_state = 1234)
# 设置待测试的不同k值
K = np.arange(1,np.ceil(np.log2(ccpp.shape[0])))
# 构建空的列表，用于存储平均MSE
mse = []
for k in K:
# 使用10重交叉验证的方法，比对每一个k值下KNN模型的计算MSE
cv_result = model_selection.cross_val_score(neighbors.KNeighborsRegressor(n_neighbors = int(k), weights = 'distance'),
X_train, y_train, cv = 10, scoring='neg_mean_squared_error')
mse.append((-1*cv_result).mean())
# 从k个平均MSE中挑选出最小值所对应的下标
arg_min = np.array(mse).argmin()
# 绘制不同K值与平均MSE之间的折线图
plt.plot(K, mse)
# 添加点图
plt.scatter(K, mse)
# 添加文字说明
plt.text(K[arg_min], mse[arg_min] + 0.5, '最佳k值为%s' %int(K[arg_min]))
# 显示图形
plt.show()
# 最佳k值为7
# 重新构建模型，并将最佳的近邻个数设置为7
knn_reg = neighbors.KNeighborsRegressor(n_neighbors = 7, weights = 'distance')
# 模型拟合
knn_reg.fit(X_train, y_train)
# 模型在测试集上的预测
predict = knn_reg.predict(X_test)
# 计算MSE值
metrics.mean_squared_error(y_test, predict)
# MSE=12.81
输出：
12.814094947334913
# 对比真实值和实际值
pd.DataFrame({'Real':y_test,'Predict':predict}, columns=['Real','Predict']).head(10)
# KNN的MSE=12.81，接下来我们使用随机森林进行预测，来对比一下
# 导入第三方模块
from sklearn import tree
# 预设各参数的不同选项值
max_depth = [19,21,23,25,27]
min_samples_split = [2,4,6,8]
min_samples_leaf = [2,4,8,10,12]
parameters = {'max_depth':max_depth, 'min_samples_split':min_samples_split, 'min_samples_leaf':min_samples_leaf}
# 网格搜索法，测试不同的参数值
grid_dtreg = model_selection.GridSearchCV(estimator = tree.DecisionTreeRegressor(), param_grid = parameters, cv=10)
# 模型拟合
grid_dtreg.fit(X_train, y_train)
# 返回最佳组合的参数值
grid_dtreg.best_params_
# 结果集为21,10,6
输出：
{'max_depth': 21, 'min_samples_leaf': 10, 'min_samples_split': 6}
# 构建用于回归的决策树
CART_Reg = tree.DecisionTreeRegressor(max_depth = 21, min_samples_leaf = 10, min_samples_split = 6)
# 回归树拟合
CART_Reg.fit(X_train, y_train)
# 模型在测试集上的预测
pred = CART_Reg.predict(X_test)
# 计算衡量模型好坏的MSE值
metrics.mean_squared_error(y_test, pred)
# 决策树的MSE=16.18，对比KNN的MSE=12.81大了一些，说明在这个模型上，KNN较好一点
输出：
16.175987702804072
12、朴素贝叶斯
思想
通过已知类别的训练数据集，计算样本的先验概率，然后利用贝叶斯概率公式测算未知类别样本属于某个类别的后验概率，最终以最大后验概率所对应的类别作为样本的预测值。
12.1、高斯贝叶斯分类器
# 导入第三方包
import pandas as pd
# 读入数据
skin = pd.read_excel(r"D:\study\pylearn\朴素贝叶斯data\Skin_Segment.xlsx")
# 设置正例和负例，把2设成0，把1设成1
skin.y = skin.y.map({2:0,1:1})
skin.y.value_counts()
# 输出显示，0为非人类面部皮肤有194198个，1为人类面部皮肤有50859个
输出：
0 194198
1 50859
Name: y, dtype: int64
# 导入第三方模块
from sklearn import model_selection
# 样本拆分
X_train,X_test,y_train,y_test = model_selection.train_test_split(skin.iloc[:,:3], skin.y, test_size = 0.25, random_state=1234)
# 导入第三方模块
from sklearn import naive_bayes
# 调用高斯朴素贝叶斯分类器的“类”
gnb = naive_bayes.GaussianNB()
# 模型拟合
gnb.fit(X_train, y_train)
# 模型在测试数据集上的预测
gnb_pred = gnb.predict(X_test)
# 各类别的预测数量
pd.Series(gnb_pred).value_counts()
# 结果显示，预测为负例一共有50630个，预测为正例有10635个
输出：
0 50630
1 10635
dtype: int64
# 为检验预测效果，构建混淆矩阵和绘制ROC曲线
# 构建混淆矩阵
# 导入第三方包
from sklearn import metrics
import matplotlib.pyplot as plt
import seaborn as sns
# 构建混淆矩阵
cm = pd.crosstab(gnb_pred,y_test)
# 绘制混淆矩阵图
sns.heatmap(cm, annot = True, cmap = 'GnBu', fmt = 'd')
# 去除x轴和y轴标签
plt.xlabel('Real')
plt.ylabel('Predict')
# 显示图形
plt.show()
print('模型的准确率为：\n',metrics.accuracy_score(y_test, gnb_pred))
print('模型的评估报告：\n',metrics.classification_report(y_test, gnb_pred))
# 准确率为92.30%，评估报告也可以得到每个类型的准确率
# 绘制ROC曲线
# 计算正例的预测概率，用于生成ROC曲线的数据
y_score = gnb.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# 得到AUC的值为0.94，超过评价模型好坏的阈值0.8，认为非常理想。
12.2、多项式贝叶斯分类器
# 导入第三方包
import pandas as pd
# 读取数据
mushrooms = pd.read_csv(r"D:\study\pylearn\朴素贝叶斯data\mushrooms.csv")
# 数据的前5行
mushrooms.head()
# 所有变量均为字符型的离散值，需要把字符转换为数值
# 将字符型数据作因子化处理，将其转换为整数型数据
columns = mushrooms.columns[1:]
for column in columns:
mushrooms[column] = pd.factorize(mushrooms[column])[0]
mushrooms.head()
from sklearn import model_selection
# 将数据集拆分为训练集合测试集
Predictors = mushrooms.columns[1:]
X_train,X_test,y_train,y_test = model_selection.train_test_split(mushrooms[Predictors], mushrooms['type'],
test_size = 0.25, random_state = 10)
from sklearn import naive_bayes
from sklearn import metrics
import seaborn as sns
import matplotlib.pyplot as plt
# 构建多项式贝叶斯分类器的“类”
mnb = naive_bayes.MultinomialNB()
# 基于训练数据集的拟合
mnb.fit(X_train, y_train)
# 基于测试数据集的预测
mnb_pred = mnb.predict(X_test)
# 构建混淆矩阵
cm = pd.crosstab(mnb_pred,y_test)
# 绘制混淆矩阵图
sns.heatmap(cm, annot = True, cmap = 'GnBu', fmt = 'd')
# 去除x轴和y轴标签
plt.xlabel('Real')
plt.ylabel('Predict')
# 显示图形
plt.show()
# 模型的预测准确率
print('模型的准确率为：\n',metrics.accuracy_score(y_test, mnb_pred))
print('模型的评估报告：\n',metrics.classification_report(y_test, mnb_pred))
# 从结果看，比较理想
from sklearn import metrics
# 计算正例的预测概率，用于生成ROC曲线的数据
y_score = mnb.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test.map({'edible':0,'poisonous':1}), y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# AUC=0.94，超过阈值0.8
12.3、伯努利贝叶斯分类器
import pandas as pd
# 读入评论数据
evaluation = pd.read_excel(r"D:\study\pylearn\朴素贝叶斯data\Contents.xlsx",sheet_name=0)
# 查看数据前10行
evaluation.head(10)
# 运用正则表达式，将评论中的数字和英文去除
evaluation.Content = evaluation.Content.str.replace('[0-9a-zA-Z]','')
evaluation.head()
# 使用正则表达式去除数字英文等内容
# 导入第三方包
import jieba
# 加载自定义词库
jieba.load_userdict(r"D:\study\pylearn\朴素贝叶斯data\all_words.txt")
# 读入停止词
with open(r"D:\study\pylearn\朴素贝叶斯data\mystopwords.txt", encoding='UTF-8') as words:
stop_words = [i.strip() for i in words.readlines()]
# 构造切词的自定义函数，并在切词过程中删除停止词
def cut_word(sentence):
words = [i for i in jieba.lcut(sentence) if i not in stop_words]
# 切完的词用空格隔开
result = ' '.join(words)
return(result)
# 对评论内容进行批量切词
words = evaluation.Content.apply(cut_word)
# 前5行内容的切词效果
words[:5]
# 统计频次
# 导入第三方包
from sklearn.feature_extraction.text import CountVectorizer
# 计算每个词在各评论内容中的次数，并将稀疏度为99%以上的词删除
counts = CountVectorizer(min_df = 0.01)
# 文档词条矩阵
dtm_counts = counts.fit_transform(words).toarray()
# 矩阵的列名称
columns = counts.get_feature_names()
# 将矩阵转换为数据框--即X变量
X = pd.DataFrame(dtm_counts, columns=columns)
# 情感标签变量
y = evaluation.Type
X.head()
# 训练预测，评估效果
from sklearn import model_selection
from sklearn import naive_bayes
from sklearn import metrics
import matplotlib.pyplot as plt
import seaborn as sns
# 将数据集拆分为训练集和测试集
X_train,X_test,y_train,y_test = model_selection.train_test_split(X,y,test_size = 0.25, random_state=1)
# 构建伯努利贝叶斯分类器
bnb = naive_bayes.BernoulliNB()
# 模型在训练数据集上的拟合
bnb.fit(X_train,y_train)
# 模型在测试数据集上的预测
bnb_pred = bnb.predict(X_test)
# 构建混淆矩阵
cm = pd.crosstab(bnb_pred,y_test)
# 绘制混淆矩阵图
sns.heatmap(cm, annot = True, cmap = 'GnBu', fmt = 'd')
# 去除x轴和y轴标签
plt.xlabel('Real')
plt.ylabel('Predict')
# 显示图形
plt.show()
# 模型的预测准确率
print('模型的准确率为：\n',metrics.accuracy_score(y_test, bnb_pred))
print('模型的评估报告：\n',metrics.classification_report(y_test, bnb_pred))
# 准确率为84.64%，效果还不错
# 计算正例Positive所对应的概率，用于生成ROC曲线的数据
y_score = bnb.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test.map({'Negative':0,'Positive':1}), y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# ROC曲线值为0.92，说明模型效果还是非常理想的
13、SVM
超平面
线性可分，非线性可分
13.1、SVM原理
超平面的理解：
寻找超平面：
由于第三步不可能寻找所有的分割线，所以这只是一个原理
分隔带
分隔带越宽越好
目标函数：
函数间隔：
几何间隔：
目标函数（等价变换等内容见PPT）：
13.2、线性可分SVM
拉格朗日乘子法
基于拉格朗日乘子法的目标函数
13.3、非线性可分SVM
将不可以线性划分的空间提升维度，找到超平面，实现数据分割
13.4、SVM实战-分类-手写字母分类
# 导入第三方模块
from sklearn import svm
import pandas as pd
from sklearn import model_selection
from sklearn import metrics
# 读取外部数据
letters = pd.read_csv(r"D:\study\pylearn\svmdata\letterdata.csv")
# 数据前5行
letters.head()
# letter为因变量，其余均为字母的信息（宽高边际等）
# 将数据拆分为训练集和测试集
predictors = letters.columns[1:]
X_train,X_test,y_train,y_test = model_selection.train_test_split(letters[predictors], letters.letter,
test_size = 0.25, random_state = 1234)
# 1.使用线性SVM
# 使用网格搜索法，选择线性可分SVM“类”中的最佳C值
C=[0.05,0.1,0.5,1,2,5]
parameters = {'C':C}
grid_linear_svc = model_selection.GridSearchCV(estimator = svm.LinearSVC(),param_grid =parameters,scoring='accuracy',cv=5,verbose =1)
# 模型在训练数据集上的拟合
grid_linear_svc.fit(X_train,y_train)
# 返回交叉验证后的最佳参数值
grid_linear_svc.best_params_, grid_linear_svc.best_score_
# 经过5次交叉验证后，发现最佳惩罚系数是0.1，模型在训练集的准确率有69.18%
输出：
({'C': 0.1}, 0.6918)
# 模型在测试集上的预测
pred_linear_svc = grid_linear_svc.predict(X_test)
# 模型的预测准确率
metrics.accuracy_score(y_test, pred_linear_svc)
# 模型在测试集的准确率有71.58%，均不高
输出：
0.7158
# 2.使用非线性SVM
# 使用网格搜索法，选择非线性SVM“类”中的最佳C值
kernel=['rbf','linear','poly','sigmoid']
C=[0.1,0.5,1,2,5]
parameters = {'kernel':kernel,'C':C}
grid_svc = model_selection.GridSearchCV(estimator = svm.SVC(),param_grid =parameters,scoring='accuracy',cv=5,verbose =1)
# 模型在训练数据集上的拟合
grid_svc.fit(X_train,y_train)
# 返回交叉验证后的最佳参数值
grid_svc.best_params_, grid_svc.best_score_
# 经过5重交叉验证后，发现最佳惩罚系数为5，最佳核函数为径向基核函数，在训练集的准确率达到了95.17%，非常好
输出：
({'C': 5, 'kernel': 'rbf'}, 0.9516666666666665)
# 模型在测试集上的预测
pred_svc = grid_svc.predict(X_test)
# 模型的预测准确率
metrics.accuracy_score(y_test,pred_svc)
# 测试集的准确率也达到了95.96%，非常好
输出：
0.9596
13.5、SVM实战-预测-预测森林火灾
# 读取外部数据
forestfires = pd.read_csv(r"D:\study\pylearn\svmdata\forestfires.csv")
# 数据前5行
forestfires.head()
# area为因变量，其余变量为火灾发生的位置，时间，各项指标天气，气温，湿度，风力等信息
# 删除day变量(因为day可能大概率与结果无关)
forestfires.drop('day',axis = 1, inplace = True)
# 将月份作数值化处理(考虑到月份也可能是火灾的一个因素，所以保留，并且转换为数值类型)
forestfires.month = pd.factorize(forestfires.month)[0]
# 预览数据前5行
forestfires.head()
# 导入第三方模块
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import norm
# 绘制森林烧毁面积的直方图
sns.distplot(forestfires.area, bins = 50, kde = True, fit = norm, hist_kws = {'color':'steelblue'}, kde_kws = {'color':'red', 'label':'Kernel Density'}, fit_kws = {'color':'black','label':'Nomal', 'linestyle':'--'})
# 显示图例
plt.legend()
# 显示图形
plt.show()
# 数据呈现严重的右偏，所以不能直接使用该变量，一般都会将数据做对数处理
# 导入第三方模块
from sklearn import preprocessing
import numpy as np
from sklearn import neighbors
# 对area变量作对数变换
y = np.log1p(forestfires.area)
# 将X变量作标准化处理
predictors = forestfires.columns[:-1]
X = preprocessing.scale(forestfires[predictors])
# 重新绘制森林烧毁面积的直方图观察
sns.distplot(y, bins = 50, kde = True, fit = norm, hist_kws = {'color':'steelblue'}, kde_kws = {'color':'red', 'label':'Kernel Density'}, fit_kws = {'color':'black','label':'Nomal', 'linestyle':'--'})
# 显示图例
plt.legend()
# 显示图形
plt.show()
# 发现好很多
# 将数据拆分为训练集和测试集
X_train,X_test,y_train,y_test = model_selection.train_test_split(X, y, test_size = 0.25, random_state = 1234)
# 构建默认参数的SVM回归模型
svr = svm.SVR()
# 模型在训练数据集上的拟合
svr.fit(X_train,y_train)
# 模型在测试上的预测
pred_svr = svr.predict(X_test)
# 计算模型的MSE
metrics.mean_squared_error(y_test,pred_svr)
# MSE=1.93
输出：
1.9268192310372867
# 使用网格搜索法，选择SVM回归中的最佳C值、epsilon值和gamma值
epsilon = np.arange(0.1,1.5,0.2)
C= np.arange(100,1000,200)
gamma = np.arange(0.001,0.01,0.002)
parameters = {'epsilon':epsilon,'C':C,'gamma':gamma}
grid_svr = model_selection.GridSearchCV(estimator = svm.SVR(),param_grid =parameters,
scoring='neg_mean_squared_error',cv=5,verbose =1, n_jobs=2)
# 模型在训练数据集上的拟合
grid_svr.fit(X_train,y_train)
# 返回交叉验证后的最佳参数值
print(grid_svr.best_params_, grid_svr.best_score_)
# 经过4重交叉验证，最佳惩罚系数为300，最佳epsilon为1.1，最佳gamma为0.001，MSE=-1.995
输出：
{'C': 300, 'epsilon': 1.1000000000000003, 'gamma': 0.001} -1.9946668196316373
# 模型在测试集上的预测
pred_grid_svr = grid_svr.predict(X_test)
# 计算模型在测试集上的MSE值
metrics.mean_squared_error(y_test,pred_grid_svr)
# 测试集的MSE=1.75，相对降低了一点
输出：
1.7455012238826775
14、GBDT
14.1、提升树-Adaboost算法：
生成每一颗子树的时候都依赖上一颗数，而随机森林的每一颗子树是随机生成的
对于Adaboost算法而言，每一棵基础决策树都是基于前一棵基础决策树的分类结果对样本点设置不同的权重，如果在前一棵基础决策树中将某样本点预测错误，就会增大该样本点的权重，否则会相应降低样本点的权重，进而再构建下一棵基础决策树，更加关注权重大的样本点
所以，AdaBoost算法需要解决三大难题，即样本点的权重w_mi如何确定、基础决策树f(x)如何选择以及每一棵基础决策树所对应的权重α_m如何计算。
Adaboost应用-信用卡违约分类
# 导入第三方包
import pandas as pd
import matplotlib.pyplot as plt
# 读入数据
default = pd.read_excel(r'D:\study\pylearn\gbdtdata\default of credit card clients.xls')
# 数据集中是否违约的客户比例
# 为确保绘制的饼图为圆形，需执行如下代码
plt.axes(aspect = 'equal')
# 中文乱码和坐标轴负号的处理
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 统计客户是否违约的频数
counts = default.y.value_counts()
# 绘制饼图
plt.pie(x = counts, # 绘图数据
labels=pd.Series(counts.index).map({0:'不违约',1:'违约'}), # 添加文字标签
autopct='%.1f%%' # 设置百分比的格式，这里保留一位小数
)
# 显示图形
plt.show()
# 违约占22.1%，不违约占77.9%
# 将数据集拆分为训练集和测试集
# 导入第三方包
from sklearn import model_selection
from sklearn import ensemble
from sklearn import metrics
# 排除数据集中的ID变量和因变量，剩余的数据用作自变量X
X = default.drop(['ID','y'], axis = 1)
y = default.y
# 数据拆分
X_train,X_test,y_train,y_test = model_selection.train_test_split(X,y,test_size = 0.25, random_state = 1234)
# 构建AdaBoost算法的类
AdaBoost1 = ensemble.AdaBoostClassifier()
# 算法在训练数据集上的拟合
AdaBoost1.fit(X_train,y_train)
# 算法在测试数据集上的预测
pred1 = AdaBoost1.predict(X_test)
# 返回模型的预测效果
print('模型的准确率为：\n',metrics.accuracy_score(y_test, pred1))
print('模型的评估报告：\n',metrics.classification_report(y_test, pred1))
# 得到模型准确率为81.25%，违约1的准确率为68%，不违约0的准确率为83%
# 计算客户违约的概率值，用于生成ROC曲线的数据
y_score = AdaBoost1.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# ROC曲线下的面积为0.78，接近0.8,效果不是很好
# 接下来，通过交叉验证，选择合理的参数值
# 自变量的重要性排序
importance = pd.Series(AdaBoost1.feature_importances_, index = X.columns)
importance.sort_values().plot(kind = 'barh')
plt.show()
# 查看重要的变量
# 取出重要性比较高的自变量建模
predictors = list(importance[importance>0.02].index)
predictors
# 通过网格搜索法选择基础模型所对应的合理参数组合
# 导入第三方包
from sklearn.model_selection import GridSearchCV
from sklearn.tree import DecisionTreeClassifier
max_depth = [3,4,5,6]
params1 = {'base_estimator__max_depth':max_depth}
base_model = GridSearchCV(estimator = ensemble.AdaBoostClassifier(base_estimator = DecisionTreeClassifier()),
param_grid= params1, scoring = 'roc_auc', cv = 5, n_jobs = 4, verbose = 1)
base_model.fit(X_train[predictors],y_train)
# 返回参数的最佳组合和对应AUC值
base_model.best_params_, base_model.best_score_
# 经过5重交叉验证，对于基础模型CART来说，最大树深度为3层比较合理，预测准确率为74.13%
输出：
({'base_estimator__max_depth': 3}, 0.7413109515305258)
# 通过网格搜索法选择提升树的合理参数组合
# 导入第三方包
from sklearn.model_selection import GridSearchCV
n_estimators = [100,200,300]
learning_rate = [0.01,0.05,0.1,0.2]
params2 = {'n_estimators':n_estimators,'learning_rate':learning_rate}
adaboost = GridSearchCV(estimator = ensemble.AdaBoostClassifier(base_estimator = DecisionTreeClassifier(max_depth = 3)),
param_grid= params2, scoring = 'roc_auc', cv = 5, n_jobs = 4, verbose = 1)
adaboost.fit(X_train[predictors] ,y_train)
# 返回参数的最佳组合和对应AUC值
adaboost.best_params_, adaboost.best_score_
# 进过5重交叉验证后，得知AdaBoost的最佳基础模型是300个，学习率为0.05,准确率为76.94%
输出：
({'learning_rate': 0.05, 'n_estimators': 100}, 0.7693742071267697)
# 使用最佳的参数组合构建AdaBoost模型
AdaBoost2 = ensemble.AdaBoostClassifier(base_estimator = DecisionTreeClassifier(max_depth = 3),
n_estimators = 300, learning_rate = 0.01)
# 算法在训练数据集上的拟合
AdaBoost2.fit(X_train[predictors],y_train)
# 算法在测试数据集上的预测
pred2 = AdaBoost2.predict(X_test[predictors])
# 返回模型的预测效果
print('模型的准确率为：\n',metrics.accuracy_score(y_test, pred2))
print('模型的评估报告：\n',metrics.classification_report(y_test, pred2))
# 模型准确率为81.29%，对比原先的81.25，仅提高了0.04%
# 计算正例的预测概率，用于生成ROC曲线的数据
y_score = AdaBoost2.predict_proba(X_test[predictors])[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# ROC的面积还是0.78
14.2、梯度提升树-GBDT算法
梯度提升树算法实际上是提升算法的扩展版，在原始的提升算法中，如果损失函数为平方损失或指数损失，求解损失函数的最小值问题会非常简单，但如果损失函数为更一般的函数（如绝对值损失函数或Huber损失函数等），目标值的求解就会相对复杂很多。
梯度提升算法，是在第m轮基础模型中，利用损失函数的负梯度值作为该轮基础模型损失值的近似，并利用这个近似值构建下一轮基础模型。利用损失函数的负梯度值近似残差的计算就是梯度提升算法在提升算法上的扩展，这样的扩展使得目标函数的求解更为方便。
信用卡违约GBDT
# 运用网格搜索法选择梯度提升树的合理参数组合
learning_rate = [0.01,0.05,0.1,0.2]
n_estimators = [100,300,500]
max_depth = [3,4,5,6]
params = {'learning_rate':learning_rate,'n_estimators':n_estimators,'max_depth':max_depth}
gbdt_grid = GridSearchCV(estimator = ensemble.GradientBoostingClassifier(),
param_grid= params, scoring = 'roc_auc', cv = 5, n_jobs = 4, verbose = 1)
gbdt_grid.fit(X_train[predictors],y_train)
# 返回参数的最佳组合和对应AUC值
gbdt_grid.best_params_, gbdt_grid.best_score_
输出：
({'learning_rate': 0.1, 'max_depth': 4, 'n_estimators': 100},
0.7746592679803446)
# 基于最佳参数组合的GBDT模型，对测试数据集进行预测
pred = gbdt_grid.predict(X_test[predictors])
# 返回模型的预测效果
print('模型的准确率为：\n',metrics.accuracy_score(y_test, pred))
print('模型的评估报告：\n',metrics.classification_report(y_test, pred))
# 计算违约客户的概率值，用于生成ROC曲线的数据
y_score = gbdt_grid.predict_proba(X_test[predictors])[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
14.3、GBDT的改进XGBoost
XGBoost是由传统的GBDT模型发展而来的，GBDT模型在求解最优化问题时应用了一阶导技术，而XGBoost则使用损失函数的一阶和二阶导，而且可以自定义损失函数，只要损失函数可一阶和二阶求导。
XGBoost算法相比于GBDT算法还有其他优点，例如支持并行计算，大大提高算法的运行效率；XGBoost在损失函数中加入了正则项，用来控制模型的复杂度，进而可以防止模型的过拟合；XGBoost除了支持CART基础模型，还支持线性基础模型；XGBoost采用了随机森林的思想，对字段进行抽样，既可以防止过拟合，也可以降低模型的计算量。
信用卡欺诈数据
# 读入数据
creditcard = pd.read_csv(r"D:\study\pylearn\gbdtdata\creditcard.csv")
# 为确保绘制的饼图为圆形，需执行如下代码
plt.axes(aspect = 'equal')
# 统计交易是否为欺诈的频数
counts = creditcard.Class.value_counts()
# 绘制饼图
plt.pie(x = counts, # 绘图数据
labels=pd.Series(counts.index).map({0:'正常',1:'欺诈'}), # 添加文字标签
autopct='%.2f%%' # 设置百分比的格式，这里保留一位小数
)
# 显示图形
plt.show()
# 欺诈数据仅占0.17%，数据比例严重失衡，因为预测非欺诈数据几乎为99%，预测为欺诈数据为0%左右
# 所以通过SMOTE算法转换为相对平衡的数据
# 将数据拆分为训练集和测试集
# 删除自变量中的Time变量
X = creditcard.drop(['Time','Class'], axis = 1)
y = creditcard.Class
# 数据拆分
X_train,X_test,y_train,y_test = model_selection.train_test_split(X,y,test_size = 0.3, random_state = 1234)
# 导入第三方包
from imblearn.over_sampling import SMOTE
# 运用SMOTE算法实现训练数据集的平衡
over_samples = SMOTE(random_state=1234)
over_samples_X,over_samples_y = over_samples.fit_resample(X_train, y_train)
#over_samples_X, over_samples_y = over_samples.fit_sample(X_train.values,y_train.values.ravel())
# 重抽样前的类别比例
print(y_train.value_counts()/len(y_train))
# 重抽样后的类别比例
print(pd.Series(over_samples_y).value_counts()/len(over_samples_y))
# 使用SMOTE算法重抽样后，比例由0.998:0.002达到平衡0.5:0.5
输出：
0 0.998239
1 0.001761
Name: Class, dtype: float64
0 0.5
1 0.5
Name: Class, dtype: float64
安装xgboost
下载xgboost版本：https://www.lfd.uci.edu/~gohlke/pythonlibs/#xgboost，
选择对应的python版本下载，然后
pip install xgboost-1.5.1-cp37-cp37m-win32.whl
# 导入第三方包
# 安装xgboost方法：https://zhuanlan.zhihu.com/p/55394252
import xgboost
import numpy as np
# 构建XGBoost分类器
xgboost1 = xgboost.XGBClassifier()
# 使用重抽样后的数据，对其建模
xgboost1.fit(over_samples_X,over_samples_y)
# 将模型运用到测试数据集中
resample_pred = xgboost1.predict(np.array(X_test))
# 返回模型的预测效果
print('模型的准确率为：\n',metrics.accuracy_score(y_test, resample_pred))
print('模型的评估报告：\n',metrics.classification_report(y_test, resample_pred))
# 预测准确率超过99%，欺诈预测准确率为81%，非欺诈预测100%
# 计算欺诈交易的概率值，用于生成ROC曲线的数据
y_score = xgboost1.predict_proba(np.array(X_test))[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# ROC曲线面积高达97%，非常good
# 在非平衡数据（抽样前数据）的预测效果
# 构建XGBoost分类器
xgboost2 = xgboost.XGBClassifier()
# 使用非平衡的训练数据集拟合模型
xgboost2.fit(X_train,y_train)
# 基于拟合的模型对测试数据集进行预测
pred2 = xgboost2.predict(X_test)
# 混淆矩阵
pd.crosstab(pred2,y_test)
# 返回模型的预测效果
print('模型的准确率为：\n',metrics.accuracy_score(y_test, pred2))
print('模型的评估报告：\n',metrics.classification_report(y_test, pred2))
# 准确率依旧99%
# 计算欺诈交易的概率值，用于生成ROC曲线的数据
y_score = xgboost2.predict_proba(X_test)[:,1]
fpr,tpr,threshold = metrics.roc_curve(y_test, y_score)
# 计算AUC的值
roc_auc = metrics.auc(fpr,tpr)
# 绘制面积图
plt.stackplot(fpr, tpr, color='steelblue', alpha = 0.5, edgecolor = 'black')
# 添加边际线
plt.plot(fpr, tpr, color='black', lw = 1)
# 添加对角线
plt.plot([0,1],[0,1], color = 'red', linestyle = '--')
# 添加文本信息
plt.text(0.5,0.3,'ROC curve (area = %0.2f)' % roc_auc)
# 添加x轴与y轴标签
plt.xlabel('1-Specificity')
plt.ylabel('Sensitivity')
# 显示图形
plt.show()
# ROC曲线面积为0.97
15、Kmeans聚类
K均值聚类
15.1、Kmeans原理
原理
- 从数据中随机挑选k个样本点作为原始的簇中心
- 计算剩余样本与簇中心的距离，并把各样本标记为离k个簇中心最近的类别
- 重新计算各簇中样本点的均值，并以均值作为新的k个簇中心
- 不断重复第二步和第三步，直到簇中心的变化趋于稳定，形成最终的k个簇
15.2、最佳k值的确定
（1）拐点法找k值
簇内离差平方和拐点法的思想很简单，就是在不同的k值下计算簇内离差平方和，然后通过可视化的方法找到“拐点”所对应的k值。当折线图中的斜率由大突然变小时，并且之后的斜率变化缓慢，则认为突然变化的点就是寻找的目标点，因为继续随着簇数k的增加，聚类效果不再有大的变化。
# 随机生成数据，并绘制散点图查看
# 导入第三方包
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn import metrics
# 随机生成三组二元正态分布随机数
np.random.seed(1234)
mean1 = [0.5, 0.5]
cov1 = [[0.3, 0], [0, 0.3]]
x1, y1 = np.random.multivariate_normal(mean1, cov1, 1000).T
mean2 = [0, 8]
cov2 = [[1.5, 0], [0, 1]]
x2, y2 = np.random.multivariate_normal(mean2, cov2, 1000).T
mean3 = [8, 4]
cov3 = [[1.5, 0], [0, 1]]
x3, y3 = np.random.multivariate_normal(mean3, cov3, 1000).T
# 绘制三组数据的散点图
plt.scatter(x1,y1)
plt.scatter(x2,y2)
plt.scatter(x3,y3)
# 显示图形
plt.show()
# 虚拟数据呈现3个簇
# 使用拐点法，绘制簇的个数与总的簇内离差平方和之间的折线图
# 构造自定义函数，用于绘制不同k值和对应总的簇内离差平方和的折线图
def k_SSE(X, clusters):
# 选择连续的K种不同的值
K = range(1,clusters+1)
# 构建空列表用于存储总的簇内离差平方和
TSSE = []
for k in K:
# 用于存储各个簇内离差平方和
SSE = []
kmeans = KMeans(n_clusters=k)
kmeans.fit(X)
# 返回簇标签
labels = kmeans.labels_
# 返回簇中心
centers = kmeans.cluster_centers_
# 计算各簇样本的离差平方和，并保存到列表中
for label in set(labels):
SSE.append(np.sum((X.loc[labels == label,]-centers[label,:])**2))
# 计算总的簇内离差平方和
TSSE.append(np.sum(SSE))
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
# 绘制K的个数与GSSE的关系
plt.plot(K, TSSE, 'b*-')
plt.xlabel('簇的个数')
plt.ylabel('簇内离差平方和之和')
# 显示图形
plt.show()
# 将三组数据集汇总到数据框中
X = pd.DataFrame(np.concatenate([np.array([x1,y1]),np.array([x2,y2]),np.array([x3,y3])], axis = 1).T)
# 自定义函数的调用
k_SSE(X, 15)
# 图表明，当簇的个数为3时，形成了一个很明显的拐点，所以合理的k=3
（2）轮廓系数法
该方法综合考虑了簇的密集性与分散性两个信息，如果数据集被分割为理想的k个簇，那么对应的簇内样本会很密集，而簇间样本会很分散。轮廓系数的计算公式可以表示为：
其中，a(i)体现了簇内的密集性，代表样本i与同簇内其他样本点距离的平均值；b(i)反映了簇间的分散性，它的计算过程是，样本i与其他非同簇样本点距离的平均值，然后从平均值中挑选出最小值。
当S(i)接近于-1时，说明样本i分配的不合理，需要将其分配到其他簇中；当S(i)近似为0时，说明样本i落在了模糊地带，即簇的边界处；当S(i)近似为1时，说明样本i的分配是合理的。
# 使用轮廓系数法选择最佳k值
# 构造自定义函数，用于绘制不同k值和对应轮廓系数的折线图
def k_silhouette(X, clusters):
K = range(2,clusters+1)
# 构建空列表，用于存储个中簇数下的轮廓系数
S = []
for k in K:
kmeans = KMeans(n_clusters=k)
kmeans.fit(X)
labels = kmeans.labels_
# 调用字模块metrics中的silhouette_score函数，计算轮廓系数
S.append(metrics.silhouette_score(X, labels, metric='euclidean'))
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
# 绘制K的个数与轮廓系数的关系
plt.plot(K, S, 'b*-')
plt.xlabel('簇的个数')
plt.ylabel('轮廓系数')
# 显示图形
plt.show()
# 自定义函数的调用
k_silhouette(X, 15)
#图表示，当k=3时，轮廓系数最大，且比较接近1，所以最佳k=3
（3）间隔统计量法
# 自定义函数，计算簇内任意两样本之间的欧氏距离
def short_pair_wise_D(each_cluster):
mu = each_cluster.mean(axis = 0)
Dk = sum(sum((each_cluster - mu)**2)) * 2.0 * each_cluster.shape[0]
return Dk
# 计算簇内的Wk值
def compute_Wk(data, classfication_result):
Wk = 0
label_set = set(classfication_result)
for label in label_set:
each_cluster = data[classfication_result == label, :]
Wk = Wk + short_pair_wise_D(each_cluster)/(2.0*each_cluster.shape[0])
return Wk
# 计算GAP统计量
def gap_statistic(X, B=10, K=range(1,11), N_init = 10):
# 将输入数据集转换为数组
X = np.array(X)
# 生成B组参照数据
shape = X.shape
tops = X.max(axis=0)
bots = X.min(axis=0)
dists = np.matrix(np.diag(tops-bots))
rands = np.random.random_sample(size=(B,shape[0],shape[1]))
for i in range(B):
rands[i,:,:] = rands[i,:,:]*dists+bots
# 自定义0元素的数组，用于存储gaps、Wks和Wkbs
gaps = np.zeros(len(K))
Wks = np.zeros(len(K))
Wkbs = np.zeros((len(K),B))
# 循环不同的k值，
for idxk, k in enumerate(K):
k_means = KMeans(n_clusters=k)
k_means.fit(X)
classfication_result = k_means.labels_
# 将所有簇内的Wk存储起来
Wks[idxk] = compute_Wk(X,classfication_result)
# 通过循环，计算每一个参照数据集下的各簇Wk值
for i in range(B):
Xb = rands[i,:,:]
k_means.fit(Xb)
classfication_result_b = k_means.labels_
Wkbs[idxk,i] = compute_Wk(Xb,classfication_result_b)
# 计算gaps、sd_ks、sk和gapDiff
gaps = (np.log(Wkbs)).mean(axis = 1) - np.log(Wks)
sd_ks = np.std(np.log(Wkbs), axis=1)
sk = sd_ks*np.sqrt(1+1.0/B)
# 用于判别最佳k的标准，当gapDiff首次为正时，对应的k即为目标值
gapDiff = gaps[:-1] - gaps[1:] + sk[1:]
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
# 绘制gapDiff的条形图
plt.bar(np.arange(len(gapDiff))+1, gapDiff, color = 'steelblue')
plt.xlabel('簇的个数')
plt.ylabel('k的选择标准')
plt.show()
# 自定义函数的调用
gap_statistic(X)
# 图表示，x轴代表了不同的簇数k,y轴代表了k值选择的判断指标gapDiff，gapDiff首次出现正值对应的为k=3,所以最佳k=3
15.3、kmeans实战-iris鸢尾花分类
已知聚类k值为3
# iris数据集包含150个观测，3种花，每种花50个样本，已知聚类为3
# 读取iris数据集
iris = pd.read_csv(r"D:\study\pylearn\kmeansdata\iris.csv")
# 查看数据集的前几行
iris.head()
# 提取出用于建模的数据集X
X = iris.drop(labels = 'Species', axis = 1)
# 构建Kmeans模型
kmeans = KMeans(n_clusters = 3)
kmeans.fit(X)
# 聚类结果标签
X['cluster'] = kmeans.labels_
# 各类频数统计
X.cluster.value_counts()
# 各个簇的样本量分别为62,50,38
输出：
0 62
1 50
2 38
Name: cluster, dtype: int64
# 导入第三方模块
import seaborn as sns
# 三个簇的簇中心
centers = kmeans.cluster_centers_
# 绘制聚类效果的散点图
sns.lmplot(x = 'Petal_Length', y = 'Petal_Width', hue = 'cluster', markers = ['^','s','o'],
data = X, fit_reg = False, scatter_kws = {'alpha':0.8}, legend_out = False)
plt.scatter(centers[:,2], centers[:,3], marker = '*', color = 'black', s = 130)
plt.xlabel('花瓣长度')
plt.ylabel('花瓣宽度')
# 图形显示
plt.show()
# 建模后的三类如下
# 增加一个辅助列，将不同的花种映射到0,1,2三种值，目的方便后面图形的对比
iris['Species_map'] = iris.Species.map({'virginica':0,'setosa':1,'versicolor':2})
# 绘制原始数据三个类别的散点图
sns.lmplot(x = 'Petal_Length', y = 'Petal_Width', hue = 'Species_map', data = iris, markers = ['^','s','o'],
fit_reg = False, scatter_kws = {'alpha':0.8}, legend_out = False)
plt.xlabel('花瓣长度')
plt.ylabel('花瓣宽度')
# 图形显示
plt.show()
# 原始数据的三类如下
# 导入第三方模块
import pygal
# 调用Radar这个类，并设置雷达图的填充，及数据范围
radar_chart = pygal.Radar(fill = True)
# 添加雷达图各顶点的名称
radar_chart.x_labels = ['花萼长度','花萼宽度','花瓣长度','花瓣宽度']
# 绘制三个雷达图区域，代表三个簇中心的指标值
radar_chart.add('C1', centers[0])
radar_chart.add('C2', centers[1])
radar_chart.add('C3', centers[2])
# 保存图像
radar_chart.render_to_file('radar_chart.svg')
# C1类花萼长度，花瓣长度，宽度都是最大的，C2类三个值都是最小的，C3类都是中值
15.4、kmeans实战-NBA球员聚类
iris是知道k值的聚类，而NBA球员不清楚k值，需要求最佳k值
# 读取球员数据
players = pd.read_csv(r"D:\study\pylearn\kmeansdata\players.csv")
players.head()
# 绘制得分与命中率的散点图
sns.lmplot(x = '得分', y = '命中率', data = players,
fit_reg = False, scatter_kws = {'alpha':0.8, 'color': 'steelblue'})
plt.show()
# 选择最佳k值-拐点法
from sklearn import preprocessing
# 数据标准化处理
X = preprocessing.minmax_scale(players[['得分','罚球命中率','命中率','三分命中率']])
# 将数组转换为数据框
X = pd.DataFrame(X, columns=['得分','罚球命中率','命中率','三分命中率'])
# 使用拐点法选择最佳的K值
k_SSE(X, 15)
# 得到最佳k=3(可能在3,4,5之间)
# 选择最佳k值-使用轮廓系数选择最佳的K值
k_silhouette(X, 15)
# k=2的时候轮廓系数最大，最佳k=2
# 选择最佳k值-使用间隙统计量选择最佳的K值
gap_statistic(X, B = 20, K=range(1, 16))
# 横坐标k=3的时候纵坐标首次为正值，所以最终确定k=3
# 将球员数据集聚为3类
kmeans = KMeans(n_clusters = 3)
kmeans.fit(X)
# 将聚类结果标签插入到数据集players中
players['cluster'] = kmeans.labels_
# 构建空列表，用于存储三个簇的簇中心
centers = []
for i in players.cluster.unique():
centers.append(players.loc[players.cluster == i,['得分','罚球命中率','命中率','三分命中率']].mean())
# 将列表转换为数组，便于后面的索引取数
centers = np.array(centers)
# 绘制散点图
sns.lmplot(x = '得分', y = '命中率', hue = 'cluster', data = players, markers = ['^','s','o'],
fit_reg = False, scatter_kws = {'alpha':0.8}, legend = False)
# 添加簇中心
plt.scatter(centers[:,0], centers[:,2], c='k', marker = '*', s = 180)
plt.xlabel('得分')
plt.ylabel('命中率')
# 图形显示
plt.show()
# 雷达图
# 调用模型计算出来的簇中心
centers_std = kmeans.cluster_centers_
# 设置填充型雷达图
radar_chart = pygal.Radar(fill = True)
# 添加雷达图各顶点的名称
radar_chart.x_labels = ['得分','罚球命中率','命中率','三分命中率']
# 绘制雷达图代表三个簇中心的指标值
radar_chart.add('C1', centers_std[0])
radar_chart.add('C2', centers_std[1])
radar_chart.add('C3', centers_std[2])
# 保存图像
radar_chart.render_to_file('radar_chart2.svg')
# 比较各个维度的差异
16、DBSCAN与层次聚类
DBSCAN 密度聚类
16.1、原理
（1）为密度聚类算法设置一个合理的半径e以及e领域内所包含的最少样本量MinPts。
（2）从数据集中随机挑选一个样本点p，检验其在e领域内是否包含指定的最少样本量，如果包含就将其定性为核心对象，并构成一个簇C；否则，重新挑选一个样本点。
（3）对于核心对象p所覆盖的其他样本点q，如果点q对应的e领域内仍然包含最少样本量MinPts，就将其覆盖的样本点统统归于簇C。
（4）重复步骤（3），将最大的密度相连所包含的样本点聚为一类，形成一个大簇。
（5）完成步骤（4）后，重新回到步骤（2），并重复步骤（3）和（4），直到没有新的样本点可以生成新簇时算法结束。
16.2、密度聚类相比Kmeans的优势
Kmeans聚类存在两个致命缺点，一是聚类效果容易受到异常样本点的影响；二是该算法无法准确地将非球形样本进行合理的聚类。
基于密度的聚类则可以解决非球形簇的问题，“密度”可以理解为样本点的紧密程度，如果在指定的半径领域内，实际样本量超过给定的最小样本量阈值，则认为是密度高的对象，就可以聚成一个簇。
Kmeans聚类-对球形簇的情况（很好）
Kmeans聚类-对非球形簇的情况（出现了异常）
通过密度聚类-解决非球形簇的情况
比较测试
# 导入第三方模块
import pandas as pd
import numpy as np
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import cluster
# 模拟数据集
X,y = make_blobs(n_samples = 2000, centers = [[-1,-2],[1,3]], cluster_std = [0.5,0.5], random_state = 1234)
# 将模拟得到的数组转换为数据框，用于绘图
plot_data = pd.DataFrame(np.column_stack((X,y)), columns = ['x1','x2','y'])
# 设置绘图风格
plt.style.use('ggplot')
# 绘制散点图（用不同的形状代表不同的簇）
sns.lmplot('x1', 'x2', data = plot_data, hue = 'y',markers = ['^','o'],
fit_reg = False, legend = False)
# 显示图形
plt.show()
# 导入第三方模块
from sklearn import cluster
# 构建Kmeans聚类和密度聚类
kmeans = cluster.KMeans(n_clusters=2, random_state=1234)
kmeans.fit(X)
dbscan = cluster.DBSCAN(eps = 0.5, min_samples = 10)
dbscan.fit(X)
# 将Kmeans聚类和密度聚类的簇标签添加到数据框中
plot_data['kmeans_label'] = kmeans.labels_
plot_data['dbscan_label'] = dbscan.labels_
# 绘制聚类效果图
# 设置大图框的长和高
plt.figure(figsize = (12,6))
# 设置第一个子图的布局
ax1 = plt.subplot2grid(shape = (1,2), loc = (0,0))
# 绘制散点图
ax1.scatter(plot_data.x1, plot_data.x2, c = plot_data.kmeans_label)
# 设置第二个子图的布局
ax2 = plt.subplot2grid(shape = (1,2), loc = (0,1))
# 绘制散点图(为了使Kmeans聚类和密度聚类的效果图颜色一致，通过序列的map“方法”对颜色作重映射)
ax2.scatter(plot_data.x1, plot_data.x2, c=plot_data.dbscan_label.map({-1:1,0:2,1:0}))
# 显示图形
plt.show()
# Kmeans聚类（左），密度聚类（右），二者都可以很好的分类，密度聚类不仅得到了很好的效果，还发现了原理簇的异常点
# 导入第三方模块
from sklearn.datasets import make_moons
# 构造非球形样本点
X1,y1 = make_moons(n_samples=2000, noise = 0.05, random_state = 1234)
# 构造球形样本点
X2,y2 = make_blobs(n_samples=1000, centers = [[3,3]], cluster_std = 0.5, random_state = 1234)
# 将y2的值替换为2(为了避免与y1的值冲突，因为原始y1和y2中都有0这个值)
y2 = np.where(y2 == 0,2,0)
# 将模拟得到的数组转换为数据框，用于绘图
plot_data = pd.DataFrame(np.row_stack([np.column_stack((X1,y1)),np.column_stack((X2,y2))]), columns = ['x1','x2','y'])
# 绘制散点图（用不同的形状代表不同的簇）
sns.lmplot('x1', 'x2', data = plot_data, hue = 'y',markers = ['^','o','>'],
fit_reg = False, legend = False)
# 显示图形
plt.show()
# 三组数据，一组球形，两组非球形
# 构建Kmeans聚类和密度聚类
kmeans = cluster.KMeans(n_clusters=3, random_state=1234)
kmeans.fit(plot_data[['x1','x2']])
dbscan = cluster.DBSCAN(eps = 0.3, min_samples = 5)
dbscan.fit(plot_data[['x1','x2']])
# 将Kmeans聚类和密度聚类的簇标签添加到数据框中
plot_data['kmeans_label'] = kmeans.labels_
plot_data['dbscan_label'] = dbscan.labels_
# 绘制聚类效果图
# 设置大图框的长和高
plt.figure(figsize = (12,6))
# 设置第一个子图的布局
ax1 = plt.subplot2grid(shape = (1,2), loc = (0,0))
# 绘制散点图
ax1.scatter(plot_data.x1, plot_data.x2, c = plot_data.kmeans_label)
# 设置第二个子图的布局
ax2 = plt.subplot2grid(shape = (1,2), loc = (0,1))
# 绘制散点图(为了使Kmeans聚类和密度聚类的效果图颜色一致，通过序列的map“方法”对颜色作重映射)
ax2.scatter(plot_data.x1, plot_data.x2, c=plot_data.dbscan_label.map({-1:1,0:0,1:3,2:2}))
# 显示图形
plt.show()
# Kmeans聚类（左），密度聚类（右），可以看出，Kmeans聚类效果很差，而密度聚类效果就很好
层次聚类
# 构造两个球形簇的数据样本点
X,y = make_blobs(n_samples = 2000, centers = [[-1,0],[1,0.5]], cluster_std = [0.2,0.45], random_state = 1234)
# 将模拟得到的数组转换为数据框，用于绘图
plot_data = pd.DataFrame(np.column_stack((X,y)), columns = ['x1','x2','y'])
# 绘制散点图（用不同的形状代表不同的簇）
sns.lmplot('x1', 'x2', data = plot_data, hue = 'y',markers = ['^','o'],
fit_reg = False, legend = False)
# 显示图形
plt.show()
# 两个球形簇
# 设置大图框的长和高
plt.figure(figsize = (16,5))
# 设置第一个子图的布局
ax1 = plt.subplot2grid(shape = (1,3), loc = (0,0))
# 层次聚类--最小距离法
agnes_min = cluster.AgglomerativeClustering(n_clusters = 2, linkage='ward')
agnes_min.fit(X)
# 绘制聚类效果图
ax1.scatter(X[:,0], X[:,1], c=agnes_min.labels_)
# 设置第二个子图的布局
ax2 = plt.subplot2grid(shape = (1,3), loc = (0,1))
# 层次聚类--最大距离法
agnes_max = cluster.AgglomerativeClustering(n_clusters = 2, linkage='complete')
agnes_max.fit(X)
ax2.scatter(X[:,0], X[:,1], c=agnes_max.labels_)
# 设置第三个子图的布局
ax2 = plt.subplot2grid(shape = (1,3), loc = (0,2))
# 层次聚类--平均距离法
agnes_avg = cluster.AgglomerativeClustering(n_clusters = 2, linkage='average')
agnes_avg.fit(X)
plt.scatter(X[:,0], X[:,1], c=agnes_avg.labels_)
plt.show()
# 三种层次聚类效果比较，最小距离（左），最大距离（中），平均距离（右），最小距离和平均距离均表现了较好的效果，而最大距离效果很差，主要是由于模糊地带（原始数据中两个簇交界的区域）的异常点夸大了簇之间的距离
16.3、密度聚类实战-出生死亡率分类
# 读取外部数据
Province = pd.read_excel(r"D:\study\pylearn\dascandata\Province.xlsx")
Province.head()
# 绘制出生率与死亡率散点图
plt.scatter(Province.Birth_Rate, Province.Death_Rate, c = 'steelblue')
# 添加轴标签
plt.xlabel('Birth_Rate')
plt.ylabel('Death_Rate')
# 显示图形
plt.show()
# 肉眼近似发现3个簇
# 读入第三方包
from sklearn import preprocessing
# 选取建模的变量
predictors = ['Birth_Rate','Death_Rate']
# 变量的标准化处理
X = preprocessing.scale(Province[predictors])
X = pd.DataFrame(X)
# 构建空列表，用于保存不同参数组合下的结果
res = []
# 迭代不同的eps值
for eps in np.arange(0.001,1,0.05):
# 迭代不同的min_samples值
for min_samples in range(2,10):
dbscan = cluster.DBSCAN(eps = eps, min_samples = min_samples)
# 模型拟合
dbscan.fit(X)
# 统计各参数组合下的聚类个数（-1表示异常点）
n_clusters = len([i for i in set(dbscan.labels_) if i != -1])
# 异常点的个数
outliners = np.sum(np.where(dbscan.labels_ == -1, 1,0))
# 统计每个簇的样本个数
stats = str(pd.Series([i for i in dbscan.labels_ if i != -1]).value_counts().values)
res.append({'eps':eps,'min_samples':min_samples,'n_clusters':n_clusters,'outliners':outliners,'stats':stats})
# 将迭代后的结果存储到数据框中
df = pd.DataFrame(res)
# 根据条件筛选合理的参数组合
df.loc[df.n_clusters == 3, :]
# 利用上述的参数组合值，重建密度聚类算法
dbscan = cluster.DBSCAN(eps = 0.801, min_samples = 3)
# 模型拟合
dbscan.fit(X)
Province['dbscan_label'] = dbscan.labels_
# 绘制聚类聚类的效果散点图
sns.lmplot(x = 'Birth_Rate', y = 'Death_Rate', hue = 'dbscan_label', data = Province,
markers = ['*','d','^','o'], fit_reg = False, legend = False)
# 添加省份标签
for x,y,text in zip(Province.Birth_Rate,Province.Death_Rate, Province.Province):
plt.text(x+0.1,y-0.1,text, size = 8)
# 添加参考线
plt.hlines(y = 5.8, xmin = Province.Birth_Rate.min(), xmax = Province.Birth_Rate.max(),
linestyles = '--', colors = 'red')
plt.vlines(x = 10, ymin = Province.Death_Rate.min(), ymax = Province.Death_Rate.max(),
linestyles = '--', colors = 'red')
# 添加轴标签
plt.xlabel('Birth_Rate')
plt.ylabel('Death_Rate')
# 显示图形
plt.show()
# 利用最小距离法构建层次聚类
agnes_min = cluster.AgglomerativeClustering(n_clusters = 3, linkage='ward')
# 模型拟合
agnes_min.fit(X)
Province['agnes_label'] = agnes_min.labels_
# 绘制层次聚类的效果散点图
sns.lmplot(x = 'Birth_Rate', y = 'Death_Rate', hue = 'agnes_label', data = Province,
markers = ['d','^','o'], fit_reg = False, legend = False)
# 添加轴标签
plt.xlabel('Birth_Rate')
plt.ylabel('Death_Rate')
# 显示图形
plt.show()
# 导入第三方模块
from sklearn import metrics
# 构造自定义函数，用于绘制不同k值和对应轮廓系数的折线图
def k_silhouette(X, clusters):
K = range(2,clusters+1)
# 构建空列表，用于存储个中簇数下的轮廓系数
S = []
for k in K:
kmeans = cluster.KMeans(n_clusters=k)
kmeans.fit(X)
labels = kmeans.labels_
# 调用字模块metrics中的silhouette_score函数，计算轮廓系数
S.append(metrics.silhouette_score(X, labels, metric='euclidean'))
# 中文和负号的正常显示
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
# 设置绘图风格
plt.style.use('ggplot')
# 绘制K的个数与轮廓系数的关系
plt.plot(K, S, 'b*-')
plt.xlabel('簇的个数')
plt.ylabel('轮廓系数')
# 显示图形
plt.show()
# 聚类个数的探索
k_silhouette(X, clusters = 10)
# 当簇的个数为3时，轮廓系数达到最大，说明聚类为3较为合理
# 利用Kmeans聚类
kmeans = cluster.KMeans(n_clusters = 3)
# 模型拟合
kmeans.fit(X)
Province['kmeans_label'] = kmeans.labels_
# 绘制Kmeans聚类的效果散点图
sns.lmplot(x = 'Birth_Rate', y = 'Death_Rate', hue = 'kmeans_label', data = Province,
markers = ['d','^','o'], fit_reg = False, legend = False)
# 添加轴标签
plt.xlabel('Birth_Rate')
plt.ylabel('Death_Rate')
plt.show()
# 分为了三类
评论区