使用R的统计学习：算法和实践（一）：MDS（4）

###因为表示数学符号不方便，正文中一般用[ ]表示元素下标，用( )表示矩阵或函数，请大家
####按文意甄别

（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[im，jm]。

但事实上，这一点不是总可以做到的。因此，需要对距离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)#需要输入最后的迭代解

参考：
[7]Shepard, R.N. (1962). The analysis of proximities: multidimensional
scaling with an unknown distance function I. Psychometrika, 27,
125–140; II Psychometrika, 27, 219–246.]
[8] Kruskal, J.B. (1964a). Multidimensional scaling by optimizing goodness
of fit to a nonmetric hypothesis, Psychometrika, 29, 1–27.
[9] Kruskal, J.B. (1964b). Nonmetric multidimensional scaling: A numerical
method, Psychometrika, 29, 115–129.