Image Style Transfer Using Keras and Tensorflow 使用Keras和Tensorflow生成风格转移图片

fast.ai有关Image Style Transfer的介绍是我见过最好的。Image Style Transfer实质上就是生成一副新的Image,使得它与给定Content Image的内容可以尽量接近,同时又兼具Stype Image的风格。为简单考虑,我们分开考虑如何分别生成Content Image和Style Image。不过,这里所讲的“生成”并不是train一个神经网络Model的参数,相反,我们会从一副随机生成的图片开始,通过L-BFGS的迭代优化算法来最小化这幅图片和给定Content Image之间的基于像素的MSE,此时,每次迭代都会利用梯度下降更新我们生成的图片。

最初提到Image Style Transfer的论文设定我们利用卷积神经网络,那么我们这里也选取熟悉的VGG16。为了更好的保留图片信息,我们需要将原始VGG16的max pooling层替换成average pooling层。

#Max pooling after Convolution
x = MaxPooling2D((2, 2), strides=(2, 2), name=’block1_pool’)(x)

#Average pooling after Convolution
x = AveragePooling2D((2, 2), strides=(2, 2), name=’block1_pool’)(x)

下面我们考察下如何生成Content Image。
#Step 0 — Create a pre-trained VGG16 model with last 3 FC layers removed
#and pooling layer replaced with average pooling 首先当然是创建一个VGG16网络,include_top=False表示舍弃后3层全连层
model = VGG16_Avg(include_top=False)

#Step 1 — Load Content Image and Image Preprocessing 接着加载Content图片,crop剪裁并无他意,由于我们并不是用VGG16来做图片分类,因此没有必要一定遵循VGG16规定的224乘224的输入图片尺寸
img = Image.open(‘JasonRidge.jpeg’)
img = img.crop(box=(30,60,430,460))
img

Oricontent

#Image Preprocessing — Subtract mean value and convert from RGB to BGR to feed VGG16 VGG16要求输入图片都减去ImageNet的均值并转化为BGR格式
rn_mean = np.array([123.68, 116.779, 103.939], dtype=np.float32)
preproc = lambda x: (x – rn_mean)[:, :, :, ::-1]
deproc = lambda x,s: np.clip(x.reshape(s)[:, :, :, ::-1] + rn_mean, 0, 255)

img_arr = preproc(np.expand_dims(np.array(img), 0))
shp = img_arr.shape

%matplotlib inline
import matplotlib.pyplot as plt
plt.imshow(np.squeeze(img_arr))

Preprocontent.png

#Step 2 — Build Neural Network Graph (Model)
from scipy.optimize import fmin_l_bfgs_b
from scipy.misc import imsave
from keras import metrics
#We take the first Convolution output of the 5th also the last Convolution Block
#We can latet verify with earlier layer 我们选取最后一个卷积块的第一层输出作为衡量标准
layer = model.get_layer(‘block5_conv1’).output
layer_model = Model(model.input, layer)
targ = K.variable(layer_model.predict(img_arr))

#Step 3 — Build Neural Network Loss and Gradient (Optimization)
#Loss is a scalar value — MSE 注意loss是一个标量,而梯度是一个ndarrdy
loss = metrics.mse(layer, targ)
grads = K.gradients(loss, model.input)
loss_fn = K.function([model.input], [loss])
grad_fn = K.function([model.input], grads)

#Create a container class for both Loss and Gradients Evaluator类可以在优化阶段方便计算Loss和梯度
class Evaluator(object):
def init(self, loss_f, grad_f, shp):
self.loss_f, self.grad_f, self.shp = loss_f, grad_f, shp

def loss(self, x):
loss_ = self.loss_f([x.reshape(self.shp)])
return np.array(loss_).astype(np.float64)

def grads(self, x):
grad_ = self.grad_f([x.reshape(self.shp)])
return np.array(grad_).flatten().astype(np.float64)
evaluator = Evaluator(loss_fn, grad_fn, shp)

#Step 4 — Time to Optimize
#Define a method which will optimize (minimize) the Loss using L_BFGS algorithm solve_image方法会将每次迭代更新的图片保存下来
def solve_image(eval_obj, niter, x):
for i in range(niter):
x, min_val, info = fmin_l_bfgs_b(eval_obj.loss, x.flatten(),
fprime=eval_obj.grads, maxfun=20)

x = np.clip(x, -127,127)
print(‘Current loss value:’, min_val)
imsave(‘results/Content_gen_conv5_1_{i}.png’.format(i=i+1), deproc(x.copy(), shp)[0])
return x

#Generate a random ‘noise’ image which will be the starting point of the optimization iteration 这里生成我们最初的图片,其实就是随机数组成的噪音图片
rand_img = lambda shape: np.random.uniform(-2.5, 2.5, shape)/100
x = rand_img(shp)
imsave(‘results/Content_gen_conv5_1_0.png’,x[0])

#Start the optimization with 10 iterations
iterations = 10
x = solve_image(evaluator, iterations, x)

('Current loss value:', array([ 92.49624634]))
('Current loss value:', array([ 38.59172058]))
('Current loss value:', array([ 25.30461121]))
('Current loss value:', array([ 19.56064987]))
('Current loss value:', array([ 16.33964348]))
('Current loss value:', array([ 14.33158779]))
('Current loss value:', array([ 12.92133713]))
('Current loss value:', array([ 11.83940792]))
('Current loss value:', array([ 10.99672127]))
('Current loss value:', array([ 10.28565693]))

Step 5 — Visulization

from IPython.display import HTML
from matplotlib import animation, rc

fig, ax = plt.subplots()
def animate(i):ax.imshow(Image.open(‘results/Content_gen_conv5_1_{i}.png’.format(i=i)))

anim = animation.FuncAnimation(fig, animate, frames=iterations+1, interval=500)
HTML(anim.to_html5_video())

可惜wordpress没法上传视频,不然那种从无到有做转变的一定会让你惊艳到。

Image.open(‘results/Content_gen_conv5_1_10.png’)

上面左图就是最初的噪音图片,右图则是10次优化后生成的图片。

有了生成Content Image的经验,我们就来考察下如何做Style Transfer。其实核心的改动就在Step2和Step3。下面就着重介绍下要点。

#Step 2 — Build Neural Network Graph (Model)

outputs = {l.name: l.output for l in model.layers}
#Content Model和上面完全一致,只是我们选了更早的卷积层输出,以期可以保留尽量多的图片内容
content_name = ‘block3_conv2’
content_layer = outputs[content_name]
content_model = Model(model.input, content_layer)
content_targ = K.variable(content_model.predict(content_arr))
#Style Model则不同,它需要5个卷积层的输出
style_layers = [outputs[‘block{}_conv2’.format(o)] for o in range(1,6)]
style_model = Model(model.input, style_layers)
style_targs = [K.variable(o) for o in style_model.predict(style_arr)]

#Step 3 — Build Neural Network Loss and Gradient (Optimization)
#First create a fucntion to calculate the Gram Matrix for Style Loss 下面新定义的方法专门计算格拉姆矩阵
def gram_matrix(x):
#We want each row to be a channel, and the columns to be flattened x,y locations
features = K.batch_flatten(K.permute_dimensions(x, (2, 0, 1)))
#The dot product of this with its transpose shows the correlation
#between each pair of channels
return K.dot(features, K.transpose(features)) / x.get_shape().num_elements()

def style_loss(x, targ): return metrics.mse(gram_matrix(x), gram_matrix(targ))

#Style Loss is a scalar
style_wgts = [0.05,0.2,0.2,0.25,0.3]
loss_sty = sum(style_loss(l1[0], l2[0])*w for l1,l2,w in zip(style_layers, style_targs, style_wgts))

#Content Loss is a scalar
loss_con = metrics.mse(content_layer, content_targ)

#Total Loss is the linear sum of Style Loss and Content Loss
#But here the denominator 10 is just used to balance whether
#Style or Content is more emphasized 这里,我们兼顾Content Loss和Style Loss的方法就是取2者的和,Content Loss下面的分母10是用来控制2者的权重,分母越大,则Content Loss的权重就越小,相应生成的图片也会更加侧重保留风格而非内容。
loss = loss_sty + loss_con/10

#K.gradients = tf.contrib.keras.backend.gradients 有了总的Loss,梯度的处理和上面Image处理也完全一致。
grads = K.gradients(loss, model.input)

#K.function = tf.contrib.keras.backend.function
#Note that loss_fn will return a scalar Loss while grad_fn a numpy array Gradients
loss_fn = K.function([model.input], [loss])
grad_fn = K.function([model.input], grads)

好了,那我们来看下最终的成果吧。

看下帅大叔和百花图以及小萌妹的搭配是什么效果吧。

 

 

Image Classification based on VGG16 Transfer Learning using Tensorflow 使用Tensorflow打造基于VGG16的图片分类器

时下最流行的机器学习非深度学习莫属。而想要事半功倍入门深度学习就非Transfer Learning莫属。fast.ai和Udacity的深度学习入门课程不约而同都以VGG16为基础,介绍了如何使用卷积神经网络实现图片分类器。我也参照着打造了以下专门识别男女的图片分类器。

作为基础的VGG16,其模型结构如下图所示。它先后由2个2层卷积、3个3层卷积和3个全连层组成,共16层。其中第1层和第2层卷积都是64个卷积核,所以输出的维度也相应变成(224,224,64),再经由(2,2)max pooling变成(112,112,64)。以后可如此类推。

vgg16

Transfer Learning 的核心思想就是,最大限度利用已经训练好的模型,为了实现特定分类功能,只需要舍弃原模型的最后一层或若干层,然后拼上新定义层并做fine tuning以获得新定义层的weights。当然,保留下来的原模型的所有层在fine tuning过程中均保持weights不变。这样显然可以显著减少模型的训练开销,并且也可以取得不错的结果,从而达到事半功倍的效果。

在fast.ai的Keras示例中,其做法如下代码所示。

model.pop() #舍弃当前模型的最后一层
for layer in model.layers: layer.trainable=False #当前模型的所有层的weights保持不变
model.add(Dense(num, activation=’softmax’)) #拼上新定义的全连层
self.compile() #使模型的新网络结构生效

而在Udacity的Tensorflow实例中,其做法如下代码所示。

feed_dict = {input_: images} #输入图片到VGG16模型
codes = sess.run(vgg.relu6, feed_dict=feed_dict) #截取VGG16模型的第一层全连层输出

以下代码重新定义一个2层全连层网络,输入即上面截取的VGG16模型的第一层全连层输出,而输出则为我们需要的特定分类器输出标签。

inputs_ = tf.placeholder(tf.float32, shape=[None, codes.shape[1]])
labels_ = tf.placeholder(tf.int64, shape=[None, labels_vecs.shape[1]])

fc = tf.contrib.layers.fully_connected(inputs_, 256)
logits = tf.contrib.layers.fully_connected(fc, labels_vecs.shape[1], activation_fn=None)

两者相比较,示例1中fine tuning是在包含VGG16前15层和新加1层全连层的模型上开展,当然真正可调的也只是最后一层新加的全连层的参数;示例2则需要先单独截取VGG16的前14层输出,再作为新定义的2层全连层的输入做模型训练。看似做法略有不同,其实本质上思想是一致的。

下面,看下图片经过VGG16不同卷积层的输出大概是什么样子的。原始图片是图像处理领域的“明星”, Lena。

第1层max pooling输出,共64个卷积核。可以看出不同的卷积核关注的焦点也不同,得到的输出也是各有侧重。

conl1

第2层max pooling输出,共128个卷积核。

convl2

第3层max pooling输出,共256个卷积核。

convl3

第5层max pooling输出,共512个卷积核。好吧,现在的输出每张image只有7乘7的像素。肉眼已经完全“沦陷”了。

convl5

那我们看下打造的性别分类器的输出吧。

Screen Shot 2017-08-07 at 23.35.36

哈,话说金星老师你即便还是他的时候,分类器还是可以“慧眼识珠”的。

Screen Shot 2017-08-07 at 23.35.52.png

额,下面这位。。。哈哈,我已无言。

Screen Shot 2017-08-07 at 23.37.41.png

至此,简单几十行代码便可以实现图片分类,并且还可以任意推广到其他图片分类应用。看来深度学习的魅力和潜力真不是盖的。那么,就继续加油吧。

fast.ai Deep Learning Courses – Making neural nets uncool again fast.ai 深度学习课程学习自省和小结

今天很开心发现fast.ai更新了深度学习的第2期课程。看来沉寂了几个月的深度学习又可以鼓噪起来了。接下来,我会定期记录深度学习的各种重点和启示。

真心佩服和感谢课程的讲师JEREMY HOWARD和RACHEL THOMAS,堪称当代AI福音传道者。

1期课程:http://course.fast.ai/lessons/lessons.html

2期课程:http://course.fast.ai/lessons/lessons2.html

Topic Extraction on Medium.com Documents using NMF and LDA in Python使用Python NMF和LDA对Medium文章主题提取

最近发现一个满是资源和干货的博客网站Medium.com,正好也在倒腾自然语言处理,决定写个网络爬虫做点主题提取分析。Medium.com可以按topic添加自己感兴趣的文章,比如下图的人工智能topic。Medium.com不会一次性加载所有的文章列表,而是通过用户下拉经由JavaScript动态加载。简单起见,我写的爬虫也就刚好爬取首次加载的文章。
Screen Shot 2017-07-31 at 14.57.42

# Define a Web Scraper which will get documents from Medium.com for a given topic
def mediumScraper(topic):
    links = []
    content = []
    
    response = requests.get('https://medium.com/topic/'+topic)
    soup = BeautifulSoup(response.text,'lxml')
    
    for post in soup.findAll("a", { "data-action" : "open-post" }):
        links.append(post['href'].split('?')[0])
    
    for link in links:
        response2 = requests.get(link)
        soup2 = BeautifulSoup(response2.text,'lxml')
        paragraphs = ''
        for paragraph in soup2.findAll("p", { "class" : "graf graf--p graf-after--p" }):
            paragraphs += paragraph.text

        content.append(paragraphs)
        print('Found one document on topic %s' %(topic))
        
    return pd.DataFrame({'topic':topic,'links':links,'content':content},columns=['topic','links','content'])

# Scrape 8 topics and return as a data frame
topics = ['data-science','politics','technology','sexuality','artificial-intelligence','culture','software-engineering','equality']
df = mediumScraper(topics[0])
for topic in topics[1:]:
    df = pd.concat([df, mediumScraper(topic)],ignore_index=True)

由于有些文章只针对会员开放,所以content内容为空。删除2篇空文章后,最终得到的data frame共200行,既每个topic有25篇文章。

接下来就是NLP的tokenization标准流程了。得到的corpus只包含去除非字符和stopwords并取字根的转小写单词。

corpus = []
for i in range(0,200):
    # Tokenization
    # Remove any character except letters
    # Convert into lower case
    content = re.sub('[^a-zA-Z]',' ',df['content'][i]).lower()
    # Convert into stem word
    ps = PorterStemmer()
    content = ' '.join([ps.stem(word) for word in content.split(' ') if not word in set(stopwords.words('english'))])
    corpus.append(content)

紧接着就可以对corpus做数值化矩阵转化。简单来说就是词频转化。

# 简单的词频转化;适用LDA
cv = CountVectorizer(max_features=2000)
X = cv.fit_transform(corpus)
# TF-IDF转化;适用NMF
tf = TfidfVectorizer(max_features=2000)
X_2 = tf.fit_transform(corpus)

说实话,Latent Dirichlet Allocation要透彻理解真不是件容易的事。简单来说,LDA是一个三层贝叶斯概率模型,包含词、主题和文档三层结构。

Non-Negative Matrix Factorisation非负矩阵分解则要直观许多。如下图,不过这里需要做下调整,V就是X或者X_2,一行是一篇文档,而一行包含的列就是文档包含的词;W就是一会作用LDA或者NMF(transform方法)得到的矩阵,它把列从词转变成主题,显然,文档可以被划归为所在行主题权值最大的那个;H就是作用LDA或者NMF得到的矩阵componets_,我们可以取每行词权值最大的几个作为该主题的一个概要。
Screen Shot 2017-07-31 at 14.14.30

下面我们可以比较下LDA和NMF的结果。

  • 上面8组是LDA得到的主题包含的top 50词的wordcloud。下面8组是NMF得到的主题包含的top 50词的wordcloud。看起来,NMF的分组更能贴近Medium.com本身的topic [‘data-science’,‘politics’,‘technology’,‘sexuality’,‘artificial-intelligence’,‘culture’,‘software-engineering’,‘equality’]。

 

  • 下面左图是LDA,右图是NMF。图片是8行25列,既每个像素代表一个文档,且每行的25列来源于df的同一个topic。像素的颜色则代表LDA或者NMF给出的topic。若能同行每列颜色一致,两行之间颜色异同,则越能接近Medium.com自身的topic给定。显然,NMF的表现优于LDA。

     

  • 另外,NMF的计算复杂度也明显小于LDA。
    LDA done in 1.403s.
    NMF done in 0.022s.

最后,小结一下,以SVD,NMF为代表的矩阵讲解运算在机器学习,图像识别,文字挖掘等诸多领域都有广泛应用。单就文本主题提取而言,NMF在这里由于LDA。不过实际使用中,不知道应该怎么取舍。哈,我会持续关注和提高的。

 

Predict Shipment Time in Supply Chain Using Random Forest Regressor 用随机森林回归预测供应链货运时长

最近参加了公司组织的Hackthon活动。其中一项的任务就是数字化供应量的货运执行与管理。作为切入点,我选择了考察货运的运输时间及其影响因素,并利用机器学习的回归算法train了一个模型,专门用来预测货运的运输时间。下面总结下过程中遇到的问题以及解决方法。

其实,最初我们是打算做一个分类器来预测货运的及时交付与否。然而,在于供应链的一线同事交流之后才发现,实际中的运营是约20%的货运可以保证on time delivery,剩余80%的只是on time delivery。所以他们并不是特别在意货运是否一定要及时交付。既然如此,我们也就改成做货运时间的回归预测。

遇到的第一个问题,就是数据集里面并没有直接给出货运所花的时间。不过,倒是可以找出特定carton的进出站记录。例如,下面截取carton ID 为#5##74@52的记录。1月8号准备运输到1月16号送达总共花费8天。
Screen Shot 2017-07-17 at 09.24.24
#简单起见,我们先对数据集按carton ID分组再取时间戳的最大值与最小值之差作为每个carton ID的运输时间。
def getTimeNMode(group):
try:
result = pd.Series([(max(group[‘DATE’])-min(group[‘DATE’])).days,group[‘MODE’].unique()[0]],index=[‘ShippingDays’,’ShippingMode’])
return result
except:
return pd.Series([0,”],index=[‘ShippingDays’,’ShippingMode’])
#group_keys=False可以将’CONTAINER ID’,’STATUS’,’DATE’,’MODE’任然保留为columns而不是index
df1 = df[[‘CONTAINER ID’,’STATUS’,’DATE’,’MODE’]].groupby(‘CONTAINER ID’,group_keys=False).apply(getTimeNMode)

下图是ShippingDays的分布图。显然,不同时间段的ShippingDays也会略有不同。
ShippingDays

有了lable值,接下来开始选取feature。这里自然少不了对供应链物流环节的业务知识。目前,公司依赖一套OTM系统来完成货运的运输安排,对应数据集里的OTM_ROUTE_CODE。比如,
CC-SJZ,从工厂直接运到客户
COC-SLCHONGKONG,从工厂先运到就近的物流中心,再运到客户
CDC-SLCROERMOND,从工厂先运到客户就近的物流中心,再运到客户
CODC-SLCHONGKONG-SLCROERMOND,从工厂先运到就近的物流中心,再运到客户就近的物流中心,最后运到客户
#提取OTM_ROUTE_CODE的3个字段,作为模型的3个feature
def getInTransitStopOvers(row):
try:
code = row[‘OTM_ROUTE_CODE’]
#code = row
codes = []
if len(code) > 0 and code.count(‘-‘) >0:
codes = code.split(‘-‘)
if codes[0] in [‘COC’,’CDC’]:
codes.append(‘NA’)
elif codes[0] == ‘CC’:
codes[1] = ‘NA’
codes.append(‘NA’)
else:
codes = [”,”,”]
return pd.Series(codes,index=[‘RouteType’,’OriSLC’,’DesSLC’])
except Exception, e:
print code, ‘ failed in parsing In-Transit StopOvers with the error ‘+str(e)
return pd.Series([”,”,”],index=[‘RouteType’,’OriSLC’,’DesSLC’])

df[‘OTM_ROUTE_CODE’].fillna(”, inplace=True)
df[[‘RouteType’,’OriSLC’,’DesSLC’]] = df.apply(getInTransitStopOvers,axis=1)

然后考察运输承运商集运输方式。为此,考察SHIPPING_METHOD_CODE的字段。例如,
000001_DHLC_P_1D中,DHLC表示承运商为DHL,运输方式是P(普通包裹),优先级为1D(1天到达)。
#提取SHIPPING_METHOD_CODE的3个字段,作为模型的3个feature
def getDeliveryData(row):
try:
method = row[‘SHIPPING_METHOD_CODE’]
#method = row
methods = []
if method.count(‘‘) >2:
methods = method.split(‘
‘,3)
else:
methods = [”,”,”,”]
return pd.Series(methods[1:],index=[‘Carrier’,’DeliveryType’,’DeliveryPriority’])
except Exception, e:
print method + ‘ failed in parsing Delivery Data with the error ‘+str(e)
return pd.Series([”,”,”],index=[‘Carrier’,’DeliveryType’,’DeliveryPriority’])

df[‘SHIPPING_METHOD_CODE’].fillna(”, inplace=True)
df[[‘Carrier’,’DeliveryType’,’DeliveryPriority’]] = df.apply(getDeliveryData, axis=1)

接着将carton的尺寸拆分为长宽高3个地段,作为另外3个feature。
def getCarbonDimensions(row):
try:
dimension = str(row[‘CART_DIM’])
#dimension = row
dimensions = []
if dimension.count(‘X’) >1:
dimensions = [float(d.strip()) for d in dimension.split(‘X’,3)]
else:
dimensions = [0.0,0.0,0.0]
return pd.Series(dimensions,index=[‘Length’,’Width’,’Height’])
except Exception, e:
print dimension, ‘ failed in parsing dimensions with the error ‘+str(e)
return pd.Series([0.0,0.0,0.0],index=[‘Length’,’Width’,’Height’])

df[‘CART_DIM’].fillna(0.0, inplace=True)
df[[‘Length’,’Width’,’Height’]] = df.apply(getCarbonDimensions, axis=1)

最终我们选取下面12个feature。
SHIP_FROM_ORG_ID
TO_CITY
WEIGHT
ORDERED_ITEM
Width
Length
Height
RouteType
OriSLC
DesSLC
Carrier
DeliveryType

我们选取的RandomForestRegressor模型要求所有的feature必须是数值型的。为此,需要先对Categorical类型的feature进行数值编码。这里,我也遇到本次Hackthon最大的难点。我发现sklearn自带的labelEncoder虽然好用,但是,接下来的难题着实耗费了很多时间。

  • 由于我们做的是offline training,需要把模型序列化成一个可导入的文件。这样的话,每个做了转化的faeture的labelEncoder都要单独序列化,甚是麻烦。
  • 在做预测的时候,我的理解是输入的值需要经由同样的labelEncoder做转化。然而,此时labelEncoder永远输出0。一阵吐槽纠结后才发现,labelEncoder的fit或者fit_transform好像只能针对一个array起作用。换句话说,labelEncoder把预测数据当成一个新的一维数组进行重现转化,所以输出值永远为0。最终,我写了一个小方法,可能存在问题,但是能达到我需要的目的。
    from sklearn.externals import joblib
    def categoryEncoding(col,ColName):
    colDict = {}
    colSet = set([i for i in col])
    i = 1
    for j in colSet:
    colDict[j] = i
    i += 1
    joblib.dump(colDict, ‘/home/felix/jupyter/data/dict_’+ColName)
    return colDictdef encodingLookup(ColName, val):
    colDict = joblib.load(‘/home/felix/jupyter/data/dict_’+ColName)
    if val in colDict:
    return colDict[val]
    else:
    return colDict[max(colDict)] +1

接下来就是水到渠成的training和tuning。
from sklearn.cross_validation import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestRegressor
from sklearn.grid_search import GridSearchCV
cols = [‘TO_CITY’,’DeliveryType’,’RouteType’,’OriSLC’,’DesSLC’,’Carrier’,’ORDERED_ITEM’]
for colName in cols:
#le = LabelEncoder()
#df_model[colName] = le.fit_transform(df_model[colName])
colDict = categoryEncoding(df_model[colName],colName)
df_model[colName] = df_model[colName].apply(lambda x:colDict[x])

Split into Training and Testing set

train_set, test_set = train_test_split(df_model, test_size=0.2)
X_train = train_set.iloc[:,1:]
y_train = train_set.iloc[:,0]

X_test = test_set.iloc[:,1:]
y_test = test_set.iloc[:,0]

param_grid = [{‘n_estimators’: [400,600,800], ‘max_depth’:[3,5,7]}]
forest_reg = RandomForestRegressor()
grid_search = GridSearchCV(forest_reg, param_grid, cv=5)
grid_search.fit(X_train, y_train)

最终,我选取n_estimators=500,性能如下。
R2 = 0.62546409754312426
相应的Feature Importances排序如下。
TO_CITY
0.26507
WEIGHT
0.191827
ORDERED_ITEM
0.148944
Width
0.075103
Length
0.065954
Height
0.065687
OriSLC
0.050146
Carrier
0.041274
DesSLC
0.040447
SHIP_FROM_ORG_ID
0.038849
DeliveryType
0.015408
RouteType
0.001292

完整的代码参加我的GitHub。经历这次Hackthon,以下几点思考:

  • 你以为的不一定就是业界关心的问题。比如,我们从最初的分类转成回归分析。
  • 字符转数值的labelEncoder在真正的线上系统是如何应用的?
  • 随机森林回归需要考虑regularization吗?
  • 从Feature Importances来看,长宽高貌似对货运时间的影响还比较大,值得再做深入挖掘。

Cisco Encrypted Traffic Analytics 思科机密流量分析

思科近日对外宣布下一代基于可编程、自适应和高安全的革新网络(DNA+ETA)。

其中,最引人注目就是思科在业界率先利用机器学习算法实现了对网络加密流量的监督学习及分类,从而可以在无需解密的前提下判定网络攻击(malware)。具体的实现细节尚未完全披露,不过可以从公开发布的白皮书里略窥一二。

nb-09-encrytd-traf-anlytcs-wp-cte-en

Data Insights into Provisioning Orders Using R and Python 使用R和Python数说Provisioning Orders

数据源来自于Oracle DB导出的2份数据文件。其中,orderList是order流水记录,二orderRequests是每个order的详细XML格式request payload。
首先,借助Python将2份文件读入Pandas Data Frame。然后遍历orderRequest利用xml.etree.ElementTree解析order包含的具体服务类型和license,并添加到orderList作为新的列。最终汇总得到19633行119维的数据集。数据清洗的部分此处先忽略。

  • pandas.DataFrame.diff
    在分析change order中license的变化时,需要比较同一客户的service license前后变化。如果用SQL,得在做完group by service之后再重新join查询得到license,然后再写DB function实现一个group里面的前后license值变化。简直如同噩梦,果断放弃。这一功能恰好可以直接在分组上使用Pandas的diff方法,只需对差值稍作转换即可得到想要的license变化类型,可谓简洁明了。
    大体的Python代码如下。
    #根据diff差值解释license是Upsell,Downsell还是Unchanged
    def licenseUpsellDownsellUnit(x):

if x > 0:
return ‘Upsell’
elif x 1:
group_df = group.sort_values(by=’ServiceOrderID’).reset_index()
group_df[serviceName+’.License Change’] = group_df[serviceName+’.TOTAL LICENSE VOLUME’].dropna().diff().apply(licenseUpsellDownsellUnit)
group = group_df.set_index(keys=’index’)
else:
group[serviceName+’.License Change’] = ‘Only One Order Found’
return group
#按照产品(这里是SiteURL)分类
def licenseUpsellDownsellAll(df_service, serviceName):
groups = df_service[[‘ServiceOrderID’,’TransactionType’,serviceTotalLicense,’OrderRecieveDate’]].groupby(df_service[‘SiteURL’])
licenseChangeGroups = groups.apply(licenseUpsellDownsellByGroup)
return licenseChangeGroups

  • 等到图形化license change数据的时候,我犯了一个小错误。
    如果我告诉你下面2张图我的本意实际想要传递同样的信息,你是不是会会对图1表示抓狂?!
    其实图1看起来凌乱,实则包含一个很重要的信息,就是同一产品的licens change频率,这点可以从曲线的长度“直观”体现。折线每次上升或者下降即表明一次Upsell与Downsell,所以折线的一段越长,代表保持的时间越长,反之,越短则表示license很快变化。而放弃曲线改为散点的图2就丢失了这部分信息。当然,图2的纵轴可以直接看到license change的百分比。所以嘛,并不是信息量越大的可视化才叫好,相反,可视化的目标不就正是希望同图表代替数据来说话吗,越少烧脑,才能更好传达信息,不是吗?
    图1
    Plot9CMRUpsellDownsell-2
    图2
    Plot9CMRUpsellDownsell-1
  • 往R dyplyr方法内传入列名变量
    这点上Python可以碾压R了。在Python里面,把DataFrame的列名作为函数入参传入是很正常的做法,但是R里面却没那么直观了。
    比如我想定义一个对order里面的服务按周分组的function,当我把serviceName传入的时候,正常使用dyplyr里面的n()和sum()来计数是不行的。这时要改用SE function, 即summarise变成summarise_。同时,由于summarise_里面既有特殊变量serviceName也有一般的totalCount和serviceBookingRatio,参考vignette(“nse”)文档使用interp。
    #Use as.name if you have a character string that gives a variable name
    interp(~ mean(var), var = as.name(“mpg”))
    #> ~mean(mpg)
    #or supply the quoted name directly
    interp(~ mean(var), var = quote(mpg))

#将Data.Frame的列名serviceName作为入参传入
plotServiceBookingTrendByWeek<- function(instanceName,serviceName) { aggregatedServices % filter(InstanceName==instanceName, TransactionType==’Provide’|TransactionType==’Trransfer’|TransactionType==’Renewal’,ProductModule!=”) %>% group_by(year,week,ProductModule) %>% summarise_(totalCount = ~n(), serviceBookingRatio = interp(~ sum(service,na.rm=TRUE)/totalCount, service = as.name(serviceName)))
… ..

  • 绘制order在世界地图上的分布图
    最简单应该是调用GoogleVis封装好的gvisGeoChart,可惜它自身的iso2c识别针对小国家不太理想,比如新加坡和香港在生成的图里面竟然没有显示。
    最终,通过引入rworldmap和countrycode这2个库成功绘制。

#Count of unique Webex Sites per each Country
webexAnnuityOrdersAggrbyCountry % filter(InstanceName==’SITE’, CountryCode != ”) %>% group_by(CountryCode) %>% summarise(WebExSitesCount = n_distinct(SiteURL)) %>% arrange(desc(WebExSitesCount))

#使用GoogleVis封装好的gvisGeoChart
webexAnnuityOrdersAggrbyCountry$CountryName <- countrycode(webexAnnuityOrdersAggrbyCountry$CountryCode,’iso2c’, ‘country.name’)
Geo=gvisGeoChart(webexAnnuityOrdersAggrbyCountry, locationvar=”CountryName”,
colorvar=”WebExSitesCount”,
options=list(projection=”kavrayskiy-vii”))
plot(Geo)

#使用rworldmap和countrycode绘制
annuityOrders_map_poly has both world map and annuityOrdersAggrbyCountry data merged together
webeAnnuityOrders_map <- joinCountryData2Map(webexAnnuityOrdersAggrbyCountry, joinCode = “ISO2”, nameJoinColumn = “CountryCode”)
webexAnnuityOrders_map_poly <- ggplot2::fortify(webeAnnuityOrders_map)
webexAnnuityOrders_map_poly <- merge(webexAnnuityOrders_map_poly, webeAnnuityOrders_map@data, by.x=”id”, by.y=”ADMIN”, all.x=TRUE)

#Ashmore and Cartier Islands作为澳大利亚的附属岛屿可以被忽略不计
webexAnnuityOrders_map_poly$WebExSitesCount[webexAnnuityOrders_map_poly$id==’Ashmore and Cartier Islands’] <- NA
webeAnnuityOrders_map@data$WebExSitesCount[webeAnnuityOrders_map@data$ADMIN==’Ashmore and Cartier Islands’] <- NA
webeAnnuityOrders_map@data[‘Text’] <- paste(webeAnnuityOrders_map@data$ADMIN,webeAnnuityOrders_map@data$WebExSitesCount,sep = ‘:’) ggplot() + coord_map(xlim = c(-180, 180), ylim = c(-60, 75)) + geom_polygon(data = webexAnnuityOrders_map_poly, aes(long, lat, group=group,fill=WebExSitesCount), color=’white’, size = 0.3) + scale_fill_gradient2(trans= ‘log’, breaks=c(0,10,100,200,4000), labels=c(‘4000′)) + xlab(NULL) + ylab(NULL) + guides(fill = guide_legend( title=’Number of Annunity WebEx Sites Around the World \n(June 1,2016 – June 1,2017)’, keywidth = 0.7, keyheight = 0.7, reverse=F, title.position=”top”)) + theme( plot.background = element_blank() ,panel.grid.major = element_blank() ,panel.grid.minor = element_blank() ,panel.border = element_blank() ,axis.ticks = element_blank() ,axis.text.x = element_blank() ,axis.text.y = element_blank() ,legend.position = “bottom” ,legend.direction = “horizontal” ) + geom_text(data=subset(webeAnnuityOrders_map@data,WebExSitesCount>0), aes(x=LON, y=LAT,label=Text), size=2.0)
最终图如下。
Plot1WorldMapofWebexSites-1

  •  Sankey绘制
    Sankey又叫RiverFlow,是一种反应不同类型转化的有效可视化形式。考查change order里面同一service的不同type的变化情况可以使用。

MCSessionSwappedCount % filter(TransactionType==’Amend’|TransactionType==’Transfer’,WEBEX.CONFERENCING.Session.Swap==’Swapped’) %>% group_by(WEBEX.CONFERENCING.FORMER.SESSION.ID,WEBEX.CONFERENCING.SESSION.ID) %>% summarise(count=n())
#这里的一个trick就是前后的名字必须不同,否则绘制会报错。简单的处理方法就是认为添加空格制造不同。
MCSessionSwappedCount$WEBEX.CONFERENCING.FORMER.SESSION.ID[MCSessionSwappedCount$WEBEX.CONFERENCING.FORMER.SESSION.ID==3] <- ‘Pro Meeting ‘
… …

MCSessionSwappedCount$WEBEX.CONFERENCING.SESSION.ID[MCSessionSwappedCount$WEBEX.CONFERENCING.SESSION.ID==3] <- ‘Pro Meeting’
… …

sk_MC <- gvisSankey(MCSessionSwappedCount, from=”WEBEX.CONFERENCING.FORMER.SESSION.ID”, to=”WEBEX.CONFERENCING.SESSION.ID”, weight=”count”)

print(sk_MC, tag = ‘chart’)

Sankey

  • Treemap绘制
    Treemap尤其适合展示多层级关系。这里我用它来显示提交order的合作商及其对应的系统账号的分布。
    treemap(subset(webexAnnuityOrdersAggrbyPartnerAccount,partnerAccountName!=’Direct’),
    index=c(“partnerName”,”partnerAccountName”),

vSize = “CCWOrderCount”,
title=”Partner Orders”,
type=”index”,
title.legend=”number of Partner orders”,
fontsize.labels=10,
align.labels = list(c(“center”, “center”), c(“right”, “bottom”)),
overlap.labels=1,
lowerbound.cex.labels=0,
inflate.labels=FALSE,
aspRatio=1.4)

Plot2PartnerTreeMap-1

当然,完整的分析还涉及到很多metric的考察,由于涉及到商业信息,就不在一一赘述了。完整的Rmd代码上传至我的Github账号。

Build A Movie Recommendation System using Python用Python打造一个电影推荐系统




In [1]:
#数据集包含100000条用户电影评分,其中9430被抽出来做测试。
#Collaborative Filtering Systems推荐系统使用Pearson Correlation Similarity Measure寻找相似度最高的用户。
#以下我尝试使用机器学习的聚类算法对电影进行分类;特定用户对特定电影的预测评分将是该用户对该电影所在聚类的均值评分。

import pandas as pd
import numpy as np 
import math 

import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.metrics import mean_squared_error

%matplotlib inline 


#读取Movies数据
names = ['movie_id','movie_title','release_date','video_release_date',
              'IMDb_URL','unknown','Action','Adventure','Animation',
              'Children','Comedy','Crime','Documentary','Drama','Fantasy', 
              'Film-Noir','Horror','Musical','Mystery','Romance','Sci-Fi', 
              'Thriller','War','Western']
movies = pd.read_table('Recommendation_data/u.item', sep='|', header=None, names=names)
movies_genre = movies.ix[:,6:]

#计算Movies类型矩阵的特征值和特征向量
mean_vec = np.mean(movies_genre, axis=0)
cov_mat = np.cov(movies_genre.T)
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
eig_pairs = [ (np.abs(eig_vals[i]),eig_vecs[:,i]) for i in range(len(eig_vals))]

eig_pairs.sort(key = lambda x: x[0], reverse= True)

tot = sum(eig_vals)
var_exp = [(i/tot)*100 for i in sorted(eig_vals, reverse=True)] 
cum_var_exp = np.cumsum(var_exp)

#选取主要成分的维数
plt.figure(figsize=(10, 5))
plt.bar(range(18), var_exp, alpha=0.3333, align='center', label='individual explained variance', color = 'g')
plt.step(range(18), cum_var_exp, where='mid',label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.show()
MovieRecommendationPCA
In [2]:
#Movies类型矩阵的维度从18维降至10维可以仍然保留90%的信息量
#应用PCA
movies_genre_array = movies_genre.as_matrix(columns=None)
pca = PCA(n_components=10)
pca.fit(movies_genre_array)
movies_genre_pca = pca.transform(movies_genre_array)
movies_genre_pca[:5]
Out[2]:
array([[ 0.80600899,  0.41113697,  0.48581643,  0.1860538 ,  0.39006354,
        -0.67679207,  0.0087999 ,  0.39217085, -0.2794839 ,  0.3694447 ],
       [ 0.34032144, -1.30358194, -0.08145307,  0.13841112, -0.41378311,
        -0.41717544,  0.06387956, -0.06770069,  0.05048498, -0.59779835],
       [ 0.15884273, -0.60255532, -0.08195047, -0.69319401,  0.20785322,
        -0.17105468,  0.17099464,  0.05925162,  0.14633722, -0.12079749],
       [-0.00931359,  0.26418695,  0.12291535,  0.30701683, -1.0959244 ,
        -0.08000186,  0.01655945,  0.31887652, -0.24933984,  0.01903654],
       [-0.67369651, -0.36994271, -0.05823788, -0.80094756, -0.19975073,
        -0.41947847, -0.63914919, -0.27094763, -0.15380729,  0.11432883]])
In [3]:
#选取K-Means的聚类族数
x =[]
y =[]
for i in range(100):
    movies_genre_cluster = KMeans(n_clusters=i+1, random_state=0).fit(movies_genre_pca)
    x.append(i+1)
    y.append(movies_genre_cluster.inertia_)
plt.figure(figsize=(10, 5))
plt.plot(x,y, 'go')
plt.ylabel('Inertia')
plt.xlabel('# of Clusters')
plt.legend(loc='best')
plt.show()
MovieRecommendationKMeans
In [4]:
#将PCA降维得到的新的Movies类型矩阵应用K-Means划分成50个聚类
movies_genre_cluster = KMeans(n_clusters=50, random_state=0).fit(movies_genre_pca)
movies_genre_cluster.labels_[:5]
Out[4]:
array([ 9, 23,  8, 19, 36], dtype=int32)
In [26]:
#读取用户及评分数据并将其与电影数据合并
movies['cluster_id'] = movies_genre_cluster.labels_
names = ['user_id','age','gender','occupation','zip_code']
users = pd.read_table('Recommendation_data/u.user', sep='|', header=None, names=names)

names = ['user_id','movie_id','rating','timestamp']
ratings = pd.read_table('Recommendation_data/u.base', sep='\t', header=None, names=names)

ratings_full = pd.merge(pd.merge(ratings,movies, on='movie_id'), users, on='user_id')

#计算每个用户对其评价过的电影所在聚类的均值矩阵
ratings_by_cluster = ratings_full.pivot_table(values='rating', index=['user_id', 'cluster_id'], aggfunc=np.mean)
ratings_by_cluster[:5]
Out[26]:
user_id  cluster_id
1        0             4.230769
         1             3.592593
         2             4.222222
         3             3.500000
         4             1.750000
Name: rating, dtype: float64
In [27]:
names = ['user_id','movie_id','rating','timestamp']
test = pd.read_table('Recommendation_data/u.test', sep='\t', header=None, names=names)

#返回用户对指定电影所在聚类的评分均值。如果用户尚未对指定聚类的任何电影给出评分,则返回NaN
def getAvgRatingByCluster(user_id,movie_id):
    avgRatingByCluster = float((ratings_by_cluster[user_id][movies[movies.movie_id==movie_id].cluster_id]))
    return avgRatingByCluster

#对测试数据计算预测值       
test['avgRatingByCluster'] = test.apply(lambda row: getAvgRatingByCluster(row['user_id'], row['movie_id']), axis=1)
test[:5]
Out[27]:
user_id movie_id rating timestamp avgRatingByCluster
0 1 20 4 887431883 3.583333
1 1 33 4 878542699 4.000000
2 1 61 4 878542420 4.230769
3 1 117 3 874965739 2.333333
4 1 155 2 878542201 3.000000
In [28]:
#其中有2320条数据为NaN
#对预测值不为空的计算MSE
print len(test[np.isnan(test.avgRatingByCluster)])
print np.sqrt(mean_squared_error(test[~np.isnan(test.avgRatingByCluster)].rating, test[~np.isnan(test.avgRatingByCluster)].avgRatingByCluster))
2320
1.1484145721
In [29]:
#为用户推荐尚未评分的电影。推荐电影来自于用户评分电影所属的电影聚类均分前三名所包含的所有用户尚未单独评分的电影。
ratings_by_cluster_df = pd.DataFrame(ratings_by_cluster)
def recommend_movies(user_id):
    return [list(movies[movies.movie_id==i].movie_title) 
            for i in movies[movies.cluster_id.isin((ratings_by_cluster_df.sort_values(by='rating', axis=0, ascending=False)).ix[user_id,][:3].index.values)].movie_id 
            if len(ratings_full[(ratings_full.user_id==user_id) &(ratings_full.movie_id==i)]) == 0]
    
recommend_movies(10)
Out[29]:
[['Net, The (1995)'],
 ['So I Married an Axe Murderer (1993)'],
 ['Akira (1988)'],
 ['Mimic (1997)'],
 ['FairyTale: A True Story (1997)'],
 ['Sphere (1998)'],
 ['Ghost (1990)'],
 ['Close Shave, A (1995)'],
 ['Old Yeller (1957)'],
 ['Parent Trap, The (1961)'],
 ['E.T. the Extra-Terrestrial (1982)'],
 ['To Catch a Thief (1955)'],
 ['Another Stakeout (1993)'],
 ['Body Snatchers (1993)'],
 ['Secret Garden, The (1993)'],
 ['So Dear to My Heart (1949)'],
 ['Arsenic and Old Lace (1944)'],
 ['Body Snatchers (1993)'],
 ['Dark City (1998)'],
 ['Fluke (1995)'],
 ['Lawnmower Man 2: Beyond Cyberspace (1996)'],
 ['Bogus (1996)'],
 ['Unforgettable (1996)'],
 ['Island of Dr. Moreau, The (1996)'],
 ['Charade (1963)'],
 ['Head Above Water (1996)'],
 ['Shiloh (1997)'],
 ['Little Princess, A (1995)'],
 ['Kim (1950)'],
 ['Virtuosity (1995)'],
 ['Robocop 3 (1993)'],
 ['Little Princess, The (1939)'],
 ['Tainted (1998)']]
In [30]:
#与皮尔逊相关PCS的对比
def pcs(user_1,user_2):
    m1 = ratings_full[ratings_full.user_id==user_1].movie_id
    m2 = ratings_full[ratings_full.user_id==user_2].movie_id
    l = list(set(m1).intersection(set(m2)))
    
    user1_avg = float(np.mean(list(ratings_full[ratings_full.user_id==user_1].rating)))
    user2_avg = float(np.mean(list(ratings_full[ratings_full.user_id==user_2].rating)))
    
    numerator = 0.0
    denominator1 = 0.0
    denominator2 = 0.0

    for i in l:
        r_1_i = float(ratings_full[(ratings_full.user_id==user_1) & (ratings_full.movie_id ==i)].rating)
        r_2_i = float(ratings_full[(ratings_full.user_id==user_2) & (ratings_full.movie_id ==i)].rating)
        numerator = numerator + (r_1_i - user1_avg)*(r_2_i - user2_avg)
        denominator1 = denominator1 + np.square(r_1_i - user1_avg)
        denominator2 = denominator2 + np.square(r_2_i - user2_avg)
    pcs = numerator/np.sqrt(denominator1)*np.sqrt(denominator2)
    return pcs


def getAvgRatingByPCS(user_id, movie_id, top=5):
    pcs_list = []
    for i in set(ratings_full.user_id):
        if i != user_id:
            pcs_i = pcs(user_id, i)
            if not math.isnan(pcs_i):
                pcs_list.append([i, pcs_i])
                
    top_pcs = sorted(pcs_list, key=lambda x:x[1], reverse=True)[:top]

    avgRatingByPCS = 0.0
    cnt = 0
    for i in range(top):
        pcs_rating= ratings_full[(ratings_full.user_id==top_pcs[i][0]) & (ratings_full.movie_id ==movie_id)].rating
        if len(pcs_rating) >0:
            avgRatingByPCS = avgRatingByPCS + float(pcs_rating)
            cnt = cnt+1
    if cnt != 0:
        avgRatingByPCS = avgRatingByPCS/cnt
    
    return avgRatingByPCS

print getAvgRatingByPCS(1,202)
print getAvgRatingByCluster(1,202)
4.4
3.76923076923
In [ ]:
#几点小结
#1. PCA加K-Means聚类得到最终预测评分并没有好于PCS。(MSE 1.15>1.06)
#2. 原因出在PCA丢失了信息?K-Means的聚类数量不够理想?
#3. 面对2320条NaN数据,可否用PCS的预测值,综合性能是否会好于单独PCS?
#4. 可否对用户+电影+评分一起进行聚类?
#5. 可否用监督学习来预测评分?性能会更好?

 

BEIJING PM2.5 Data Analysis in R – 用R分析北京PM2.5数据

library(ggplot2)
library(corrplot)
library(dplyr)

读取文件。

bjpm25 <- read.table("/Users/edwzhang/Downloads/PRSA_data_2010.1.1-2014.12.31.csv", header=TRUE, sep=",");

生成按年,按月,按周和按日的日期。

bjpm25$month[bjpm25$month<10] = paste('0',bjpm25$month[bjpm25$month<10], sep = "")
bjpm25$day[bjpm25$day<10] = paste('0',bjpm25$day[bjpm25$day<10], sep = "")
bjpm25['date'] = as.Date(do.call(paste0, bjpm25[c(2, 3, 4)]), "%Y%m%d")
bjpm25$yearDate = as.Date(cut(bjpm25$date, breaks = 'year'))
bjpm25$monthDate = as.Date(cut(bjpm25$date, breaks = 'month'))
bjpm25$weekDate = as.Date(cut(bjpm25$date, breaks = 'week', start.on.monday = TRUE))

考虑到存在全天PM2.5数据完全缺失,用当周均值来填补NA。

indx <- which(is.na(bjpm25$pm2.5), arr.ind = TRUE)
for (id in indx){
  bjpm25[id,'pm2.5'] <- mean( subset( bjpm25,weekDate==bjpm25[id,'weekDate'] )$pm2.5, na.rm=TRUE )
} 

按照PM2.5中国国家标准对数据添加分类。如果当时PM2.5计数是0,那么也划入一级。

PM2.5中国国家标准

bjpm25['pm25Cat'] = cut(bjpm25$pm2.5, c(0,50,100,150,200,300,1000), labels=c('I','II','III','IV','V','VI'))
bjpm25$pm25Cat[bjpm25$pm2.5==0] % count(yearDate, pm25Cat) %>%    
       mutate(pct=n/sum(n), ypos = cumsum(n) - 0.5*n),  
       aes(yearDate, n/24, fill=pm25Cat)) +
  geom_bar(stat="identity") +
  scale_x_date(date_labels = "%Y") +
  scale_fill_manual(values = c("green", "yellow", "orange","red",'purple','dark red'), guide = guide_legend(title = "PM2.5 Level")) +
  ggtitle("Beijing PM 2.5 per year") +
  coord_flip() +
  labs(y="Count of days",x="Year") +
  geom_text(aes(label=paste0(sprintf("%1.1f", pct*100),"%"), y=ypos/24))

总体来看,2010到2012年三年中空气质量逐年好转,但是2013年又急速恶化。2014年略有提高。

北京PM2.5按月分布图

ggplot(bjpm25 %>% count(monthDate, pm25Cat),  
       aes(monthDate, n/24, fill=pm25Cat)) +
  geom_bar(stat="identity") +
  scale_x_date(date_labels = "%Y/%m", date_breaks = "3 months") +
  scale_fill_manual(values = c("green", "yellow", "orange","red",'purple','dark red'), guide = guide_legend(title = "PM2.5 Level")) +
  ggtitle("Beijing PM 2.5 per month") +
  coord_flip() +
  labs(y="Count of days",x="Date") 

污染爆表(6类)天多发于10,11月和1,2月。

来看下污染爆表PM2.5达到令人发指的800以上分别有哪几天。

bjpm25$date[bjpm25$pm2.5>800]

##  [1] "2010-02-14" "2012-01-23" "2012-01-23" "2013-01-12" "2013-01-12"
##  [6] "2013-01-12" "2013-01-12" "2013-01-12" "2013-01-12" "2013-01-12"
## [11] "2013-01-12"
巧合的是,2010-02-14和2012-01-23分别是2010和2012农历新年的大年初一。

让我们重点关注下春节前后PM2.5的变化。

bjpm25_newyear <- subset(bjpm25, date %in% as.Date(c('2010-02-13','2010-02-14','2011-02-02','2011-02-03','2012-01-22','2012-01-23','2013-02-09','2013-02-10','2014-01-29','2014-01-30')))
bjpm25_newyear$date <- factor(bjpm25_newyear$date)

ggplot(data=bjpm25_newyear, aes(x=hour, y=pm2.5)) +
  geom_line() +
  ggtitle("Beijing PM 2.5 during Spring Festival") +
  facet_wrap(~date, nrow = 5)

有意思的是,2010到2012这3年每逢春节跨年PM2.5都明显飙升,但是2013年有明显收敛,甚至2014年还略有下降。

PM2.5按月直方图

ggplot(data=bjpm25, aes(x=pm2.5)) +
  geom_histogram() +
  facet_wrap(~month, nrow = 2)

总体而言,1,2,3,10,11,12月份的空气质量相对较差。

PM2.5按小时直方图

ggplot(data=bjpm25, aes(x=pm2.5)) +
  geom_histogram() +
  facet_wrap(~hour, nrow = 2)

总体而言,凌晨时空气质量优于白天。

PM2.5和累计风速分类箱图

ggplot(data=bjpm25, aes(x=IwsCat, y=pm2.5)) +
  geom_boxplot() 

随着风力增强,PM2.5均值值逐步下降,整体空气质量也向好。

PM2.5和风向箱图

ggplot(data=bjpm25, aes(x=cbwd, y=pm2.5)) +
  geom_boxplot() 

显然,刮西北风的时候空气质量较好。

PM2.5,累计风速分类和风向散点图

ggplot(data=bjpm25, aes(x=date,y=pm2.5, size =IwsCat , color = cbwd)) +
  geom_point() +
  scale_x_date(date_breaks = "3 month")

再次应证沉在下面的多是风力较强的西北风,同时期的空气质量相对较好。

相关系数

corrplot(cor(bjpm25[,c('pm2.5','DEWP','TEMP','PRES','Iws','Ir','Is')]), method = "circle")

相关系数分析也显示累计风速与PM2.5呈现较强的反向关系。