标签归档:数据分析

数据分析实战教程|Pandas处理数据太慢,来试试Polars吧!

很多人在学习数据分析的时候,肯定都会用到Pandas这个库,非常的实用!

从创建数据到读取各种格式的文件(text、csv、json),或者对数据进行切片和分割组合多个数据源,Pandas都能够很好的满足。

Pandas最初发布于2008年,使用Python、Cython和C编写的。是一个超级强大、快速易于使用的Python库,用于数据分析和处理。

当然Pandas也是有不足之处的,比如不具备多处理器,处理较大的数据集速度很慢。

今天,小F就给大家介绍一个新兴的Python库——Polars。

使用语法和Pandas差不多,处理数据的速度却比Pandas快了不少。

一个是大熊猫,一个是北极熊~

GitHub地址:https://github.com/ritchie46/polars

使用文档:https://ritchie46.github.io/polars-book/

Polars是通过Rust编写的一个库,Polars的内存模型是基于Apache Arrow。

Polars存在两种API,一种是Eager API,另一种则是Lazy API。

其中Eager API和Pandas的使用类似,语法差不太多,立即执行就能产生结果。

而Lazy API就像Spark,首先将查询转换为逻辑计划,然后对计划进行重组优化,以减少执行时间和内存使用。

安装Polars,使用百度pip源。

# 安装polars
pip install polars -i https://mirror.baidu.com/pypi/simple/

安装成功后,开始测试,比较Pandas和Polars处理数据的情况。

使用某网站注册用户的用户名数据进行分析,包含约2600万个用户名的CSV文件。

import pandas as pd

df = pd.read_csv(‘users.csv’)
print(df)

数据情况如下。

此外还使用了一个自己创建的CSV文件,用以数据整合测试。

import pandas as pd

df = pd.read_csv(‘fake_user.csv’)
print(df)

得到结果如下。

首先比较一下两个库的排序算法耗时。

import timeit
import pandas as pd

start = timeit.default_timer()

df = pd.read_csv(‘users.csv’)
df.sort_values(‘n’, ascending=False)
stop = timeit.default_timer()

print(‘Time: ‘, stop – start)

————————-
Time:  27.555776743218303

可以看到使用Pandas对数据进行排序,花费了大约28s。

import timeit
import polars as pl

start = timeit.default_timer()

df = pl.read_csv(‘users.csv’)
df.sort(by_column=‘n’, reverse=True)
stop = timeit.default_timer()

print(‘Time: ‘, stop – start)

———————–
Time:  9.924110282212496

Polars只花费了约10s,这意味着Polars比Pandas快了2.7倍。

下面,我们来试试数据整合的效果,纵向连接。

import timeit
import pandas as pd

start = timeit.default_timer()

df_users = pd.read_csv(‘users.csv’)
df_fake = pd.read_csv(‘fake_user.csv’)
df_users.append(df_fake, ignore_index=True)
stop = timeit.default_timer()

print(‘Time: ‘, stop – start)

————————
Time:  15.556222308427095

使用Pandas耗时15s。

import timeit
import polars as pl

start = timeit.default_timer()

df_users = pl.read_csv(‘users.csv’)
df_fake = pl.read_csv(‘fake_user.csv’)
df_users.vstack(df_fake)
stop = timeit.default_timer()

print(‘Time: ‘, stop – start)

———————–
Time:  3.475433263927698

Polars居然最使用了约3.5s,这里Polars比Pandas快了4.5倍。

通过上面的比较,Polars在处理速度上表现得相当不错。

可以是大家在未来处理数据时,另一种选择~

当然,Pandas目前历时12年,已经形成了很成熟的生态,支持很多其它的数据分析库。

Polars则是一个较新的库,不足的地方还有很多。

如果你的数据集对于Pandas来说太大,对于Spark来说太小,那么Polars便是你可以考虑的一个选择。

本文转载自公众号【法纳斯特】

我们的文章到此就结束啦,如果你喜欢今天的 Python 教程,请持续关注Python实用宝典。

有任何问题,可以在公众号后台回复:加群,回答相应验证信息,进入互助群询问。

原创不易,希望你能在下面点个赞和在看支持我继续创作,谢谢!

给作者打赏,选择打赏金额
¥1¥5¥10¥20¥50¥100¥200 自定义

​Python实用宝典 ( pythondict.com )
不只是一个宝典
欢迎关注公众号:Python实用宝典

Python Pyecharts 可视化轻松展示交通拥堵情况

就在今天,我感受到了来自堵车的深深恶意。没有错!我今天被堵在路上近乎3个小时,美好的约会就这样化为泡影了。

我倒还真想看看这路到底能有多堵。于是,我爬取了各城市的拥堵数据,并将它们可视化:

特别说明:由于数据具有实时性,画图时已经过了高峰期,于是图上一片绿油油也并不奇怪。有感兴趣的客官,您接着往下看,待我给您慢慢分解。(ps.涉及到爬虫pyechartsflask等)

一、爬取拥堵指数

某度智慧交通提供了各个城市的拥堵指数的数据,我们只需要通过几行代码便可轻松抓取:

# 获取各城市的拥堵指数  
url = 'https://jiaotong.baidu.com/trafficindex/city/list' # 接口api  
res = requests.get(url)  
data = res.json()  

其中,url为获取数据的接口地址,通过简单的抓包分析便能知道。而data为返回后的数据,它包括很多字段,但是我们只需要提取其中的城市名拥堵指数即可:

# 提取数据  
citys = [i['cityname'for i in data['data']['list']] # 提取城市  
indexs = [float(i['index']) for i in data['data']['list']] # 提取对应的指数  

有了数据,接下来我们就可以将其可视化展示出来。

二、数据可视化

利用可视化神器pyecharts库绘制地图,并将城市以及对应的拥堵指数表示出来。其安装如下:

pip install pyecharts  

部分版本需要再安装额外的地图库,方法如下:

pip install echarts-countries-pypkg  
pip install echarts-cities-pypkg  
pip install echarts-china-provinces-pypkg   
pip install echarts-china-cities-pypkg  

首先定义地图:

geo = Geo()  
geo.add_schema(maptype = 'china'# 加入中国地图  

添加数据并进行相关设置:

geo.add('各城市拥堵指数', zip(citys,indexs), type_ = 'effectScatter'# 设置地图类型及数据  
geo.set_series_opts(label_opts = opts.LabelOpts(is_show = False))  #设置是否显示标签  

根据拥堵指数的大小进行分类,分别为畅通、缓行、拥堵、严重拥堵:

geo.set_global_opts(visualmap_opts = opts.VisualMapOpts(  
                    #max_ = 2.5, # 用于连续表示  
                    is_piecewise = True# 是否分段  
                    pieces = [{'min':1.0,'max':1.5,'label':'畅通','color':'#16CE95'},  
                              {'min':1.5,'max':1.8,'label':'缓行','color':'#F79D06'},  
                              {'min':1.8,'max':2.0,'label':'拥堵','color':'#D80304'},  
                              {'min':2.0,'max':2.5,'label':'严重拥堵','color':'#8F0921'}])) # 设置图例显示  

最后将地图保存在本地:

geo.render(path='各城市拥堵指数.html')  

到这里,我们就得到了文章一开始看到的那张图~

然而,由于拥堵数据是实时变化的,如果我每次都要去运行一次代码岂不是很麻烦?

很显然,机智的社会主义青年是不会这么做的,您接着往下看。

三、搭建展示网站

为了更加方便地将各城市拥堵情况展示出来,我决定搭建一个用于展示的网站。方法可以是各式各样的,在这里我选择了利用flask框架,简单快捷,完整代码回复堵车获得:

代码中,get_geo()为获取地图的函数,返回了pyecharts绘制的地图。在当前目录下创建templates文件夹,并创建模块文件geo.html,如下:

<!DOCTYPE html>  
<html>  
  
<head>  
    <meta charset="utf-8">  
    <title>各城市交通拥堵指数</title>  
</head>  
  
<body>  
  {{mygeo|safe}}  
</body>  
  
</html>  

至此,访问网站地址即可看到绘制的拥堵情况地图~

本文转自快学Python.

我们的文章到此就结束啦,如果你喜欢今天的 Python 教程,请持续关注Python实用宝典。

有任何问题,可以在公众号后台回复:加群,回答相应验证信息,进入互助群询问。

原创不易,希望你能在下面点个赞和在看支持我继续创作,谢谢!

给作者打赏,选择打赏金额
¥1¥5¥10¥20¥50¥100¥200 自定义

​Python实用宝典 ( pythondict.com )
不只是一个宝典
欢迎关注公众号:Python实用宝典

如何使用Python预测机票价格

印度的机票价格基于供需关系浮动,很少受到监管机构的限制。因此它通常被认为是不可预测的,而动态定价机制更增添了人们的困惑。

我们的目的是建立一个机器学习模型,根据历史数据预测未来航班的价格,这些航班价格可以给客户或航空公司服务提供商作为参考价格。

1.准备

开始之前,你要确保Python和pip已经成功安装在电脑上,如果没有,请访问这篇文章:超详细Python安装指南 进行安装。

如果你用Python的目的是数据分析,可以直接安装Anaconda:Python数据分析与挖掘好帮手—Anaconda,它内置了Python和pip.

此外,推荐大家用VSCode编辑器来编写小型Python项目:Python 编程的最好搭档—VSCode 详细指南

Windows环境下打开Cmd(开始—运行—CMD),苹果系统环境下请打开Terminal(command+空格输入Terminal),输入命令安装依赖:

pip install pandas
pip install numpy
pip install matplotlib
pip install seaborn
pip install scikit-learn

2.导入相关数据集

本文的数据集是 Data_Train.xlsx,首先看看训练集的格式:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')


flights = pd.read_excel('./Data_Train.xlsx')
flights.head()

可见训练集中的字段有航空公司(Airline)、日期(Date_of_Journey)、始发站(Source)、终点站(Destination)、路线(Route)、起飞时间(Dep_Time)、抵达时间(Arrival_Time)、历经时长(Duration)、总计停留站点个数(Total_Stops)、额外信息(Additional_Info),最后是机票价格(Price)。

与其相对的测试集,除了缺少价格字段之外,与训练集的其他所有字段均一致。

访问:https://pythondict.com/download/predict-ticket/

或在Python实用宝典后台回复:预测机票,下载完整数据源和代码。

3.探索性数据分析

3.1 清理缺失数据

看看所有字段的基本信息:

flights.info()

其他的非零值数量均为10683,只有路线和停靠站点数是10682,说明这两个字段缺少了一个值。

谨慎起见,我们删掉缺少数据的行:

# clearing the missing data
flights.dropna(inplace=True)
flights.info()

现在非零值达到一致数量,数据清理完毕。

3.2 航班公司分布特征

接下来看看航空公司的分布特征:

sns.countplot('Airline', data=flights)
plt.xticks(rotation=90)
plt.show()

前三名的航空公司分别是 IndiGo, Air India, JetAirways.

其中可能存在廉价航空公司。

3.3 再来看看始发地的分布

sns.countplot('Source',data=flights)
plt.xticks(rotation=90)
plt.show()

某些地区可能是冷门地区,存在冷门机票的可能性比较大。

3.4 停靠站点的数量分布

sns.countplot('Total_Stops',data=flights)
plt.xticks(rotation=90)
plt.show()

看来大部分航班在飞行途中只停靠一次或无停靠。

会不会某些停靠多的航班比较便宜?

3.5 有多少数据含有额外信息

plot=plt.figure()
sns.countplot('Additional_Info',data=flights)
plt.xticks(rotation=90)

大部分航班信息中都没有包含额外信息,除了部分航班信息有:不包含飞机餐、不包含免费托运。

这个信息挺重要的,是否不包含这两项服务的飞机机票比较便宜?

3.6 时间维度分析

首先转换时间格式:

flights['Date_of_Journey'] = pd.to_datetime(flights['Date_of_Journey'])
flights['Dep_Time'] = pd.to_datetime(flights['Dep_Time'],format='%H:%M:%S').dt.time

接下来,研究一下出发时间和价格的关系:

flights['weekday'] = flights[['Date_of_Journey']].apply(lambda x:x.dt.day_name())
sns.barplot('weekday','Price',data=flights)
plt.show()

大体上价格没有差别,说明这个特征是无效的。

那么月份和机票价格的关系呢?

flights["month"] = flights['Date_of_Journey'].map(lambda x: x.month_name())
sns.barplot('month','Price',data=flights)
plt.show()

没想到4月的机票价格均价只是其他月份的一半,看来4月份是印度的出行淡季吧。

起飞时间和价格的关系

flights['Dep_Time'] = flights['Dep_Time'].apply(lambda x:x.hour)
flights['Dep_Time'] = pd.to_numeric(flights['Dep_Time'])
sns.barplot('Dep_Time','Price',data=flights)
plot.show()

可以看到,红眼航班(半夜及早上)的机票比较便宜,这是符合我们的认知的。

3.7 清除无效特征

把那些和价格没有关联关系的字段直接去除掉:

flights.drop(['Route','Arrival_Time','Date_of_Journey'],axis=1,inplace=True)
flights.head()

4.模型训练

接下来,我们可以准备使用模型来预测机票价格了,不过,还需要对数据进行预处理和特征缩放。

4.1 数据预处理

将字符串变量使用数字替代:

from sklearn.preprocessing import LabelEncoder
var_mod = ['Airline','Source','Destination','Additional_Info','Total_Stops','weekday','month','Dep_Time']
le = LabelEncoder()
for i in var_mod:
    flights[i] = le.fit_transform(flights[i])
flights.head()

对每列数据进行特征缩放,提取自变量(x)和因变量(y):

flights.corr()
def outlier(df):
    for i in df.describe().columns:
        Q1=df.describe().at['25%',i]
        Q3=df.describe().at['75%',i]
        IQR= Q3-Q1
        LE=Q1-1.5*IQR
        UE=Q3+1.5*IQR
        df[i]=df[i].mask(df[i]<LE,LE)
        df[i]=df[i].mask(df[i]>UE,UE)
    return df
flights = outlier(flights)
x = flights.drop('Price',axis=1)
y = flights['Price']

划分测试集和训练集:

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=101)

4.2 模型训练及测试

使用随机森林进行模型训练:

from sklearn.ensemble import RandomForestRegressor
rfr=RandomForestRegressor(n_estimators=100)
rfr.fit(x_train,y_train)

在随机森林中,我们有一种根据数据的相关性来确定特征重要性的方法:

features=x.columns
importances = rfr.feature_importances_
indices = np.argsort(importances)
plt.figure(1)
plt.title('Feature Importances')
plt.barh(range(len(indices)), importances[indices], color='b', align='center')
plt.yticks(range(len(indices)), features[indices])
plt.xlabel('Relative Importance')

可以看到,Duration(飞行时长)是影响最大的因子。

对划分的测试集进行预测,得到结果:

predictions=rfr.predict(x_test)
plt.scatter(y_test,predictions)
plt.show()

这样看不是很直观,接下来我们要数字化地评价这个模型。

4.3 模型评价

sklearn 提供了非常方便的函数来评价模型,那就是 metrics :

from sklearn import metrics
print('MAE:', metrics.mean_absolute_error(y_test, predictions))
print('MSE:', metrics.mean_squared_error(y_test, predictions))
print('RMSE:', np.sqrt(metrics.mean_squared_error(y_test, predictions)))
print('r2_score:', (metrics.r2_score(y_test, predictions)))
MAE: 1453.9350628905618
MSE: 4506308.3645551
RMSE: 2122.806718605135
r2_score: 0.7532074710409375

这4个值中你可以只关注R2_score,r2越接近1说明模型效果越好,这个模型的分数是0.75,算是很不错的模型了。

看看其残差直方图是否符合正态分布:

sns.distplot((y_test-predictions),bins=50)
plt.show()

不错,多数预测结果和真实值都在-1000到1000的范围内,算是可以接受的结果。其残差直方图也基本符合正态分布,说明模型是有效果的。

部分译自 https://www.kaggle.com/harikrishna9/how-to-predict-flight-ticket-price/notebook,有较多的增删。

我们的文章到此就结束啦,如果你喜欢今天的 Python 教程,请持续关注Python实用宝典。

有任何问题,可以在公众号后台回复:加群,回答相应验证信息,进入互助群询问。

原创不易,希望你能在下面点个赞和在看支持我继续创作,谢谢!

给作者打赏,选择打赏金额
¥1¥5¥10¥20¥50¥100¥200 自定义

​Python实用宝典 ( pythondict.com )
不只是一个宝典
欢迎关注公众号:Python实用宝典

Python 利用股市找到那些存在问题的上市企业

当某个股票保持下跌的时候,你可以肯定该股票一定有什么地方不对,要么是它的市场,要么是它的公司。——利弗莫尔

因此结合自己最近的研究,提出了这个比较新颖的想法:利用股票市场找到已有或潜在食品安全问题的企业。

为什么会有这个想法呢?我们知道食品安全事故发生的时间点和媒体报道的时间点之间实际上是有一个间隔的,然而相关利益人士、内部人士却能提前知道事件的发生。

比如说,2012年11月19日,酒鬼酒被国家质检总局爆出塑化剂超标247%,但是在19号之前,酒鬼酒的收盘价却神奇地从2012年11月2日的55元回落到2012年11月16号的46元。然而国家质检总局的消息是在2012年11月19号才发布的。

而且这样的下跌,明显违背该股票当时上涨的趋势,在10月底时,它的股价已突破周K的压力线,却很不自然地下跌了?而且还违背了许多技术指标,如1号、14号、16号,KDJ和BOLL指标明显提示上涨,在股市这种情绪化的市场中,有人却不为所动,仍然售出大量股票导致其不正常下跌。

2012年11月19号消息公布时紧急停牌,复牌后有三个跌停,股价暴跌48%,但我相信某些人已经成功躲开了这场股灾。

找到已有或潜在食品安全问题的企业的重点在于两个方面:

1.该股票在该板块或者该股强走势的情况下却发生连续多日的下跌

2.不正常的跌停板

如果我们在消息公布/(不公布)前能提前捕捉到这个异常信息,我们就能提前捕捉到某个食品的安全问题,而且也能规避投资风险。当然,出现这种异常的可能性非常多,这种异常只能作为参考。

下面让我们尝试用Python来找到第一种趋势的企业,当然我最后不会公布结果,大家感兴趣可以自己试试:

首先利用tushare找到和食品安全相关的上市企业:

import tushare as ts
def food_codes():
    data = ts.get_industry_classified()
    print data[data.c_name.isin(['食品行业','农药化肥','酿酒行业'])]

得到结果:

均线是我们获得该股票趋势的基础,下面我们编写均线函数:

def get_ma(code,start='',end=datetime.date.today().strftime("%Y-%m-%d")):
    data = ts.get_k_data(code)
    data = data.sort_index(ascending=False)
    data['ma2'] = data['close'].rolling(2).mean().shift(-1)
    data['ma5'] = data['close'].rolling(5).mean().shift(-4)
    data['ma10'] = data['close'].rolling(10).mean().shift(-9)
    data['ma20'] = data['close'].rolling(20).mean().shift(-19)
    data['ma60'] = data['close'].rolling(60).mean().shift(-59)
    data['ma240'] = data['close'].rolling(240).mean().shift(-239)
    data['date'] = pd.to_datetime(data['date'])
    if start == '':
        return data
    start = pd.to_datetime(start)
    end = pd.to_datetime(end)
    if data['date'][len(data) - 1] < start:
        return 0
    while data.loc[data.date == start].empty:
        start = start + dateutil.relativedelta.relativedelta(days=1)
    while data.loc[data.date == end].empty:
        end = end - dateutil.relativedelta.relativedelta(days=1)
    return data.loc[(data.date >= start) & (data.date <= end)]

我们只需要确定两点:

1. 以20个交易日为窗口期,其内最高价和最低价的差距大于其最高价的15%。

2. 该股票处于涨势。

def analyzeOne(code):
    rng = pd.date_range('2018-1-1', datetime.date.today().strftime("%Y-%m-%d"), freq='D')
    # 获得日期
    flag = 0
    for i in range(20,len(rng)):
        data = get_ma(code,rng[i-20],rng[i])
        count = 0
        data = data.sort_index(ascending=True)
        # 升序
        max = data['close'][data['close'].argmax()]
        min = data['close'][data['close'].argmin()]
        for j in range(len(data)):
            if data.iloc[j]['ma5'] >= data.iloc[j]['ma20'] and data.iloc[j]['ma20'] >= data.iloc[j]['ma60'] and \
                    data.iloc[j]['ma60'] >= data.iloc[j]['ma240']:
                count = count + 1
                if count >= 5 and (max-min) - max*0.15 > 0:
                    print 'Code: ' + str(code) +', Problem:' + str(rng[i-20])+ ' ' + str(rng[i])
            else:
                count = 0

如果存在这样的趋势,就输出这个趋势的窗口时间段,当然这只是一个粗略的模型,结果中可能会有很多意外的情况,但是我觉得足够启发大家了。

想要应用于所有食品相关股票:

def find_down():
    for i in food_codes().code:
        analyzeOne(i)

通过这样的操作,我找到了两支类似的股票:

但这些股票都没有新闻报道发生了问题,最后也让我怀疑自己的模型的准确性,到底是市场正常波动,还是这个企业存在问题。

因此,这个模型存在许多的优化空间,如何通过股市最准确地找出那些出问题的企业,还有待进一步研究,本实验只是提出一个初步的模型和一些想法。

个人愚见,欢迎讨论。

我们的文章到此就结束啦,如果你喜欢今天的Python 实战教程,请持续关注Python实用宝典。

有任何问题,可以在公众号后台回复:加群,回答相应红字验证信息,进入互助群询问。

原创不易,希望你能在下面点个赞和在看支持我继续创作,谢谢!

点击下方阅读原文可获得更好的阅读体验

Python实用宝典 (pythondict.com)
不只是一个宝典
欢迎关注公众号:Python实用宝典

Python 2015-2020年美国警察致命枪击案EDA分析

2014年在密苏里州一名叫做弗格森(Ferguson)的警察杀害了迈克尔·布朗(Michael Brown)后,美国黑人开始了一场抗议警察暴力对待黑人的运动—Black Lives Matter(黑人的命也是命,简称BLM)。

2020年,在明尼阿波利斯警察Derek Chauvin杀害乔治·弗洛伊德(George Floyd)之后,BLM运动再次成为头条新闻,引起国际社会的进一步关注。

自2015年1月1日起,《华盛顿邮报》一直在整理一个数据库,其中记录了值班警员在美国发生的每起致命枪击事件。这个数据库里包含了死者的种族,年龄和性别,该人是否有武器,以及受害人是否正在遭受精神健康危机。

此外,还有四个其他数据集。有关贫困率,高中毕业率,家庭收入中位数和种族人口统计数据的美国人口普查数据。

下面就让我们来使用这些数据集来进行数据分析。

1.准备

开始之前,你要确保Python和pip已经成功安装在电脑上,如果没有,请访问这篇文章:超详细Python安装指南 进行安装。如果你用Python的目的是数据分析,可以直接安装Anaconda:Python数据分析与挖掘好帮手—Anaconda,它内置了Python和pip.

此外,推荐大家用VSCode编辑器,因为它可以在编辑器下方的终端运行命令安装依赖模块:Python 编程的最好搭档—VSCode 详细指南。

本文具备流程性,建议使用 VSCode 的 Jupiter Notebook 扩展,新建一个名为 test.ipynb 的文件,跟着教程一步步走下去。

Windows环境下打开Cmd(开始—运行—CMD),苹果系统环境下请打开Terminal(command+空格输入Terminal),准备开始输入命令安装依赖。

所需依赖:

pip install numpy
pip install pandas
pip install plotly
pip install seaborn

本文译自:https://www.kaggle.com/edoardo10/fatal-police-shooting-eda-plotly-seaborn/data,如需数据请在公众号后台回复:警察枪击EDA

2.代码与分析

首先,引入我们分析所需要使用的模块:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from datetime import datetime
import plotly.express as px
import plotly.graph_objects as go
import warnings
import plotly.offline as pyo
pyo.init_notebook_mode()
warnings.filterwarnings('ignore')
pd.set_option('display.max_columns', 500)
sns.set_style('white')
%matplotlib inline

打开需要分析的数据集:

df = pd.read_csv('./PoliceKillingsUS.csv', encoding='cp1252')
df.head()

时间特征

从这6年的月度数据来看,我们可以看到,在2015年上半年、2018年初和2020年第一季度,我们达到了每月超过100起致命事故的高峰。从月度来看,这种现象不具备明显的季节性。

df['date'] = df['date'].apply(lambda x: pd.to_datetime(x))
df['date'].groupby(df.date.dt.to_period('M')).count().plot(kind='line')

看看警察枪击案的事故是否具有周末特征:

count = df['date'].apply(lambda x: 'Weekday' if x.dayofweek < 5 else 'Weekend').value_counts(normalize=True)
f, ax = plt.subplots(1,1)
sns.barplot(x=count.index, y=count.values, ax=ax, palette='twilight')

显然,我们没有证据表明周末会发生更多的案件。

不过,如果细化到星期里的每一天,我们会发现周中发生案件的概率较高:

count = df['date'].apply(lambda x: x.dayofweek).value_counts(normalize=True).sort_index()
count.index = ['Mon','Tue','Wed','Thu','Fri','Sat','Sun']
f, ax = plt.subplots(1,1)
sns.barplot(x=count.index, y=count.values, ax=ax, palette='twilight')
ax.set_title('Cases (%) for each day of the week');

接下来看看以下4个特征的分布:

signs_of_mental_illness:是否精神不稳定
threat_level:威胁等级
body_camera:警察是否带了随身摄像头
manner_of_death:死亡方式
count_1 = df['signs_of_mental_illness'].value_counts(normalize=True)
count_2 = df['threat_level'].value_counts(normalize=True)
count_3 = df['body_camera'].value_counts(normalize=True)
count_4 = df['manner_of_death'].value_counts(normalize=True)
fig, axes = plt.subplots(2, 2, figsize=(8, 8), sharey=True)
sns.barplot(x=count_1.index, y=count_1.values, palette="rocket", ax=axes[0,0])
axes[0,0].set_title('Signs of mental illness (%)')
sns.barplot(x=count_2.index, y=count_2.values, palette="viridis", ax=axes[0,1])
axes[0,1].set_title('Threat level (%)')
sns.barplot(x=count_3.index, y=count_3.values, palette="nipy_spectral", ax=axes[1,0])
axes[1,0].set_title('Body camera (%)')
sns.barplot(x=count_4.index, y=count_4.values, palette="gist_heat", ax=axes[1,1])
axes[1,1].set_title('Manner of death (%)');

我们可以看到,只有20%的案例受害者有精神不稳定的迹象;

只有10%的警察有随身摄像头;

70%的情况被宣布为危险状况;

死亡方式似乎不是一个有趣的变量,因为大多数案件都是“枪毙”;

美国的警察是否具有种族主义倾向?

count = df.race.value_counts(normalize=True)
count.index = ['White', 'Black', 'Hispanic', 'Asian', 'Native American', 'Other']

f, ax = plt.subplots(1,1, figsize=(8,6))
sns.barplot(y=count.index, x=count.values, palette='Reds_r')
ax.set_title('Total cases for each race (%)');

从上图我们知道,大部分致命的枪击事件中,涉及最多的是白人,其次是黑人和西班牙裔。

但这个图表并没有考虑人种比例。参考2019年美国的种族比例,我们可以看到,美国黑人受害者的比例更高:

数据来源:https://data.census.gov/cedsci/table?q=Hispanic%20or%20Latino&tid=ACSDP1Y2019.DP05&hidePreview=false

share_race_usa_2019 = pd.Series([60.0, 12.4, 0.9, 5.6, 18.4, 2.7], index=['White','Black','Native American','Asian','Hispanic','Other'])

count_races = count / share_race_usa_2019
count_races = count_races.sort_values(ascending=False)
f, ax = plt.subplots(1,1, figsize=(8,6))
sns.barplot(y=count_races.index, x=count_races.values, palette='Greens_r')
ax.set_title('Total cases for each race on total USA race percentage rate');

受害者的年龄

sns.set_style('whitegrid')
fig, axes = plt.subplots(1, 1, figsize=(10, 8))
axes.xaxis.set_ticks(np.arange(0,100,10))

sns.kdeplot(df[df.race == 'N'].age, ax=axes, shade=True, color='#7FFFD4')
sns.kdeplot(df[df.race == 'O'].age, ax=axes, shade=True, color='#40E0D0')
sns.kdeplot(df[df.race == 'B'].age, ax=axes, shade=True, color='#00CED1')
sns.kdeplot(df[df.race == 'H'].age, ax=axes, shade=True, color='#6495ED')
sns.kdeplot(df[df.race == 'A'].age, ax=axes, shade=True, color='#4682B4')
sns.kdeplot(df[df.race == 'W'].age, ax=axes, shade=True, color='#008B8B')


legend = axes.legend_
legend.set_title("Race")
for t, l in zip(legend.texts,("Native", "Other", 'Black', 'Hispanic', 'Asian', 'White')):
    t.set_text(l)

由这些叠加的密度图可以看出:

对于亚裔和白人来说,大多数案件的受害者年龄都在30岁左右。

对于其他和印第安人来说,在大多数案件中,受害者大约28岁。

对于西班牙裔和黑人来说,大多数案件的受害者年龄都在25岁左右。

所以我们可以说,西班牙裔和黑人的年轻人,是被警察开枪射击的高危群体。

受害者性别比例

按常理,这种暴力事件的受害者一般都为男性,看看是不是这样:

fig = px.pie(values = df.gender.value_counts(normalize=True).values, names=df.gender.value_counts(normalize=True).index, title='Total cases gender (%)')
fig.update(layout=dict(title=dict(x=0.5),autosize=False, width=400, height=400))
fig.show()

果然如此,超过95%的受害者都为男性。

简单的EDA分析就是这些,作者还分享了许多深层次的分析,不过并没有将数据分享出来,这里就不展示了。

我们的文章到此就结束啦,如果你喜欢今天的 Python 教程,请持续关注Python实用宝典。

有任何问题,可以在公众号后台回复:加群,回答相应验证信息,进入互助群询问。

原创不易,希望你能在下面点个赞和在看支持我继续创作,谢谢!

给作者打赏,选择打赏金额
¥1¥5¥10¥20¥50¥100¥200 自定义

​Python实用宝典 ( pythondict.com )
不只是一个宝典
欢迎关注公众号:Python实用宝典

Python 预测银行信用卡客户是否流失

1.题目背景

越来越多的客户不再使用信用卡服务,银行的经理对此感到不安。如果有人能为他们预测哪些客户即将流失,他们将不胜感激,因为这样他们可以主动向客户提供更好的服务,并挽回这些即将流失的客户。

2.数据集

该数据集由10,000个客户组成,其中包含了他们的年龄,工资,婚姻状况,信用卡限额,信用卡类别等。

不过,这里面只有16%的客户是流失的,因此拿来预测客户是否会流失有点难度。

在Python实用宝典后台回复 预测客户流失 下载这份数据和源代码。

译自kaggle并对原文进行了修改和补充,感谢原作者:

https://www.kaggle.com/thomaskonstantin/bank-churn-data-exploration-and-churn-prediction/

3.代码与分析

开始之前,你要确保Python和pip已经成功安装在电脑上,如果没有,请访问这篇文章:超详细Python安装指南 进行安装。如果你用Python的目的是数据分析,可以直接安装Anaconda:Python数据分析与挖掘好帮手—Anaconda,它内置了Python和pip.

此外,推荐大家用VSCode编辑器,因为它可以在编辑器下方的终端运行命令安装依赖模块:Python 编程的最好搭档—VSCode 详细指南。

本文具备流程性,建议使用 VSCode 的 Jupiter Notebook 扩展,新建一个名为 test.ipynb 的文件,跟着教程一步步走下去。

Windows环境下打开Cmd(开始—运行—CMD),苹果系统环境下请打开Terminal(command+空格输入Terminal),准备开始输入命令安装依赖。

所需依赖:

pip install numpy
pip install pandas
pip install plotly
pip install scikit-learn
pip install scikit-plot

# 这个需要conda
conda install -c conda-forge imbalanced-learn

3.1 导入需要的模块

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as ex
import plotly.graph_objs as go
import plotly.figure_factory as ff
from plotly.subplots import make_subplots
import plotly.offline as pyo
pyo.init_notebook_mode()
sns.set_style('darkgrid')
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split,cross_val_score
from sklearn.ensemble import RandomForestClassifier,AdaBoostClassifier
from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import f1_score as f1
from sklearn.metrics import confusion_matrix
import scikitplot as skplt

plt.rc('figure',figsize=(18,9))
%pip install imbalanced-learn
from imblearn.over_sampling import SMOTE

遇到任何 No module named “XXX” 都可以尝试pip install一下。

如果pip install没解决,可以谷歌/百度一下,看看别人是怎么解决的。

3.2 加载数据

c_data = pd.read_csv('./BankChurners.csv')
c_data = c_data[c_data.columns[:-2]]
c_data.head(3)

这里去掉了最后两列的朴素贝叶斯分类结果。

显示前三行数据, 可以看到所有的字段:

3.3 探索性数据分析

下面看看这20+列数据中,哪一些是对我们有用的。

首先,我想知道数据集中的客户年龄分布:

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Customer_Age'],name='Age Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Customer_Age'],name='Age Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of Customer Ages")
fig.show()

可以看到,客户的年龄分布大致遵循正态分布,因此使用可以在正态假设下进一步使用年龄特征。

同样滴,我想知道性别分布如何:

ex.pie(c_data,names='Gender',title='Propotion Of Customer Genders')

可见,在我们的数据集中,女性的样本比男性更多,但是差异的百分比不是那么显著,所以我们可以说性别是均匀分布的。

每个客户的家庭人数的分布怎么样?

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Dependent_count'],name='Dependent count Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Dependent_count'],name='Dependent count Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of Dependent counts (close family size)")
fig.show()

可见,它也是大致符合正态分布的,偏右一点。

客户的受教育水平如何?

ex.pie(c_data,names='Education_Level',title='Propotion Of Education Levels')

假设大多数教育程度不明(Unknown)的顾客都没有接受过任何教育。我们可以指出,超过70%的顾客都受过正规教育,其中约35%的人受教育程度达到硕士以上水平,45%的人达到本科以上水准。

他们的婚姻状态如何?

ex.pie(c_data,names='Marital_Status',title='Propotion Of Different Marriage Statuses')

看来,这家银行几乎一半的客户都是已婚人士,有趣的是,另一半客户几乎都是单身人士,另外只有7%的客户离婚了。

看看收入分布和卡片类型的分布:

ex.pie(c_data,names='Income_Category',title='Propotion Of Different Income Levels')
ex.pie(c_data,names='Card_Category',title='Propotion Of Different Card Categories')

可见大部分人的年收入处于60K美元以下。

在持有的卡片的类型上,蓝卡占了绝大多数。

每月账单数量有没有特征?

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Months_on_book'],name='Months on book Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Months_on_book'],name='Months on book Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of months the customer is part of the bank")
fig.show()

可以看到中间的峰值特别高,显然这个指标不是正态分布的。

每位客户持有的银行业务数量有没有特征呢?

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Total_Relationship_Count'],name='Total no. of products Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Total_Relationship_Count'],name='Total no. of products Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of Total no. of products held by the customer")
fig.show()

基本上都是均匀分布的,显然这个指标对于我们而言也没太大意义。

用户不活跃月份数量有没有好用的特征?

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Months_Inactive_12_mon'],name='number of months inactive Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Months_Inactive_12_mon'],name='number of months inactive Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of the number of months inactive in the last 12 months")
fig.show()

这个似乎有点用处,会不会越不活跃的用户越容易流失呢?

信用卡额度的分布如何?

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Credit_Limit'],name='Credit_Limit Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Credit_Limit'],name='Credit_Limit Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of the Credit Limit")
fig.show()

大部分人的额度都在0到10k之间,这比较正常,暂时看不出和流失有什么关系。

客户总交易额的分布怎么样?

fig = make_subplots(rows=2, cols=1)

tr1=go.Box(x=c_data['Total_Trans_Amt'],name='Total_Trans_Amt Box Plot',boxmean=True)
tr2=go.Histogram(x=c_data['Total_Trans_Amt'],name='Total_Trans_Amt Histogram')

fig.add_trace(tr1,row=1,col=1)
fig.add_trace(tr2,row=2,col=1)

fig.update_layout(height=700, width=1200, title_text="Distribution of the Total Transaction Amount (Last 12 months)")
fig.show()

这个有点意思,总交易额的分布体现出“多组”分布,如果我们根据这个指标将客户聚类为不同的组别,看他们之间的相似性,并作出不同的画线,也许对我们最终的流失分析有一定的意义。

接下来,最重要的流失用户分布

ex.pie(c_data,names='Attrition_Flag',title='Proportion of churn vs not churn customers')

我们可以看到,只有16%的数据样本代表流失客户,在接下来的步骤中,我将使用SMOTE对流失样本进行采样,使其与常规客户的样本大小匹配,以便给后面选择的模型一个更好的机会来捕捉小细节。

3.4 数据预处理

使用SMOTE模型前,需要根据不同的特征对数据进行One Hot编码:

c_data.Attrition_Flag = c_data.Attrition_Flag.replace({'Attrited Customer':1,'Existing Customer':0})
c_data.Gender = c_data.Gender.replace({'F':1,'M':0})
c_data = pd.concat([c_data,pd.get_dummies(c_data['Education_Level']).drop(columns=['Unknown'])],axis=1)
c_data = pd.concat([c_data,pd.get_dummies(c_data['Income_Category']).drop(columns=['Unknown'])],axis=1)
c_data = pd.concat([c_data,pd.get_dummies(c_data['Marital_Status']).drop(columns=['Unknown'])],axis=1)
c_data = pd.concat([c_data,pd.get_dummies(c_data['Card_Category']).drop(columns=['Platinum'])],axis=1)
c_data.drop(columns = ['Education_Level','Income_Category','Marital_Status','Card_Category','CLIENTNUM'],inplace=True)

显示热力图:

sns.heatmap(c_data.corr('pearson'),annot=True)

3.5 SMOTE模型采样

SMOTE模型经常用于解决数据不平衡的问题,它通过添加生成的少数类样本改变不平衡数据集的数据分布,是改善不平衡数据分类模型性能的流行方法之一。

oversample = SMOTE()
X, y = oversample.fit_resample(c_data[c_data.columns[1:]], c_data[c_data.columns[0]])
usampled_df = X.assign(Churn = y)
ohe_data =usampled_df[usampled_df.columns[15:-1]].copy()
usampled_df = usampled_df.drop(columns=usampled_df.columns[15:-1])
sns.heatmap(usampled_df.corr('pearson'),annot=True)

3.6 主成分分析

我们将使用主成分分析来降低单次编码分类变量的维数,从而降低方差。同时使用几个主成分而不是几十个单次编码特征将帮助我构建一个更好的模型。

N_COMPONENTS = 4

pca_model = PCA(n_components = N_COMPONENTS )

pc_matrix = pca_model.fit_transform(ohe_data)

evr = pca_model.explained_variance_ratio_
cumsum_evr = np.cumsum(evr)

ax = sns.lineplot(x=np.arange(0,len(cumsum_evr)),y=cumsum_evr,label='Explained Variance Ratio')
ax.set_title('Explained Variance Ratio Using {} Components'.format(N_COMPONENTS))
ax = sns.lineplot(x=np.arange(0,len(cumsum_evr)),y=evr,label='Explained Variance Of Component X')
ax.set_xticks([i for i in range(0,len(cumsum_evr))])
ax.set_xlabel('Component number #')
ax.set_ylabel('Explained Variance')
plt.show()
usampled_df_with_pcs = pd.concat([usampled_df,pd.DataFrame(pc_matrix,columns=['PC-{}'.format(i) for i in range(0,N_COMPONENTS)])],axis=1)
usampled_df_with_pcs

特征变得越来越明显:

sns.heatmap(usampled_df_with_pcs.corr('pearson'),annot=True)

4.模型选择及测试

选择出以下特征划分训练集并进行训练:

X_features = ['Total_Trans_Ct','PC-3','PC-1','PC-0','PC-2','Total_Ct_Chng_Q4_Q1','Total_Relationship_Count']

X = usampled_df_with_pcs[X_features]
y = usampled_df_with_pcs['Churn']

train_x,test_x,train_y,test_y = train_test_split(X,y,random_state=42)

4.1 交叉验证

分别看看随机森林、AdaBoost和SVM模型三种模型的表现如何:

rf_pipe = Pipeline(steps =[ ('scale',StandardScaler()), ("RF",RandomForestClassifier(random_state=42)) ])
ada_pipe = Pipeline(steps =[ ('scale',StandardScaler()), ("RF",AdaBoostClassifier(random_state=42,learning_rate=0.7)) ])
svm_pipe = Pipeline(steps =[ ('scale',StandardScaler()), ("RF",SVC(random_state=42,kernel='rbf')) ])


f1_cross_val_scores = cross_val_score(rf_pipe,train_x,train_y,cv=5,scoring='f1')
ada_f1_cross_val_scores=cross_val_score(ada_pipe,train_x,train_y,cv=5,scoring='f1')
svm_f1_cross_val_scores=cross_val_score(svm_pipe,train_x,train_y,cv=5,scoring='f1')
plt.subplot(3,1,1)
ax = sns.lineplot(x=range(0,len(f1_cross_val_scores)),y=f1_cross_val_scores)
ax.set_title('Random Forest Cross Val Scores')
ax.set_xticks([i for i in range(0,len(f1_cross_val_scores))])
ax.set_xlabel('Fold Number')
ax.set_ylabel('F1 Score')
plt.show()
plt.subplot(3,1,2)
ax = sns.lineplot(x=range(0,len(ada_f1_cross_val_scores)),y=ada_f1_cross_val_scores)
ax.set_title('Adaboost Cross Val Scores')
ax.set_xticks([i for i in range(0,len(ada_f1_cross_val_scores))])
ax.set_xlabel('Fold Number')
ax.set_ylabel('F1 Score')
plt.show()
plt.subplot(3,1,3)
ax = sns.lineplot(x=range(0,len(svm_f1_cross_val_scores)),y=svm_f1_cross_val_scores)
ax.set_title('SVM Cross Val Scores')
ax.set_xticks([i for i in range(0,len(svm_f1_cross_val_scores))])
ax.set_xlabel('Fold Number')
ax.set_ylabel('F1 Score')
plt.show()

看看三种模型都有什么不同的表现:

看得出来随机森林 F1分数是最高的。

4.2 模型预测

对测试集进行预测,看看三种模型的效果:

rf_pipe.fit(train_x,train_y)
rf_prediction = rf_pipe.predict(test_x)

ada_pipe.fit(train_x,train_y)
ada_prediction = ada_pipe.predict(test_x)

svm_pipe.fit(train_x,train_y)
svm_prediction = svm_pipe.predict(test_x)

print('F1 Score of Random Forest Model On Test Set - {}'.format(f1(rf_prediction,test_y)))
print('F1 Score of AdaBoost Model On Test Set - {}'.format(f1(ada_prediction,test_y)))
print('F1 Score of SVM Model On Test Set - {}'.format(f1(svm_prediction,test_y)))

4.3 对原始数据(采样前)进行模型预测

​接下来对原始数据进行模型预测:

ohe_data =c_data[c_data.columns[16:]].copy()
pc_matrix = pca_model.fit_transform(ohe_data)
original_df_with_pcs = pd.concat([c_data,pd.DataFrame(pc_matrix,columns=['PC-{}'.format(i) for i in range(0,N_COMPONENTS)])],axis=1)

unsampled_data_prediction_RF = rf_pipe.predict(original_df_with_pcs[X_features])
unsampled_data_prediction_ADA = ada_pipe.predict(original_df_with_pcs[X_features])
unsampled_data_prediction_SVM = svm_pipe.predict(original_df_with_pcs[X_features])

效果如下:

F1最高的随机森林模型有0.63分,偏低,这也比较正常,毕竟在这种分布不均的数据集中,查全率是很难做到很高的。

4.4 结果

让我们看看最终在原数据上使用随机森林模型的运行结果:

ax = sns.heatmap(confusion_matrix(unsampled_data_prediction_RF,original_df_with_pcs['Attrition_Flag']),annot=True,cmap='coolwarm',fmt='d')
ax.set_title('Prediction On Original Data With Random Forest Model Confusion Matrix')
ax.set_xticklabels(['Not Churn','Churn'],fontsize=18)
ax.set_yticklabels(['Predicted Not Churn','Predicted Churn'],fontsize=18)

plt.show()

可见,没有流失的客户命中了7709人,未命中791人。

流失客户命中了1130人,未命中497人。

整体而言,是一个比较优秀的模型了。

我们的文章到此就结束啦,如果你喜欢今天的 Python 教程,请持续关注Python实用宝典。

有任何问题,可以在公众号后台回复:加群,回答相应验证信息,进入互助群询问。

原创不易,希望你能在下面点个赞和在看支持我继续创作,谢谢!

给作者打赏,选择打赏金额
¥1¥5¥10¥20¥50¥100¥200 自定义

​Python实用宝典 ( pythondict.com )
不只是一个宝典
欢迎关注公众号:Python实用宝典