多元统计是统计学的一个重要分支，研究同时发生的多重随机变量的统计性质及各变量间的相互关系。多元统计的应用称多元分析，传统的多元分析可以看作是各个应用学科对多变量问题的研究而发展起来的各种方法的总和。其中大部分来自社会科学领域，特别是教育学和心理学，比如因子分析，主成分分析，对应分析，典型相关分析和多维标度法。还有一些方法来自早期统计学的应用实践，比如线性判别分析来源于分类学，而多元方差分析来源于农作物生长的随机试验，可以追溯到R.A.Fisher在农业试验站的工作。回归分析则起源于遗传学和优生学的研究。这种多样化的来源与实践是多元统计在很长时间内保持发展活力的重要原因。

     但是在今天，随着数据规模越来越巨大，数据分析问题的复杂程度越来越深以及对数据的储存分析所需的计算能力要求越来越高，传统的基于正态假设的，以研究数据协方差结构为主的多元统计遇到了越来越多的挑战。大部分传统的多元统计方法创造于较早的年代，能处理的数据规模有限，在计算能力上也有很大的制约。

      在传统上应用多元统计的各个学科按照各自面对的问题，给出了很多的新的思想和方法，特别是在计算机和信息科学领域，系统化的发展起了数据挖掘和机器学习这两个新的分支学科。这两个学科发展了很多传统的多元统计方法，也增加了很多新方法。这些新的发展促使多元统计扩展自己的框架以适应新的挑战，这就是作者Izenman写作本书的出发点。

      这本书是作者按照统计学习观点重新审视多元统计的结果。本书基本可分为四个部分：多元随机变量的统计性质；回归和有监督学习方法；聚类；非线性维数缩减。包括一些传统的如回归分析这样的主题，还有一些如流形学习这样比较新颖的主题。内容上，突出了各种方法的背景和算法，并给出了很多案例。虽称不上经典，可也是一本内容丰富很值得一读的著作。
我打算以这本书的框架写一些读书笔记。一方面用来串一下自己的统计和机器学习的知识，
另一方面可以由此系统一下R在多元统计和统计学习的应用。

——————————————— 书评的分割线——————————————————

   第一篇是一个相对（回归聚类什么的）冷僻的方法：多维标度法（MDS）。（本书第十三章）

1.MDS简介
多维标度法（Multidimensional Scaling），是一种维数缩减方法，把高维的数据点映射到一个低维的流形上；同时也是一种可视化方法，实践中通常利用2D或3D的MDS 结果观察（投影后）点的分布和聚集来研究数据的性质。


      简单来说，MDS要处理的问题是：由n个指标（变量）反应的实体（entity），仅知它们之间的某种距离（相异度）或相似度，如何在较低维的流形中推测实体间的原始距离，以反映这n个实体的真实结构关系。

   MDS起源于心理测量，具有广泛的应用，比如：
      心理学:研究不同类别的心理刺激（如人格特质，性别角色）或物理刺激（如面孔，声音，颜色，味道）的认知的潜在结构，并绘制这些刺激的“感知图”（perceptual map）
     市场研究：研究消费者的产品选择和产品偏好，可以识别产品间的联系。
     社交网络：识别大型网络的集群。
以及地理学，生态学，天文学，分子生物学，计算化学，图形学甚至流行音乐研究等等。

    MDS不是一个单独的方法，而是有相似思想的不同算法的集合。常用的MDS为距离标度（distance scaling），可以分为度量标度（包括经典标度和最小二乘标度等）和非度量标度。MDS方法和自组织映射（SOM）有着相似的目标。
MDS的思想方法在随机森林中也起着重要的作用，同时启发了非线性流形学习。

一个简单的例子和一些术语

一个经典的例子。关于MDS的易于理解的经典案例往往来自地理方面(本例最早来自文献[4]）：
关于美国十个城市间飞行里程的数据。这十个城市是：亚特兰大（Atlanta）、芝加哥（Chicago）、丹佛（Denver）、休斯顿（Houston）、洛杉矶（Los Angeles）、迈阿密（Miami）、纽约（New York）、旧金山（Sanfran）、西雅图（Seattle）和华盛顿（Wash•DC）。用MDS分析这个数据集，绘出一张“真实”的地图。数据取自http://lem.gdcc.edu.cn/jiaox/epd/contents/softstudy/sas/jlfx01-left1.HTM。
把数据按照表格中的形式整理为csv文件。使用R函数cmdscale（stats包，使用经典标度）

city<-read.csv('airline.csv',header=TRUE)
city1<-city[,-1]#这个数据集的第一列是名字，先把它去掉
for (i in 1:9)
for(j in(i+1):10)
city1[i,j]=city1[j,i] #把上三角部分补足
rownames(city1)<-colnames(city1) #再把行名加回来
city2<-as.dist(city1, diag = TRUE, upper = TRUE)#转换为dist类型
city3<-as.matrix(city2) #转化为矩阵
citys<-cmdscale(city3,k=2) #计算MDS，为可视化，取前两个主坐标plot(citys[,1],citys[,2],type='n') #绘图
text(citys[,1],citys[,2],labels(city2),cex=.7)#标上城市名字
看看图和真实的地图方向(东西海岸)是反的阿，没关系，修改一下：
plot(-citys[,1],citys[,2],type='n')
text(-citys[,1],citys[,2],labels(city2),cex=.7)
现在就是真实地图的结构了。

下面用ggmap包画张地图，对比一下：

MDS建立于实体(entity)的成对比较之上。实体的接近程度称为“邻近”（Proximity），
邻近可以是定量的，如真实实体的物理距离，对顺序尺度（ordinal scale）的主观判断，也可以是定性的，如在心理测量中关于相似或者相异程度的主观评分。
设δ[i，j]表示i，j两个实体的相异程度，对全部的n个实体X[1]…X[n]，我们可以得到矩阵
Δ =（δ[i，j]），这称为距离矩阵（Proximity Matrices）。这个矩阵是通常写为一个下三角的非负的矩阵（主对角线一般为0）。在可度量（metric）距离之下，满足三角不等式。（注：在某些心理测量实验中，距离矩阵可能是非对称的，一种调整的方法是取对角线两侧的对应元素均值，使之对称）。关于距离矩阵的构建在某些问题里可能相当复杂。可参考[2]。

MDS的目的是寻找一个低维（ t维）空间，其上对应n个实体的点为Y[1]，…,Y[n]，
他们的距离d[i,j]与δ[i，j]尽量接近。Y 称为MDS的一个解，在进行可视化时也可以称为拟合构图。所以，MDS是在寻找关于解的优化算法。

2.经典标度法（Classical Scaling）

经典标度法是MDS系列算法中最早出现的一种。尽管经典标度法可视为最小二乘标度法在某种情形下的一个特例，探讨经典标度法的算法仍有助对MDS的理解。
关于经典标度法的算法如下：（读书笔记嘛，直接截取书里边的内容）

关于算法的细节可以参考[1]，[3]。
说明：B = HAH，来源于对观察距离δ[i，j]的平方和展开。这个算法就是特征值分解，其思想和主成分分析（PCA）非常相似。在欧氏距离之下，这种算法求得的主坐标（principal coordinates）相当于PCA中前t个主成分的得分。

关于维度的选取
一种常用方法是按照矩阵B的特征值的结构来选择维数。（类似于PCA对主成分的选择）把正特征值按从大到小排列，变得稳定即可。当然处于可视化的考虑，一般都是选择维数为2或3。

下面继续美国城市航线的例子，看看经典标度的效果：
a<-cmdscale(city3, k = 2, eig = TRUE)#结果是一个列表
#观察其特征值（矩阵B的特征值）：
a$eig
 [1] 9.509732e+06 1.719867e+06 4.151052e+05 2.950924e+03 1.560345e+03 -1.164153e-10 -2.628176e+03 -3.497864e+03 -2.268600e+04 -5.233140e+05
#绘原始距离和估计距离的散点图
city2d<-data.frame(a［1］)#把前两个主坐标取出
city2d11<-dist(city2d, method = "euclidean",diag=FALSE)#生成距离矩阵
city2dv1<-as.vector(city2d11)
city2<-as.dist(city1, diag = FALSE)
city2v<-as.vector(city2)
c<-data.frame(city2v,city2dv1)
library(ggplot2)
p<-ggplot(data=c,aes(x=city2v,y=city2dv1))+
geom_point(pch=19,colour='red',size=2,alpha=.6)+
stat_smooth(method=lm)
update_labels( p,list(x="Observed Distance",y="Estimated Distance"))
#绘残差图
res<-c[,1]-c[,2]
d=data.frame(res,x=c(1:45))
ggplot(data=d,aes(y=res,x=x))+
geom_point(pch=19,colour='red',size=2,alpha=.6)+
geom_text(data=d,colour="black",aes(x=x,y=res+10,label=x,fontface=1),size=3)

得到原始距离和由前2个主坐标得到的拟合距离的散点图和残差图,拟合的效果还可以。明显的过拟合和拟合不足的点基本都和丹佛有关，哪位熟悉美国航线的给解释一下？

3.距离标度(Distance Scaling)

距离标度中,降维之后的距离d[i,j]≈f(δ[i，j]),其中f是一个单调增函数。
距离标度可分为度量距离标度（metric distance scaling）：相异度是定量的（比例或区间）和非度量距离标度（nonmetric distance scaling）：相异度是定性的（比如顺序）。
        
（1）度量距离标度

度量距离标度

在度量距离标度中，单调函数f是一个参数化的函数，比如有这样的形式：
f(δ[i，j])=a+bδ[i，j],(a,b是参数)。对这样的形式，很自然的会采用最小二乘法使用{d[i,j]}来对f(δ[i，j])进行拟合。这就是度量最小二乘标度（Metric Least-Squares Scaling）。通常使用的是加权的最小二乘的形式，寻找t(
 (即MDS的一个解)，使用加权的最小二乘回归作为拟合优度的准则：

，
其平方根

称为压力函数（stress function）。其中W=（w（i，j））为权重矩阵。
极小化压力函数，可得到MDS的优化解。这样就可以保证尽可能的保持每对观测之间的距离。求解一般采用梯度下降算法。

经典标度和最小二乘标度的关系
当采用欧氏距离且f为恒等函数时，经典标度可视为最小二乘标度的一个特殊情况。
从最小二乘的观点，经典标度事实上采用了中心化的内积（参见上一小节经典标度算法的第一步），从而存在一个可以用特征向量表示的显示解。

Sammon映射

最小二乘标度的一个变种是Sammon映射(Sammon, 1969)（文献在这里可以找到http://en.wikipedia.org/wiki/Sammon_mapping），其中权重取为

且单调函数f为恒等函数。Sammon映射可以保持较小的距离δ[i，j]，在拟合中给以较大的权重。这样在识别数据的聚类时有更好的效果。求解Sammon mapping采用迭代数值方法（拟牛顿迭代，pseudo-Newton iterative procedure），通常把经典标度解作为初始迭代值。[6]

度量标度的例：劳埃德银行的雇员（Lloyds Bank Employees）数据

这个例子最早来自[5],该文的作者研究的是类似劳埃德银行这样的英格兰银行是如何从1900年代早期的稳定的员工安排变化成如今高度动态的职业类型。

本书提供的这个数据集可由http://astro.temple.edu/~alan/MMST/datasets.html下载。因为原数据集的记录庞大结构复杂，本书选取的是1905-1909以及1925-1929两个时间段入职的劳埃德员工的职业记录。关于数据结构，在下载页有文件说明，把这个数据集转化为距离矩阵采用的是文献[5]中的匹配方法。关于这个方法，在本章最后我专文讨论。在这里先用本书作者给出的距离矩阵samp05d（1905-1909时间段）看一下度量标度的效果。

s05d<-read.csv('samp05d.csv',head=TRUE)
s05d1<-s05d[,-1]
colnames(s05d1)<-c(1:80)
rownames(s05d1)<-colnames(s05d1)
s05dmds<-cmdscale(s05d1,k=2)#求经典标度

下面用MASS包的sammon函数求sammon映射
s05d1<-as.matrix(s05d1, diag = TRUE, upper = TRUE)#sammon要求对称阵或dist class
for(i in(1:80))
{for(j in(1:80))
{if((i!=j)&(s05d1[i,j]==0))
s05d1[i,j]=0.00000001}}#sammon要求输入矩阵非对角线不能为0
library(MASS)
s05dsa<-sammon(s05d1,y=s05dmds,k = 2, niter = 100, trace = TRUE,magic = 0.2, tol = 1e-4) #求sammon映射，使用默认参数

画图比较结果：
s05dsad<-s05dsa［1］#双层中括号
par(mfrow=c(1,2))
plot(s05dmds[,1],s05dmds[,2],type='p',col='blue',pch=19) #绘图text(s05dmds[,1],s05dmds[,2]+.1,labels(colnames(s05d1)),cex=.7) plot(s05dsad[,1],s05dsad[,2],type='p',col='blue',pch=19)
text(s05dsad[,1],s05dsad[,2]+.1,labels(colnames(s05d1)),cex=.7)

从结果看，比对雇员数据，这些银行雇员可分为两个集群右图（sammon映射的结果）数据的结构情况看得更明显，其中29，30，57这三个点是比较明显的离群点，对比数据，他们分别为1587, 1590, 3240号观测，其中1587号（29）在银行工作了59年（最长），而1590（30），3240（57）号只在银行工作了2年（最短）。

图有点小，各自来一张。

（2）非度量距离标度

Shepard diagram

当相异度是定性的次序量,采用非度量标度。此时，就是寻找一个t维空间上的结构
Y={Y[1]，…,Y[n]},用它们的内点距离的大小次序来拟合原始相异度的大小次序。
设对称的相异度矩阵为（δ[i，j]），（其主对角线为0），把其中的相异度按从小到大的顺序排列：
δ[i1，j1]<…<δ[im，jm], (m=n(n-1)/2)
对给定的t，希望t维空间上点的结构Y的内点距离（不一定是欧氏的）能保持原r维空间的相异次序。即d[i1，j1]<…

但事实上，这一点不是总可以做到的。因此，需要对距离d进行修正，变成距离d^,使d[i,j]≈d^[i，j],使d^[ik，jk],k=1，…，m单调不减。这一种近似的替代，称为“不一致”（disparities）。
因此，事实上，我们需要作的是找到一个保持单调性的任意函数f，f满足：
当δ[i，j]< δ[k，l] 时， f(δ[i，j])小于等于f(δ[k，l])
使得d^[i，j]=f(δ[i，j])。
这些disparities叠加在针对次序的结构距离图上是，连起来形成一条曲线，称为Shepard diagram[7]。

进行计算的方法常见的是保序回归（isotonic regression，Kruskal，1964[8，9]）。
关于保序回归，在此不赘述，请参考http://en.wikipedia.org/wiki/Isotonic_regression。

在R的stats包有函数isoreg，可以计算保序回归，看一个简单的例子：
a<-c(2.3,2.7,8.1,5.7,6.2,8.1,8.6,7.7,6.8,9.3,10.5,9.8,10,12.6,12.8)
a.ir<-isoreg(a)
plot(a.ir,col='blue')

图上蓝色的点为初始的向量a,连线即为保序回归拟合的结果。经过保序回归之后，原来不单调的序列a,变成了单调不减的序列：
> a.ir$yf
 [1] 2.300000 2.700000 6.666667 6.666667 6.666667 7.800000 7.800000 7.800000 7.800000 9.300000 10.100000 10.100000 10.100000 12.600000
[15] 12.800000


非度量标度的求解

Shepard diagram的残差平方和：

因残差平方和不能对结构{Y}的尺度伸缩保持不变，因此需要对其进行修正，采用加权形式：

权重常取为：

即有

这称为Kruskal压力公式。也即非度量标度的拟合优度的损失函数。
压力函数还有其他的形式。

利用这个损失函数，一般采用基于梯度的优化方法求非度量标度的解。
算法如下：

但需要注意的是，基于梯度的方法并不一定得到全局的最优解。因此，实践上应尝试不同的初始迭代结构，以检验算法的收敛性质。

关于MDS解的评价

Kruskal给出了关于评价MDS解的一个经验方法：按照stress的值来进行评价：

但是，这个方法并非适用所有的问题。

关于维数的选择

MDS本质上是个维数缩减方法，对于变化的t,t越大，stress越小。
选择维数的一个方法是碎石图（scree Plot）,类似于PCA或者CA选择主成分（因子）的方法，把对应每个t值的最小stress画出来，当stress值变化变平缓的时候，可视为合适的t值。
当然，一般出于可视化的考虑，t常取为2或3。

R的MASS包有函数isoMDS，采用上面所述的Kruskal方法求解非度量标度。
例：先看一个简单的例子，前面我们用过的美国城市航线的例子。
非度量标度当然也可用于度量标度的情形，只要把实体的连续的距离看作表示差异的次序值。
整理数据的过程略，得到dist类city2

library(MASS)
city.nmds<-isoMDS(city2, y = cmdscale(city2, 2), k = 2, maxit = 50, trace = TRUE,
       tol = 1e-3, p = 2)#采用默认参数

和isoMDS相伴随的还有一个表示保序回归拟合的函数Shepard

Shepard(city2, x=city.nmds$points, p = 2)#需要输入最后的迭代解