Author Archive for 写长城的诗

社会科学的代码和数据工作指南

偶尔在知乎上看到有人推荐了一本小册子:《Code and Data for the Social Sciences:A Practitioner’s Guide》。专门讲非计算机背景的分析研究人员如何归整自己的分析代码和研究数据。看下来还是总结得非常好,很有益于创建高效的工作规范和流程。将其中一些基本的规则摘要如下:

Automate
(A) Automate everything that can be automated.
(B) Write a single script that executes all code from beginning to end.

Version Control
(A) Store code and data under version control.
(B) Run the whole directory before checking it back in.

Directories
(A) Separate directories by function.
(B) Separate files into inputs and outputs.
(C) Make directories portable.

\input CSV
\code R, SQL
\output pic, ppt
\temp
readme.txt

Keys
(A) Store cleaned data in tables with unique, non-missing keys.
(B) Keep data normalized as far into your code pipeline as you can

Abstraction
(A) Abstract to eliminate redundancy.
(B) Abstract to improve clarity.
(C) Otherwise, don’t abstract.

Documentation
(A) Don’t write documentation you will not maintain.
(B) Code should be self-documenting.

Management
(A) Manage tasks with a task management system.
(B) E-mail is not a task management system.

python的决策树和随机森林

Python的决策树和随机森林

决策树模型是一种简单易用的非参数分类器。它不需要对数据有任何的先验假设,计算速度较快,结果容易解释,而且稳健性强,对噪声数据和缺失数据不敏感。下面示范用titanic中的数据集为做决策树分类,目标变量为survive。

第一步:读取数据

In [2]:
%pylab inline
import pandas as pd
df = pd.read_csv('titanic.csv')
df.head()
#df.info()

Populating the interactive namespace from numpy and matplotlib

Out[2]:
survived pclass name sex age sibsp parch ticket fare cabin embarked
0 0 3 Braund, Mr. Owen Harris male 22 1 0 A/5 21171 7.2500 NaN S
1 1 1 Cumings, Mrs. John Bradley (Florence Briggs Th… female 38 1 0 PC 17599 71.2833 C85 C
2 1 3 Heikkinen, Miss. Laina female 26 0 0 STON/O2. 3101282 7.9250 NaN S
3 1 1 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35 1 0 113803 53.1000 C123 S
4 0 3 Allen, Mr. William Henry male 35 0 0 373450 8.0500 NaN S

第二步:数据整理

  • 只取出三个自变量
  • 将将age缺失值进行补全
  • 将pclass变量转为三个哑变量
  • 将sex转为0-1变量
In [3]:
subdf = df[['pclass','sex','age']]
y = df.survived
# sklearn中的Imputer也可以
age = subdf['age'].fillna(value=subdf.age.mean())
# sklearn OneHotEncoder也可以
pclass = pd.get_dummies(subdf['pclass'],prefix='pclass')
sex = (subdf['sex']=='male').astype('int')
X = pd.concat([pclass,age,sex],axis=1)
X.head()
Out[3]:
pclass_1 pclass_2 pclass_3 age sex
0 0 0 1 22 1
1 1 0 0 38 0
2 0 0 1 26 0
3 1 0 0 35 0
4 0 0 1 35 1

第三步:建模

  • 数据切分为train和test
In [4]:
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=33)
  • 使用决策树观察在检验集表现
In [6]:
from sklearn import tree
clf = tree.DecisionTreeClassifier(criterion='entropy', max_depth=3,min_samples_leaf=5)
clf = clf.fit(X_train,y_train)
print "准确率为:{:.2f}".format(clf.score(X_test,y_test))

准确率为:0.83

  • 观察各变量的重要性
In [7]:
clf.feature_importances_
Out[7]:

array([ 0.08398076, 0. , 0.23320717, 0.10534824, 0.57746383])

  • 使用更多指标来评估模型
In [17]:
from sklearn import metrics
def measure_performance(X,y,clf, show_accuracy=True,
show_classification_report=True,
show_confusion_matrix=True):
y_pred=clf.predict(X)
if show_accuracy:
print "Accuracy:{0:.3f}".format(metrics.accuracy_score(y,y_pred)),"\n"

if show_classification_report:
print "Classification report"
print metrics.classification_report(y,y_pred),"\n"

if show_confusion_matrix:
print "Confusion matrix"
print metrics.confusion_matrix(y,y_pred),"\n"

measure_performance(X_test,y_test,clf, show_classification_report=True, show_confusion_matrix=True)

Accuracy:0.834

Classification report
precision recall f1-score support

0 0.85 0.88 0.86 134
1 0.81 0.76 0.79 89

avg / total 0.83 0.83 0.83 223


Confusion matrix
[[118 16]
[ 21 68]]


  • 使用交叉验证来评估模型
In [8]:
from sklearn import cross_validation
scores1 = cross_validation.cross_val_score(clf, X, y, cv=10)
scores1
Out[8]:

array([ 0.82222222, 0.82222222, 0.7752809 , 0.87640449, 0.82022472,
0.76404494, 0.7752809 , 0.76404494, 0.83146067, 0.78409091])

第三步:决策树画图

  • 需要安装GraphViz’
In [9]:
import pydot,StringIO
dot_data = StringIO.StringIO()
In [14]:
tree.export_graphviz(clf, out_file=dot_data, feature_names=['age','sex','1st_class','2nd_class','3rd_class']) 
dot_data.getvalue()
pydot.graph_from_dot_data(dot_data.getvalue())
graph = pydot.graph_from_dot_data(dot_data.getvalue())
#graph.write_png('titanic.png')
#from IPython.core.display import Image
#Image(filename='titanic.png')

第四步:使用随机森林进行比较

In [11]:
from sklearn.ensemble import RandomForestClassifier
clf2 = RandomForestClassifier(n_estimators=1000,random_state=33)
clf2 = clf2.fit(X_train,y_train)
scores2 = cross_validation.cross_val_score(clf2,X, y, cv=10)
clf2.feature_importances_
Out[11]:

array([ 0.05526809, 0.02266161, 0.08156048, 0.46552672, 0.37498309])

In [12]:
scores2.mean(), scores1.mean()
Out[12]:

(0.81262938372488946, 0.80352769265690616)

python的数据科学资源

python和R是数据科学家手中两种最常用的工具,R已经介绍的太多了,后续我们来玩玩python吧。从出身来看,R是统计学家写的,python是计算机科学家写的,两者的出生背景不一样,随着数据爆发,python也慢慢发展,逐渐在数据科学中找到了一席之地。

包:
python也有非常多的扩展包,不过用于数据分析的并不象R那么品种繁多。常用的:
numpy:提供最基本的数值计算,使向量化计算成为可能。
scipy:提供了包括最优化在内的科学计算函数,不用自己写啦。
pandas:提供了类似dataframe的数据结构,处理表格数据非常方便。
matplotlib:画图必备,但用起来感觉不如ggplot2啊。
statsmodel:提供包括回归、检验等多种统计分析函数,python也能干R的活。
sklearn:数据挖掘必备,各种函数非常丰富,文档齐全,看得出CS出品就是不一样啊。

书:
python的数据方面书还不算很多,不过很有CS的味道,就是用show me the code,公式不多,这点我很喜欢。

现有可以找到的书基本上分为三类,一类是用基本语法实现统计分析和科学计算,例如下面的:
Think Stats
Think Bayes
A Primer on Scientific Programming with Python

另一类是以介绍一些包为目的,带有一些案例,例如:
Introduction to Python for Econometrics, Statistics and Data Analysis
numpy begin guide
Python for Data Analysis
matplotlib cookbook
Learning scikit learn Machine Learning
python text processing with nltk 2.0 cookbook
Social Network Analysis for Startups

最后一类是专门讲数据挖掘、机器学习的书:
programming collective intelligence:不用numpy包,只用基本语法实现一些算法
Machine Learning in Action:使用了numpy包,介绍了如何实现大部分算法
Machine learning an algorithmic perspective:体系非常完善,而且示例代码中使用了类

Notebook:
python的一大妙处就是ipython notebook,它可以把代码及其结果都存在一个网页上,方便分享学习。网上有非常多的notebook,其中成体系而又精彩的有下面三个:
https://github.com/datadave/GADS9-NYC-Spring2014-Lectures
http://blog.yhathq.com/posts/data-science-in-python-tutorial.html
http://slendermeans.org/pages/will-it-python.html

have fun!

python读入csv的三种方式

读数据到python有好几种方法,我们以读取iris.csv为例,将其中的数值部分提取出来。第一种方法是列表理解,文件读取到lines之后用一个嵌套的列表理解就可以将数值存为一个list。第二种方法是使用numpy库,它内带的loadtxt函数,读取的数据都认作是字符串,所以在第二行取我们需要的部分,并转为数值array。第三种方法是使用pandas库,它内带read_csv函数,读取数据会自动判断数值还是字符串,而且会自动保存好变量名,只需要用ix方法就可以类似R一样取出需要的子集,它存为dataframe对象。

这三种方法中最后一种最简单,不过花费时间比较长一点,第一种最麻烦,不过用时最短。这个可以通过ipython中的magic函数%%timeit来看。

python和ggplot2

python有个非常强大的工具,那就是ipython notebook。用户可以在浏览器中直接编写python脚本,并立即得到输出结果。这类文档可以存为ipynb分享给其它人,也可以存为html直接放在网站上,非常有利于学习交流。

在R语言方面就缺乏这类工具,不过ipython有一种“魔法”,可以在ipython中运行其它语言。在数据分析时,可以将python和R代码混编,充分利用两种语言的优势。以可视化为例,R的ggplot2图形语法可谓是独步江湖,python中虽然已经有不少优秀的绘图库。但总不及ggplot2用得习惯。下面的小例子就是示范在ipython notebook中画ggplot2。

首先是用numpy库建立两个向量,再用%load_ext建立python和R的连接机制。之后在ipython notebook的一个cell里面就可以使用%R后面接R代码行,或者使用%%R使用代码块。如果要在R代码块中读入python的对象,需要使用-i参数。

其它例子可以参见这个,python中也有人复制了ggplot语法,可参见这里

理解MCMC

在贝叶斯统计中,经常需要计算后验概率,概率计算就涉及到积分问题。一种解决方法是用解析式得到后验概率直接计算,另一种是利用统计模拟来计算近似值。 考虑一个简单问题,我们对一个硬币反复投掷,对于出现正面的概率theta先主观设定为一个均匀分布,然后实际投掷14次,得到11次正面,要根据这个信息data来更新后验概率。下面用统计模拟来计算。 因为实际的数据表现是正面出现次数较多,所以后验概率会做出相应调整。这个统计模拟的例子非常简单,容易实施。但如果涉及的参数个数很多时,计算量就会非常大。此时就可以尝试另一种统计模拟方法,即MCMC(Markov chain Monte Carlo)。我们先考虑一个思想实验:

一个岛国下面有7个岛,各岛人数不一样。旅行者已经在其中一个岛上,他要去各岛游历,准备在人口多的岛上游玩的时间长一些。但他并不知道具体的人口数,可以询问本岛和相邻岛的市长。他的旅行计划如下:
首先每天会扔一枚硬币,如果正面向上,就计划去左岛,如果反面向上,就计划去右岛。
然后比较当前所在岛人口数a和计划去的目的岛人口数b,若b小于a,则实际去的概率为b/a,若b大于a,则实际去的概率为1。所以实际去的概率归纳为min(1,b/a)。
就这样,旅行者每天不断的在不同岛间移动,最终他在各岛上呆的时间之比,就会收敛为各岛人口之比,这样就通过一个随机模拟得到了一个概率分布。

参考资料:Doing Bayesian Data Analysis with R