落园 |专注经济视角下的互联网

统计学习精要(The Elements of Statistical Learning)课堂笔记(二十五):降维和PCA

 降维

降维完全属于unsupervised learning了,即给定数据集\{x_{1},...,x_{n}\},x_{i}\in\mathbb{R}^{p},我们希望降到q维的\{z_{1},...,z_{n}\}。从这个角度来讲,降维和聚类还是有相通之处的,都是对于特征的提取。只是一个从行的角度出发,一个对列操作的感觉。

PCA(主成分分析,Principle Component Analysis)

个人觉得这也是起名字起的比较好的模型之一...乍一听起来很有用的感觉 -_-||

1. 求u_{1},\left\Vert u_{1}\right\Vert =1使得z_{i}=u_{i}x_{i},且\sum z_{i}^{2}最大。

PCA

直觉上来讲,就是想寻找一个主方向。

这样,求解问题为:

\max\sum z_{i}^{2}=\sum(u_{1}x_{i})^{2}=\sum u_{1}'x_{i}x_{i}'u_{1}=u_{1}'\left(\sum x_{i}x_{i}'\right)u_{1}。所以我们只需要求一阶导数\frac{\partial\sum z_{i}^{2}}{\partial u_{1}}即可。

设A为对称矩阵,则存在正交阵u'u=I使得A=u\Lambda u',其中\Lambda=\left[\begin{array}{ccc}\lambda_{1}\\ & \ddots\\ & & \lambda_{n}\end{array}\right]为A的特征值矩阵,故Au_{i}=\lambda_{i}u_{i}(列向量为特征向量)。不失一般性,我们可以排序使得\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{n}(从大到小排序)。

最大特征值: \max_{\left\Vert x\right\Vert =1}x'Ax=\lambda_{1}

同时C=\frac{1}{n}\sum x_{i}x_{i}'为x的相关矩阵,\{u_{1},u_{2},...,u_{q}\}=U_{q},从而z_{i}=U_{q}U_{q}'x_{i},\forall i

2. 找到\{u_{1},u_{2},...,u_{q}\}=U_{q}(q维的子空间)

x_{i}投影到该q维空间,这样z_{i}=U_{q}U_{q}'x_{i}=U_{q}\left[\begin{array}{c}u_{1}'\\\vdots\\u_{q}'\end{array}\right]x_{i}=U_{q}\left[\begin{array}{c}u_{1}'x_{i}\\\vdots\\u_{q}'x_{i}\end{array}\right]=\left[u_{1},u_{2},...,u_{q}\right]\left[\begin{array}{c}u_{1}'x_{i}\\\vdots\\u_{q}'x_{i}\end{array}\right]=(u_{1}'x_{i})u_{1}+...+(u_{q}'x_{i})u_{q},且\left\Vert u_{q}u_{q}'x-x\right\Vert _{F}^{2}最小。

A矩阵的范数:\left\Vert A\right\Vert _{F}=\sum_{i}\sum_{j}a_{ij}^{2}=tr(AA')
tr表示矩阵的迹(对角线元素和)。

则上述问题等价于,求(u_{1},...,u_{q})=U_{q}使得tr[(u_{q}u_{q}'x-x)'(u_{q}u_{q}'x-x)]最小。

=tr(x'U_{q}U_{q}'U_{q}U_{q}'x-x'U_{q}U_{q}'x-x'U_{q}U_{q}'x+x'x)=tr(-x'U_{q}U_{q}'x+x'x)最小。

即使得tr(x'U_{q}U_{q}'x)最大(注意没有负号)。

S=x'x称为数据的相似矩阵=(x_{i}'x_{j})

AA'A'A均为对称阵,且两个阵有相同的特征值。记r为A的秩,AA'的特征向量u_{1},...,u_{m}\in\mathbb{R}^{m},A'A的特征向量v_{1},...,v_{m}\in\mathbb{R}^{n},则A'u_{i}=\lambda_{i}u_{i}A'v_{i}=\lambda_{i}v_{i}。做奇异值分解,则A=U\Lambda V'.

由此,tr(x'U_{q}U_{q}'x)求得的和前述结果等价。

回到PCA。如果降维后需要重构,则u_{i}z_{i}=\hat{x_{i}},解\min\sum_{i}\left\Vert x_{i}-\hat{x_{i}}\right\Vert 即可。

3. 对偶PCA。如果p\geq n即数据非常高的时候,可以转置后再做。

4. KPCA (kernel)PCA也可以先用核函数\Phi(\cdot),即实现非线性的降维。需要注意,降维的过程需要保持可逆。

---------------

PS. PCA不适合解决overfitting的问题。如果需要解决,加regularization项。


统计学习精要(The Elements of Statistical Learning)课堂笔记(二十二):核函数和核方法

补上笔记。这节课讲的就是大名鼎鼎的Kernel Method...

核函数(正定)

定义 K(x,y), x,y\in\mathbb{R}满足:

1) 对称: K(x,y)=K(y,x)

2) 正定: n个观测x_{1},x_{2},...,x_{n}\in\mathbb{R}^{p}K_{n}=\left[\begin{array}{ccc}K(x_{1},x_{1}) & \cdots & K(x_{1},x_{n})\\\vdots & \ddots & \vdots\\K(x_{n},x_{1}) & \cdots & K(x_{n},x_{n})\end{array}\right] 正定(或者非负定)。

K(x,y)举例:

  • 常数——K(x,y)=C\Rightarrow\sum_{j}\sum_{j}cu_{i}u_{j}=c\left\Vert u_{i}\right\Vert
  • 内积—— K(x,y)=\sum x_{i}y_{i},或广义下K(x,y)=(\Phi(x),\Phi(y)),其中\Phi(x):x\rightarrow X,从\mathbb{R}^{p}\rightarrow\mathbb{R}^{q}

性质:

1. 封闭性

1) K(x,y)正定,\alpha>0,则\alpha K(x,y)正定。

2) K_{1}(x,y)正定,K_{2}(x,y)正定,则K_{1}(x,y)+K_{2}(x,y)正定,K_{1}(x,y)\cdot K_{2}(x,y)正定。

3) \{K(x,y)\}正定,K_{n}(x,y)\rightarrow K(x,y),则K(x,y)正定。

4) (1+(x,y))^{k}正定

5) \exp(-\frac{\left\Vert x-y\right\Vert }{2\sigma^{2}})正定。

2. 归一性

\bar{K}(x,y)=\frac{K(x,y)}{\sqrt{k(x,x)}\sqrt{k(y,y)}}正定,\Rightarrow\bar{K}(x,x)=1

再生核Hilbert空间(RKHS)

(走神一下:关于这个命名的吐槽猛击 -> 翻译版、 英文原版Normal Deviate

1. Hilbert空间:完备内积空间,可以视作欧氏空间的推广。H=\{x,y,z,...\}

在这个空间中,我们定义:

  • 加法:x+y
  • 数乘:\alpha x, \alpha\in\mathbb{R}
  • 内积(x,y):对称性(x,y)=(y,x);线性 (x_{1}+x_{2},y)=(x_{1},y)+(x_{2},y)\alpha(x,y)=(\alpha x,y).
  • 零元素:若(x,x)=0,则x=\phi定义为零元素。
  • 完备性:如果x_{n}\rightarrow xx_{n}\in H,则x\in H。(收敛到该空间内)。

2. 再生核Hilbert空间

给定K(x,y)正定,可以构造Hilbert空间H使得K(\cdot,y)\in H(K(\cdot,y),K(\cdot,z))=K(y,z);且构造一个\Phi(x):\mathbb{R}^{p}\rightarrow H,使得K(x,y)=(\Phi(x),\Phi(y)),即核函数可以写成内积形式。

这样对于\forall f\in H(f,K(\cdot,x))=f(x)

核方法

1. 基本思想

将线性模型推广到非线性模型的方法(其中较为简单的一种)

x\underrightarrow{\Phi(x)}\tilde{X}=\Phi(x),从\mathbb{R}^{p}\mathbb{R}^{q}(H)的一个映射。举例:\Phi(x)=(x,x^{2}),这样就可以拓展为广义线性模型。

2. SVM

\min\frac{1}{2}\left\Vert w\right\Vert ^{2}+C\sum_{i}\xi_{i}

s.t.\, y_{i}(w_{i}+b)\geq1-\xi_{i},\forall i

可以转化为:

\min-g(\lambda,\mu)=\frac{1}{2}\sum_{i}\sum_{j}(\lambda_{i}y_{i})(\lambda_{j}y_{j})(x_{i}'x_{j})-\sum_{i}\lambda_{i}

s.t.\sum_{i}\lambda_{i}y_{i}=0

0\leq\lambda_{i}\leq C

\mu_{i}\geq0

w=\sum\lambda_{i}x_{i}y_{i}b=y_{i}-w'x_{i},则f(x)=sign(\sum\lambda_{i}y_{i}(x'x_{i})+b)

非线性变换之后,

\min\frac{1}{2}\left\Vert w\right\Vert ^{2}+C\sum_{i}\xi_{i}

s.t.\, y_{i}(w_{i}+b)\geq1-\xi_{i},\forall i

注意此时w的维数有变化(p\rightarrow q)。

---------------------

如果各位更关心SVM后面的直觉,还是去看看Andrew Ng的相关课程吧...这里推导太多,直觉反而丢了一些。


从R里面底层操纵Excel/xlsx(自动化报告福音)

好吧,我在eBay折腾的最多的就是生成自动化报告时候各种软件之间的相互调用,什么R啊,SAS啊,Teradata啊,Excel啊,Python啊,反正基本都有机会相互调用一下。每到此时我就深深感慨选择一个library丰富的工具是多么的重要!You could hardly expect what you colleagues are handy with!(P.s. 不要跟我提VBA这种逆天存在的东西。有哪个时间研究它你学点啥别的不好...)

今天忍无可忍+心情大好的折腾了一下R和excel。这个不是简单的从R里面读写excel数据,而是真心用R去操纵excel里面的单元格(cell),除了读写数据之外还要定义样式什么的。excel作为一个奇葩的软件,you may never expect where people would paste data to! 然后他们再自定义一堆样式(我恨这种点点鼠标就能改的东西,你丫又不是Photoshop...)。

但是没办法,人家定义好的“高端洋气”的报表姿态你不能轻易动啊。只能乖乖的往里面paste数据。这件事虽说一次两次手动也就罢了,三五次真的是要疯掉的。anyway,万事总有解决的途径...

很久以前从Yixuan 的博客上得知有xlsx这么个包,当时只记得这东西可以读写xlsx...直到后面折腾了一下才知道这货底层居然调用的是java的xlsx API,也就是说不用写Java也可以操作xlsx了,yeah!

为了生成excel格式的自动化报告(不要问我为啥不用knitr,不用***,说起来都是泪呀!),我主要需要解决的就是:

  • 读取原有xlsx文件,保持格式、附加新格式。
  • 在相应的位置粘进去新的数据。(当然如果只有这么一个需求可以通过ODBC来做...)

第一个倒是满简单的,就是较之yixuan代码里面的createWorkbook(),改成loadWorkbook()就可以了。然后就是找到相应的sheet,这个也满简单的,一行getSheets搞定。

然后第二步建议不要去操作cell(太没效率了),直接操作cellblock。CellBlock()可以用来定义一个新的CellBlock,然后灵活运用CB.setBorder()和CB.setColData()就可以先增加边框、然后一列列填充数据。这里使用按列填充数据主要是因为R里面的Data Frame是一列一个数据格式的,一下子把一块儿都paste到excel的cellblock里面的话,会报错...BTW为了定义边框的样式,需要用到Border()。类似的还可以定义Fill和Font这些。

同上,最好不要直接用addDataFrame()来直接贴数据...格式不能覆盖。如果是要在一个新的sheet上贴数据,那么就write.xlsx(sheetName="newsheet",append=T)好了。不需要通过上述底层的API折腾了。

最后还有一个比较有用的函数,autoSizeColumn()可以用来自动调整列宽。全鼓捣完之后saveWorkbook()保存就可以啦。

最后的最后,一个珍贵的建议——都在R里面把数据整理好再去想输出到excel里面(什么reshape2啊,data.table啊,plyr啊,该上的一起上啊!),千万别手贱在excel里面改一点点小东西...每一次都手动改一下下你的时间就被白白浪费了好几分钟!珍爱生命,远离excel...

附上一段我最后搞定自动化报告的代码:

 



长安一片月,万户捣衣声——西安速记

不知道为什么,每次西安都是来去匆匆。没有呆过超过三天的。于是乎,对于这座我最喜欢的北方城市,印象总是断断续续的,一丝一缕的。

DSC07959-r
长安万户

一直在想,我可能再也不会经历早起赶飞机、暮迟宿夜车这样的旅行了,毕竟已经年老体衰、精神不比当年了。然而主观臆断就是用来被打破的。几句电话、几往短信,然后我就背着背包出现在了大雨滂沱的古城西安(小轩哥霸气~~~一硕下了班就直接来了,随性至此,V5~~~)。那天的雨着实太大(两朵乌云着实不能出现在同一座城市),一切的记忆都是湿漉漉的。那种酣畅淋漓的大雨,正是北方下雨的记忆——或许慢慢习惯了南方的阴雨绵绵,这样久违的畅快一下子勾起了北方的亲切。

“少小离家老大回,乡音未改鬓毛衰。”

南方实在是缺乏一点气势。总是柔柔的,柔柔的。时间久了,就觉得整个人都靡软了,缺乏奋斗的力气了。然而原本,尚年少,可疯狂。

长安是一座有脾气的城市。火车站正对着的城墙,相比于苏州火车站正对着的园林,气势上远远压过。一下子,就不知道自己身处的年代——为什么岁月可以这样变化,穿越古今?骑行在宽阔的城墙上,正如丞相当年那般气宇轩昂——两侧的车水马龙凡人喧嚣,与我何干?沏一壶清茶,落座于伞下。携三五好友,谈天论地,好不畅快!

DSC07933-r
冲天香阵透长安,满城尽带黄金甲。

DSC08070-r
长安陌上无穷树,唯有垂杨管别离。

金戈铁马,唤起的是对于远古的铁血记忆。许多奇迹,长埋地底。

DSC07881-r
秋风吹渭水,落叶满长安。

偶尔,月夜朦胧,恍若江南。停辙,起伞,落座。

DSC07868-r
长相思,在长安

DSC07846-r
曲江烟雨

回程之前,随意漫步。东北角楼,斜阳无数。偶逢高人,少聊摄影。略略指点,方知陕西之博大,意欲驾车重游。

DSC08083-r
古城落日

来长安,只要吃到皮脆里嫩的肉夹馍,以及酸酸辣辣的凉皮,人生就满足了。太解馋了。(战绩:5.5个肉夹馍,吼吼)

DSC08000-e
一只边走边吃的吃货的幸福

DSC07986-r
曲径求“道”

长安漫游路线图:钟楼/鼓楼(有定时钟鼓表演) -> 曲江池遗迹 -> 兵马俑/秦始皇陵 -> 城墙(骑行,南门有定时士兵表演) -> 书院门漫步 -> 东北角楼落日。下面是随手拍的照片,随便看看,权作回忆。
[......]

Read more


数据分析职业病

分析是种职业病,融贯在每一分血液里,每一分骨髓里...去参加个Qcon看看人家的创意网站,然后心里各种痒痒,拉着思喆哥、堰平兄饭局讨论实现的原理也就罢了,最近只要一出门就习惯性的开始思考某些稍微“违背常理”或者“略显聪明”的现象...

比如这次去西安,在去的航班上,就开始思考起来了“航空公司的数据挖掘”....
-------------------回忆的分割线-------------
有些事情纯属职业病。这次上海飞西安坐的是南航的航班,一路折腾到飞机上就已经疲惫不已了,直接睡过去。

后面的一切毫无波折,如果不是临下飞机十几分钟的一段对话,我估计会对这段航行毫无记忆。只是突然间空姐的一句问候,"您是***女士么(我的本名)",让我第一反应是我不是你们的金银卡啊,这也开始问候了?...然后笑意盈盈的递给我一张会员申请表。"您虽然还不是我们的会员,但是您是我们的潜在会员。所以请您加入我们南航明珠俱乐部"...

我第一反应是,"潜在会员"这个是怎么分析出来的?目测我大概有一年的时间没有飞过南航(在过去的半年时间我似乎都完全没有飞过),难道他们有个"沉睡用户苏醒监测"?要不就是,正巧这次我累计乘坐南航的次数达到了他们分析的阈值(比如10次)?要不就是,每次坐南航我都累积东航,让他们终于忍无可忍了...还是说,他们真有一个customer life value model,算出来我对他们的潜在价值了?

蛮有意思的是,我曾经有段时间周周飞南航,他们却从来对我爱理不理,所以我猜测他们的模型里面对于"reactivated"的顾客有个很高的权重。

到最后,南航猜的准吗?准,有可能是我确实在东航还有一些里程可以挥霍。不准,则是我现在飞行大都是私人旅行了,基本不可能象以前做咨询时候出差那般频率了,所以我的潜在价值肯定没有模型上估计出来的高。如果这个模型进一步分析"公务旅行"还是"私人旅行",怕是效果会更好吧。不知道能不能通过机票代理来区分这两类客源...所以,准,却有点晚了。

当然对于垄断国企来说,这个CRM并不是那么重要,反正利润会一直在那里,国内的常客计划也发展不起来。只是感慨一下,这样的分析要是想做好,绝对离不开自己对于这项服务的亲身体验。只有飞的多了,才知道常客计划的最大引力和最关键时点。而后的分析,才会水到渠成吧?

职业病发作完毕...