http://www.ma-xy.com
目录
第一章 统计基础 1
1.1 引例：身高问题 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 单总体参数估计与检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 正态总体方差已知的单总体均值的估计与检验 . . . . . . . . . . . . . . . . 3
1.2.2 正态总体方差的估计与检验 . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.3 正态总体方差未知的单总体均值的估计与检验 . . . . . . . . . . . . . . . . 9
1.2.4
两总体方差大小的估计与检验
. . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 方差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.1 单因素方差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.2 方差分析的多重比较 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.3.3 方差齐性检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.4 MATLAB 应用实例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 分布检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.4.1 正态分布检验方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.4.2 适用所有分布的检验方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.5 非参数统计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.5.1 单总体推断与检验 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5.2 两总体位置与尺度推断与检验 . . . . . . . . . . . . . . . . . . . . . . . . . 29
1.5.3 非参数方差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.6 相关性分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.6.1 相关性简介 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
1.6.2 连续变量对连续变量的相关性 . . . . . . . . . . . . . . . . . . . . . . . . . 31
1.6.3 分类变量对分类变量的相关性 . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.6.4 Copula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
1
http://www.ma-xy.com
第一章 统计基础
1.1 引例：身高问题
我们想要了解安徽财经大学男生身高的规律。在查阅相关资料后，得知身高服从正态分布
(或者说大部分人是中等身高，只有小部分人的身高偏高或偏低)。设男生身高为 x，x ∼ N (µ, σ
2
)，
现在的问题是：参数 µ(男生的平均身高)，参数 σ
2
(男生身高的方差) 是多少？
当然，我们可以把安财所有男生的身高记录下来，然后求平均以得到 µ(同样的方法可以得
到方差)，这样得到的 µ 虽然准确，但会耗费很大的人力物力。我们自然会想到：能否只测量部
分男生的身高，然后由部分来推断总体。我们可以从安财男生中随机抽取部分样本 (不同的采样
方法会得到不同的均值，这里采用随机采样)，设样本数为 n，抽取的样本为 x
1
, x
2
, . . . , x
n
，且样
本独立。可以想象的是，只要样本数量 n 给出，我们可以抽很多次样本，用 x = {x
1
, x
2
, . . . , x
n
}
表示 x
1
, x
2
, . . . , x
n
的具体某一次抽样，例如：x = {1.75, 1.73, 1.70, 1.72, 1.90, . . . }。现在，我们
要用样本 x
1
, x
2
, . . . , x
n
来估计 (推断) 总体参数 µ, σ。
方法一：首先，我们来介绍第一种估计思想 - 矩估计。我们很自然想到用下面的计算方法来
估计 µ, σ
ˆµ =
1
n
n
i=1
x
i
ˆ
σ
2
=
1
n
n
i=1
(x
i
− ˆµ)
2
我们说 ˆµ 表示 µ 的估计，
ˆ
σ
2
表示 σ
2
的估计，这是一种直接估计法，直接用样本均值估计 (代
替) 总体均值，样本方差估计总体方差。由此推广，我们可以得到总体 k 阶原点矩 A
k
和中心矩
B
k
(k = 1 , 2 . . . ) 的估计
ˆ
A
k
=
1
n
n
i=1
(x
i
)
k
ˆ
B
k
=
1
n − 1
n
i=1
(x
i
− ˆµ)
k
方法二：这里，我们来介绍第二种估计思想 - 最小二乘估计。学过 OLS(最小二乘) 的人会
想到下面这样的估计：N(µ, σ
2
) 是一个函数，含参数 µ, σ。我们通过给出的点 x
1
, x
2
, . . . , x
n
来
进行拟合 f ≡ N ，找到最好的 µ, σ
2
，使得 f (x|µ, σ) 对 x 的拟合效果最好 (即离差平方和最小)。
1
http://www.ma-xy.com
1.1 引例：身高问题 第一章 统计基础
在给出 x
1
, x
2
, . . . , x
n
后，用 m
i
表示落在 x
i
邻域 x
i
±h 的样本个数，p
i
=
m
i
n
表示落在
x
i
邻域 x
i
± h 的频率 (可能/概率)，令 h → 0，p
i
是 x
i
处的真实值，f (x
i
|µ, σ) 就是要拟合 (逼
近/估计)p
i
，如图 (??) 所示
图 1.1: 身高分布拟合示意图
很自然的有 OLS 估计：我们使离差平方和最小，有
min
µ,σ
J(µ, σ) =
n
i=1
(f(x
i
) − p
i
)
2
其中：在 h 给定后，样本给定后，p
i
是一常量，而 f (x|µ, σ) =
1
√
2πσ
2
e
−
(x−µ)
2
2σ
2
。我们可以用最简
单的 (不用推导的：例如 GA,PSO) 方法求解 µ, σ
2
，易解。但是，在 OLS 中，我们可能不能求
其显式估计量表达式。
方法三：下面，我们来介绍第三种估计思想 - 极大似然估计。在已知总体分布 f : x ∼ N (µ, σ
2
)，
但不知分布参数 µ, σ 时，我们从安财男生中开始 (有放回) 抽取样本 x
1
, x
2
, . . . , x
n
，设第 i 个
样本为 x
i
，则 x
i
被抽中的概率为 p(x
i
)。注意，这里的样本 x
i
是身高。如果给出的 µ, σ
2
比较
好，则每次样本出现的概率都很大，比如我们给定了 µ, σ，并且这个恰好就是真实的样本分布，
那么，我们从这个样本分布中采样的话，概率大的样本更容易被挑中。极端地，可能 n 次都挑中
x
i
= 1 .75(如果 1.75 是分布均值的话)。优良的 µ, σ 可以使样本的出现概率 (联合概率) 最大，假
设样本独立同分布，有
max
µ,σ
J(µ, σ) ≜ L(x
1
, x
2
, . . . , x
n
)
= p(x
1
)p(x
2
) . . . p(x
n
)
=
n
i=1
p(x
i
)
=
n
i=1
1
√
2πσ
2
e
−
(x
i
−µ)
2
2σ
2
在 x
1
, x
2
, . . . , x
n
给出后，L 仅是 µ, σ 的函数，我们在参数 µ, σ 的参数空间 Θ 中求最好的
µ, σ 使得目标 L 最大 max L。由于有连乘
，我们将 L 取对数 log L，然后再关于 µ, σ 求导 (求
极值)，有
ˆµ =
1
n
x
i
ˆσ
2
=
1
n
(x
i
− ˆµ)
2
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.2 单总体参数估计与检验
问：1. 同一参数有不同的估计方法，那么那种方法好？即如何评价估计量的质量。2. 估计量 ˆµ 是
样本的函数，每个一组样本 x
1
, x
2
, . . . , x
n
，都会有一个估计值，样本是随机的，所以估计量 ˆµ 也
是随机的，那么随机变量的分布 (统计特性) 如何？
1.2 单总体参数估计与检验
1.2.1 正态总体方差已知的单总体均值的估计与检验
下面，我们来具体的处理一下安财男生身高问题。我们假设身高 x 服从正态分布 N (µ, σ
2
)，
并且假设分布的方差 σ
2
已知。给出具体样本 x
1
, x
2
, . . . , x
n
后，有
ˆµ =
1
n
n
i=1
x
i
前面，我们提到过估计量 ˆµ 是一个随机变量，下面，我们来研究一下 ˆµ 的分布情况。
∵ x
i
∼ N(µ, σ
2
) iid
∴
n
i=1
x
i
∼ N(nµ, nσ
2
) 正态分布可加性
也即
nˆµ ∼ N(nµ, nσ
2
)
于是有
E(nˆµ) = nµ
⇒ E(ˆµ) = µ
V ar(nˆµ) = nσ
2
⇒ n
2
V ar(ˆµ) = nσ
2
⇒ V ar(ˆµ) =
σ
2
n
即 ˆµ = N(µ,
σ
2
N
)，或者归一化写为
ˆµ − µ
σ
2
/n
∼ N(0, 1)
统计量 ˆµ 的分布如图 (1.2) 所示
图
1.2:
估计量分布图
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.2 单总体参数估计与检验 第一章 统计基础
参数 µ 的区间估计：在上面的估计 µ 的过程中，我们对总体均值 µ 给出了点估计，即估计 µ
是一个具体的值，例如：ˆµ = 1.75。我们也可以给出 µ 的一个估计区间，比如我们说安财男生的
平均身高在 [1.70, 1.80] 之间，我们可以用一个区间来估计 µ，我们记 µ 的估计区间为 [ˆµ
L
, ˆµ
K
]。
从 ˆµ 的分布图 (1.2) 中，我们可以看到，ˆµ 包含在估计区间 [ˆµ
L
, ˆµ
K
] 的概率为
P =
ˆµ
K
ˆµ
L
p(ˆµ)dµ
其中：p(ˆµ) 为 µ 的密度函数。概率 P 为多次采样 x 后估计区间 [ˆµ
L
(x), ˆµ
K
(x)] 包含 µ 的概率，
其示意图如图 (1.3) 所示
图 1.3: 区间估计示意图
图 (1.3) 的横轴为总体参数 µ，每一个纵向线段为一次区间估计，当取 n 次区间估计时，这
n 次中会有 m = P n 次包含 µ。
参数 µ 的假设检验：在 σ
2
= 1 已知的情况下，我们假设安财男生身高 x 的平均值为 µ = 0.75。
现在我们来考虑这样的假设是否合理？我们知道，在原假设成立的情况下，将样本 x 带入估计量
ˆµ 当中，有
ˆµ(x) − 0.75
1/n
= 0 .6
其中：x 和 n 皆已知。由
ˆµ−µ
√
σ
2
/n
∼ N(0, 1)，我们可以找到 0.6 对应的概率 f(0.6)，这里的 f 就
是标准正态分布 N (0, 1)。f(0.6) 即为样本出现的概率，即抽一次样本，样本为 x 的概率，如图
(1.4) 所示
图 1.4: 假设检验示意图 1
样本出现的可能性太小，即拒绝 µ = 0.75 假设 (原假设)，因为 µ 有更大的可能取其它值。
为此，我们让 µ 从 0.5 到 1.85 遍历取值，哪一个 µ 值使样本出现的概率最大，我们就选择这个
假设，这也对应了极大似然估计的思想，只有当 ˆµ −µ = 0 时，样本出现的概率才最大。我们可
以在检验统计量分布中设置界限 (阈值)α，如图 (1.5) 所示
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.2 单总体参数估计与检验
图 1.5: 假设检验示意图 2
当
ˆµ(x)−µ
√
1/n
= ∆(其中 µ = 1.75，∆ > Z
1−
α
2
) 时，我们不能接受原假设 µ = 1.75，称在 α 显
著水平下，不能接受原假设。这里的 α 是一个概率，称为显著水平，Z
1−
α
2
是一个样本点 (横坐
标 z(z =
ˆµ−µ
√
1/n
) 上的一个点)，P {Z
α
2
< z < Z
1−
α
2
} = 1 − α。
1.2.2 正态总体方差的估计与检验
我们假设 x ∼ N(µ, σ
2
)，σ
2
未知。有样本 x = {x
1
, x
2
, . . . , x
n
}，下面考虑如何依据样本来
估计总体的分布参数 σ
2
？
1、在总体参数 µ 已知的情况下，同前，有
ˆσ
2
=
1
n
n
i=1
(x
i
− µ)
2
2、如果 µ 未知，我们需要用样本均值 (估计量)ˆµ 来代替 µ，于是有
ˆσ
2
=
1
n − 1
n
i=1
(x
i
− ˆµ)
2
≜ s
2
我们知道 s
2
是一个随机变量
¬
，并且知道 ˆµ ∼ N (µ, σ
2
/n)，下面，我们来求解随机变量 s
2
的分布
。我们将样本均值 ˆµ 记为 ¯x，有
s
2
=
1
n − 1
n
i=1
(x
i
− ¯x)
2
即
(n − 1)s
2
=
n
i=1
(x
i
− ¯x)
2
=
n
i=1
(x
2
i
+ ¯x
2
− 2x
i
¯x)
=
n
i=1
x
2
i
− (
√
n¯x)
2
这里我们简单的设想一下：
n
i=1
x
2
i
∼ χ
2
(n)，(
√
n¯x)
2
∼ χ
2
(1)，那么，s
2
∼ χ
2
(n − 1)。我
们来详细看一下刚才的设想，x
i
∼ N (µ, σ
2
)，¯x ∼ N(µ,
σ
2
n
)，
√
n¯x ∼ N(
√
nµ, σ
2
)，我们需要 x
2
i
的分布。
¬
本书中的 :=、≜、≡ 均表示等价
可以参考《数理统计引论》P4 或者《数理统计学》P36 定理 1.3.3
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.2 单总体参数估计与检验 第一章 统计基础
引理 (卡方分布) 假设 x
1
, x
2
, . . . , x
n
iid
∼ N (0, 1)，则
X =
n
i=1
x
2
i
∼ χ
2
(n)
其中：n 为卡方分布 χ
2
的自由度。
下面，我们来看一下卡方分布 χ
2
(n) 的分布函数和密度函数。设 k(x|n) 为密度函数，K(x|n)
为分布函数，当 x ⩽ 0 时，有 K(x|n) = 0；当 x > 0 时，有
K(x|n) = P (X < x)
= P
n
i=1
x
2
i
< x
= P
x
T
x < x
= ( −2π)
−
n
2
···
B
exp
−
1
2
x
T
x
dx (1.1)
其中：B 为 R
n
的球体 {x : x
T
x < x}，而 dx = dx
1
. . . dx
n
。作为示例，我们先来看一个二维情
况：
x
2
1
+x
2
2
<x
e
−
1
2
(x
2
1
+x
2
2
)
dx
1
dx
2
令
x
1
= r cos θ
x
2
= r sin θ
于是有
2π
0
dθ
√
x
0
re
−
1
2
r
2
dr
将 (1.1) 转化为球面坐标，不难看出 x 被积函数将变为 D(θ
1
, . . . , θ
n−1
)e
−r
2
/2
r
n−1
的形状，
且 ( θ
1
, . . . , θ
n−1
) 的积分范围与 r 无关，于是有
K(x|n) = C
n
√
x
0
e
−r
2
/2
r
n−1
dr
其中：C
n
只与 n 有关。为求 C
n
，由分布 K(x → ∞) = 1，令 x → ∞，有
1 = C
n
∞
0
e
−r
2
/2
r
n−1
dr
= C
n
2
n
2
−1
Γ
n
2
由此得到 C
n
C
n
=
1
2
n
2
−1
Γ
n
2
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.2 单总体参数估计与检验
然后带入 K(x|n) 中，有
K(x|n) =
1
2
n
2
−1
Γ
n
2
√
x
0
e
−
1
2
r
2
r
n
−
1
dr
两边对 x 求导，有
k(x|n) =
1
2
n
2
−1
Γ
n
2
e
−
x
2
x
n
2
−1
卡方分布 χ
2
(x|n) 中，n 为自由度参数，当 n 取不同值时，卡方分布的图像如图 (1.6) 所示
®¯
图 1.6: 卡方分布
推论 (卡方分布的一般形式 1) 在前面的卡方分布中，我们假设 x
1
, x
2
, . . . , x
n
独立同分布于
N(0, 1)，现在，我们来把假设放宽。假设 x
1
, x
2
, . . . , x
n
独立同分布于 N(0, σ
2
)，令 X =
n
i=1
x
2
i
，
则
X/σ
2
∼ χ
2
(n)
推论 (卡方分布的一般形式 2) 假设 x
1
, x
2
, . . . , x
n
独立同分布于 N(µ, σ
2
)，令 X =
n
i=1
x
2
i
，
则
X/σ
2
∼ χ
2
(n, r)
其中：r = n
µ
2
2σ
2
为非中心参数。
推论 (卡方分布的一般形式 3) 假设 x
1
, x
2
, . . . , x
n
独立，x
i
∼ N(µ
i
, σ
2
), i = 1, 2, . . . , n，令
X =
n
i=1
x
2
i
，则
X/σ
2
∼ χ
2
(n, δ)
其中：
δ =
n
i=1
µ
2
i
σ
1/2
为非中心参数。
®
https://en.wikipedia.org/wiki/Chi-squared_distribution
¯
图中的 k 即为 n
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.2 单总体参数估计与检验 第一章 统计基础
下面，我们继续来求解 s
2
的分布。根据正态分布标准化，我们有
x
i
∼ N(µ, σ
2
)
x
i
− µ
σ
∼ N(0, 1)
√
n¯x ∼ N (
√
nµ, σ
2
)
√
n¯x −
√
nµ
σ
∼ N(0, 1)
由卡方分布，我们有
n
i=1
(x
i
− µ)
2
σ
2
∼ χ(n)
[
√
n(¯x − µ)]
2
σ
2
∼ χ(1)
并且
n
i=1
(x
i
− µ)
2
− [
√
n(¯x − µ)]
2
=
x
2
i
+ nµ
2
− 2µ
x
i
− n(¯x − µ)
2
=
x
2
i
− µ
2
− (
√
n¯x)
2
+ µ
2
=
x
2
i
− (
√
n¯x)
2
= : (n − 1)s
2
于是
(n − 1)s
2
σ
2
=
n
i=1
(x
i
− µ)
2
σ
2
−
[
√
n(¯x − µ)]
2
σ
2
∼ χ
2
(n) − χ
2
(1)
引理 (卡方分布的可加性) 设 x
1
∼ χ
2
(m)，x
2
∼ χ
2
(n)，则 x
1
+ x
2
∼ χ
2
(m + n)。
由上面的卡方分布的可加性，我们有
(n − 1)s
2
σ
2
∼ χ
2
(n − 1)
至此，我们得到了 s
2
的分布。利用
(n−1)s
2
σ
2
∼ χ
2
(n − 1)，我们可以对总体方差 σ
2
做假设检
验。我们 (原) 假设 H0 : σ
2
= 1 ，在 H0 成立的条件下，
T (x) =
(n − 1)s
2
(x)
σ
2
= 1
∼ χ
2
(n − 1)
将具体样本 x
1
, x
2
, . . . , x
n
带入 T ，可以得到 T 的一个具体值，例如 T
1
= 0.569。我们在卡方分
布 χ
2
(n − 1) 下找 T
1
= 0 .569 对应的概率值，例如 P
1
= 0 .23，并计算 p 值
p 值 = 1 − P (T < 0.569)
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.2 单总体参数估计与检验
如果 p 值较小，则说明样本出现的概率较小，不能接受原假设 H0。我们还可以设置规则 α，当
T
1
> χ
2
1−α
时，不能接受原假设 H0，这里的 χ
2
1−α
是 χ
2
的 1 − α 分位数。统计上称规则 α 为
显著水平，一般取值 0.1、0.05 和 0.01。注意，这里的 p 值和规则 α 都是概率。
卡方分布的性质
卡方分布具有如下性质：
1). 均值方差。若 X ∼ χ
2
(n)，则 E(X) = n，V ar(X) = 2n。
2). 若 Y
1
, Y
2
, . . . , Y
m
独立，Y
i
∼ χ
2
(n
i
, δ
i
)，则
Y =
Y
i
∼ χ
2
(n, δ)
其中：n =
n
i
，δ
2
=
δ
2
i
。
3). 若 X
1
, X
2
, . . . , X
n
∼ χ
2
(n)，则当 n → ∞ 时，有
X
n
− n
2n
L
−→ N (0, 1)
且
√
2X
n
−
√
2n
L
−→ N (0, 1)。其中：
L
−→ 表示依分布收敛。
4). 若 X
i
∼ χ
2
(n, δ
i
)，i = 1, 2，而 δ
1
> δ
2
，则 X
1
> X
2
。
1.2.3 正态总体方差未知的单总体均值的估计与检验
在前面的单总体均值估计部分 (1.2.1)，我们假设总体方差 σ
2
已知，但是在实际问题中，我
们往往只能假设安财男生身高 x ∼ N (µ, σ
2
)，并不知道总体方差 σ
2
是多少，我们只能利用样本
来推断 (估计) 总体。总体均值 µ 的估计仍然可以使用
ˆµ =
1
n
n
i=1
x
i
∼ N(µ, σ
2
/n)
但是总体均值的检验则不能使用下面的检验统计量
ˆµ − µ
σ
2
/n
∼ N(0, 1)
因为这里的总体方差 σ
2
不知道。很自然的，我们想到用样本方差 s
2
代替总体方差，那么新的检
验统计量为
ˆµ − µ
s
2
/n
接下来要讨论的问题是：新的统计量的分布是什么？我们来简单看一下随机变量
ˆµ−µ
√
s
2
/n
，分
式上面的 ˆµ − µ 是正态分布，分式下面的
s
2
/n 是近似的卡方分布 (开平方)，下面，我们来详
细推导这个随机变量的分布。
设
T =
ˆµ − µ
s
2
/n
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.2 单总体参数估计与检验 第一章 统计基础
已知
ˆµ−µ
√
σ
2
/n
∼ N(0, 1)，
(n−1)s
2
σ
2
∼ χ
2
(n − 1)，于是有
T =
ˆµ − µ
σ
2
/n
σ
2
/n
s
2
/n
=
ˆµ − µ
σ
2
/n
σ
2
s
2
=
ˆµ−µ
√
σ
2
/n
s
2
σ
2
我们下面要求
s
2
σ
2
。令 Z =
s
2
σ
2
，设其密度函数为 f
Z
(z)，分布函数为 F
Z
(z)，当 z ⩽ 0 时，
f
Z
(z) = 0，当 z > 0 时，有
F
Z
(z) = P {Z ⩽ z}
= P
s
2
σ
2
⩽ z
= P
s
2
σ
2
⩽ z
2
= P
ns
2
σ
2
< nz
2
= F
Y
(nz
2
)
上式中，我们令 Y =
ns
2
σ
2
∼ χ
2
(n)。于是有 Z 的密度函数为
f
Z
(z) = F
′
Z
(z)
= F
Y
(nz
2
)
′
z
= f
Y
(nz
2
)(2nz)
=
1
2
n
2
−1
Γ
n
2
n
n
2
z
n−1
e
−
nz
2
2
z > 0
令 X =
ˆµ−µ
√
σ
2
/n
，则 T = X/Z。现在，我们知道了 X, Z 的分布函数，要求二者商的分布函
数，引入随机变量商的分布：
引理 (随机变量商的分布) 设 x
1
∼ g(x)，x
2
∼ h(x)，且二者之间相互独立。令 z =
x
1
x
2
，则
随机变量 z 的密度函数为
f(x) =
∞
0
tg(xt)h(t)dt
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.2 单总体参数估计与检验
由上面的随机变量商的分布，我们有 T = X/Z 的分布
f
T
(t) =
∞
−∞
|z|f
Z
(z)f
X
(zt)dz
=
∞
0
z
1
√
2π
e
−
z
2
t
2
2
1
2
n
2
−1
Γ
n
2
n
n
2
z
n−1
e
−
nz
2
2
dz
=
n
n
2
√
2π2
n−1
2
Γ(
n
2
)
∞
0
z
n
e
−
z
2
2
(n+t
2
)
dz
令u=
n+t
2
2
z
2
=========
1
√
nπΓ
n
2
1 +
t
2
n
n+1
2
∞
0
u
n+1
2
−1
e
−u
du
=
Γ(
n+1
2
)
√
nπΓ(
n
2
)
1 +
t
2
n
−
n+1
2
上式中用到了 Γ 函数的表达式
Γ(α) =
∞
0
x
α−1
e
−x
dx
也是因为要拼凑 Γ 函数，所以有变换 u =
n+t
2
2
z
2
。我们用 t(n) 来表示 t 分布的密度函数 f
T
(t)，
其中，n 为自由度参数，当 n 取不同值时，t 分布密度函数图像如 T 分布密度函数图 (1.7) 所示
图 1.7: Γ 分布密度函数图
我们利用 T = X/Z ∼ t(n −1) 可以进行方差未知时总体均值 µ 的假设检验，这里不再叙述。
t 分布的性质
t 分布具有以下性质：
1). 自由度为 1 的 t 分布为柯西分布，它的期望不存在。
2). n > 1 时，t 分布的数学期望存在，并且为 0。
3). n > 2 时，t 分布的方差存在，并且为 n/(n − 2)。
4). 自由度 n 越大，t(n) 分布越接近标准正态分布。当 n → ∞ 时，t(n) 分布的极限分布为
标准正态分布。一般认为，当 n > 30 时，t(n) 可以用标准正态分布近似。
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.2 单总体参数估计与检验 第一章 统计基础
1.2.4 两总体方差大小的估计与检验
传说中，统计有三大分布：卡方分布 χ
2
、t 分布和 F 分布。前面，我们给出了 χ
2
分布和 t
分布的简单形式，下面，我们引出 F 分布。
蚌埠有 5 所大学，我们选择蚌埠医学院和安财来做研究。我们想知道 2 所高校男生身高离
散程度 (方差) 的差异，直接的说，两所高校男生身高的方差是否相同，谁大谁小？设安财男生
身高总体为 x
1
∼ N(µ
1
, σ
2
1
)，蚌医男生身高总体为 x
2
∼ N(µ
2
, σ
2
2
)，不论总体均值是否知道，在
σ
2
1
, σ
2
2
未知的情况下，问 σ
2
1
, σ
2
2
的关系如何？有样本 x
11
, x
12
, . . . , x
1n
; x
21
, x
22
, . . . , x
2m
，如果仅
考虑参数估计问题，我们用样本把 σ
2
1
, σ
2
2
估计出来 ˆσ
2
1
, ˆσ
2
2
，然后比较大小即可。如果是假设检验
问题，我们假设二者的方差相等 H0 : σ
2
1
= σ
2
2
，我们可以尝试构造下面两种检验统计量
T
1
=
ˆσ
1
2
ˆσ
2
2
T
2
= ˆσ
1
2
− ˆσ
2
2
如果假设 H0 成立，那么 T
1
取值应该在 1 左右，T
2
取值应该在 0 左右。但是，对于假设检
验，我们不仅仅要看统计量的取值，更要求解统计量的分布，对于第二个统计量 T
2
，我们知道
nˆσ
2
1
σ
2
1
∼ χ
2
(n)，
mˆσ
2
2
σ
2
2
∼ χ
2
(m)，但是由于方差 σ
2
未知，我们不能求出具体的统计量值。对于第一
个统计量 T
1
T
1
=
ˆσ
1
2
ˆσ
2
2
=
nˆσ
2
1
σ
2
1
σ
2
1
n
m
ˆ
σ
2
2
σ
2
2
σ
2
2
m
令 X =
nˆσ
2
1
σ
2
1
∼ χ
2
(n)，Y =
mˆσ
2
2
σ
2
2
∼ χ
2
(m)，则
T
1
=
X
Y
m
n
可以看出，统计量 T
1
是随机变量商的形式，我们可以用随机变量商的分布来求解。有两种
思路：1、将
X
n
和
Y
m
视为两个随机变量，先求各自的分布，然后再求二者商的分布，由于 X, Y
的相似性，我们只需要求解一个就可以了。2、将
X
Y
和
m
n
视为两部分，先求
X
Y
的商形式的分布，
再求 T
2
。
对于第一种方法：我们知道 X ∼ χ
2
(n)，Y ∼ χ
2
(m)，令 X
1
=
X
n
，其分布函数为 F
X
1
(x
1
)
其密度函数为 f
X
1
(x
1
)。
F
X
1
(x
1
) = P {X
1
< x
1
}
= P
X
n
< x
1
= P {X < nx
1
}
= F
X
(nx
1
)
=
nx
1
0
χ
2
dx
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.2 单总体参数估计与检验
其中：χ
2
为卡方分布密度函数。X
1
密度函数为
f
X
1
(x
1
) = F
X
1
(x
1
)
′
x
1
= F
X
(nx
1
)
′
x
1
= f
X
(nx
1
)n
上式最后的求导是积分限含参变量的求导，可以继续求解，这里不再叙述。
对于第二种方法：首先，我们导出 Z = X
1
/X
2
的密度函数。记 p
1
(x), p
2
(x) 为 X ∼ χ
2
(n), Y ∼
χ
2
(m) 的密度函数，根据独立随机变量上的分布的引理，我们有
f
Z
(z) =
∞
0
x
2
p
1
(zx
2
)p
2
(x
2
)dx
2
=
z
n
2
−1
Γ
n
2
Γ
m
2
2
n+m
2
∞
0
x
n+m
2
−1
2
e
−
x
2
2
(1−z)
dx
2
令u=
x
2
2
(1+z)
==========
z
n
2
−1
(1 + z)
−
n+m
2
Γ
n
2
Γ
m
2
∞
0
u
m+n
2
−1
e
−u
du
=
Γ
n+m
2
Γ
n
2
Γ
m
2
z
n
2
−1
(1 + z)
−
n+m
2
z > 0
然后，我们导出 T
1
= Z
m
n
的密度函数。对于 y > 0，有分布函数
F
T
1
(y) = P {T
1
< y}
= P
Z
m
n
< y
= P
Z <
n
m
y
= F
Z
n
m
y
分布函数对 y 求导，有
f
T
1
(y) = F
T
1
(y)
′
y
= F
Z
n
m
y
′
y
=
f
Z
n
m
y
n
m
=
Γ
n+m
2
Γ
n
2
Γ
m
2
n
m
y
n
2
−1
1 +
n
m
y
−
n+m
2
n
m
=
Γ
n+m
2
Γ
n
2
Γ
m
2
n
m
n
2
y
n
2
−1
1 +
n
m
y
−
n+m
2
我们记 F 分布为 F (n, m)，其中，n, m 是分布的自由度。
F 分布的性质
F 分布具有如下性质：
1). m > 2 时，F 分布的数学期望存在，且为 m/(m − 2)。
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.3 方差分析 第一章 统计基础
2). m > 4 时，F 分布的方差存在，且为
2m
2
(n+m−2)
n(n−2)
2
(m−4)
。
3). 若 F ∼ F (n, m)，则
1
F
∼ F (m, n)。
4). 若 t ∼ t(n)，则 t
2
∼ F (1, n)。
1.3 方差分析
1.3.1 单因素方差分析
方差分析用于研究分类变量对连续变量的影响。前面，在研究安财男生身高时，我们是从学
校男生中进行随机抽样的。现在，我们对随机抽取的样本进行标记，标记每个男生所在的学院，
假设共有 4 个学院，则对每个样本 (男生) 而言，除了有“身高”连续变量外，还有“学院”分类
变量，其数据结构如图 (1.8) 所示。当然，我们也可以直接从 4 个学院中各取相同数量的男生。
图 1.8: 安财 4 个学院男生身高数据表
现在要研究不同学院男生身高是否有差异？也即学院变量是否影响身高变量。设分类变量
(学院) 为 A(factor A)，k 为学院种类数，有 A
1
, A
2
, A
3
, A
4
；设身高变量为 x，学院 A
i
的样
本为 x
i
，样本数为 n
i
，学院 A
i
的第 j 个样本可以表示为 x
ij
，我们假设 x
i
∼ N (µ
i
, σ
2
i
)，并
假设 σ
1
= σ
2
= σ
3
= σ
4
。现在，我们可以进行假设检验，我们 (原) 假设：不同学院男生身
高均值相同，或者说所有 x
ij
均来自同一正态总体，即 H0 : µ
1
= µ
2
= µ
3
= µ
4
。逆假设为：
H1 : ∃i, j ∈ k, µ
i
= µ
j
。
下面，我们要构造检验统计量。如果 H0 成立 (为真)，各学院男生身高均值相等，那么样本
均值也应该接近 ¯x
1
≈ ¯x
2
≈ ¯x
3
≈ ¯x
4
≈
¯
¯x，这里的
¯
¯x 是所有男生身高的均值。依据各组身高方差
相等的假设，并且在各组样本量差别不大的情况下，我们有
(x
1j
− ¯x
1
)
2
≈
(x
2j
− ¯x
2
)
2
≈
(x
3j
− ¯x
3
)
2
≈
(x
4j
− ¯x
4
)
2
我们称
SSE =
k
i=1
n
i
j=1
(x
ij
− ¯x
i
)
2
为组内离差平方和，SSE 为 sum of squres of error 的缩写。称
SST =
i
j
(x
ij
−
¯
¯x)
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.3 方差分析
为总的离差平方和，这里的 SST 是 sum of squres of total 的缩写。二者相减，有
SST − SSE = SSA
=
k
i=1
n
i
(¯x
i
−
¯
¯x)
2
称 SSA 为组间离差平方和，SSA 是 sum of squres of factor A 的缩写。可以看出，在原假设 H0
成立的情况下，有 SST = SSE，即 SSA = 0。只要 SSA 不等于 0，我们就说分类变量 A 影响
身高 x，不过，我们要讨论这种影响是否显著。根据前面介绍的卡方分布，我们有
(n − 1)SST
σ
2
∼ χ
2
(n − 1)
(k − 1)SSA
σ
2
∼ χ
2
(k − 1)
(n − k)SSA
σ
2
∼ χ
2
(n − k)
但是由于 σ
2
未知，无法计算出具体的检验统计量值，所以我们构造 F 分布
F
1
=
SSE
SST
F
2
=
SSA/(k − 1)
SSE/(n − k)
∼ F (k − 1, n − k)
可以尝试证明统计量 F
2
的分布 F (k −1, n −k)。下面，我们给出分类变量 A 对 x 的影响强
度的度量
R
2
=
SSA
SST
∈ [0, 1]
当 R
2
= 1 时，表明 A 强烈影响 x；当 R
2
= 0 ，表明 A 不影响 x。
1.3.2 方差分析的多重比较
在前面的单因素方差分析当中，我们检验了各个学院男生身高均值是否相等，其原假设为
H0 : µ
1
= µ
2
= µ
3
= µ
4
。如果检验的结果没有接受原假设，那么，接下来我们就要考虑：是哪
一个学院男生身高平均值特别，还是所有学院的身高均值都不相同？对于这个问题，我们有如下
思路：
我们仍然从假设检验做起，现在的目标是找到与众不同的 µ
i
。共有 4 个学院，我们可以两
个学院两个学院的进行检验，检验 12, 13, 14, 23, 24, 34 六组的两两总体的均值是否相同。对于第
i 组和第 j 组
Step1: 原假设 H0 : µ
i
= µ
j
，备择假设 H1 : µ
i
= µ
j
。
Step2: 构建统计量。我们知道
ˆµ
i
∼ N(µ
i
, σ
2
/n
i
)
ˆµ
j
∼ N(µ
j
, σ
2
/n
j
)
ˆµ
i
− ˆµ
j
∼ N
µ
i
− µ
j
,
(n
i
+ n
j
)σ
2
n
i
n
j
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.3 方差分析 第一章 统计基础
我们令
T =
(ˆµ
i
− ˆµ
j
) − (µ
i
− µ
j
)
(n
i
+n
j
)σ
2
n
i
n
j
∼ N(0, 1)
在各组方差相等但未知的情况下 (注意：这是前面单因素方差分析的假设)，上式的 σ
2
未知，用
样本方差代替，有 T ∼ t(n
i
+ n
j
− 1)。
Step3
:
在
H
0
为真时，将得到的样本带入
T
中，即可进行检验。
在所有两两分组检验结束后，会得到一个表格 (1.1)
表 1.1: 方差分析多重比较结果表
组均值 µ
1
µ
2
µ
3
µ
4
µ
1
1 1 1 0
µ
2
1 0 0
µ
3
1 0
µ
4
1
表 (1.1) 中的 1 表示不能拒绝原假设 H0，0 表示拒绝原假设。从表中可得到的结论为
(1, 4), (2, 3), (2, 4), (3, 4) 组的均值有显著差异。存疑：1 与 2 同，1 与 3 同，2 与 3 不同？
常用的多重比较有两种方法：一个是费希尔提出的最小显著差异方法 (LSD)，一个是 N-K
检验。LSD 的基本步骤和上面的步骤相似，不过其检验统计量为
T =
ˆµ
i
− ˆµ
j
MSE
1
n
i
+
1
n
j
∼ t(n − k )
其中：MSE =
SSE
n−k
=
i
j
(x
ij
− ¯x
i
)
2
。N-K 检验是非参数统计中的方法，留在后面介绍。
1.3.3 方差齐性检验
常用的方差齐性检验有 Bartlett 检验 (1937) 和 Levene 检验 (1960)。在前面 (1.2.4) 的安财
和蚌医男生身高方差比较当中，我们曾讨论过方差是否相等这个问题，并由此导出了 F 分布。但
是，这仅是两个方差的比较，如果在多个总体的情况下，检验统计量就变得无能为力了。在单因
素方差分析中，我们曾假设各个分组的方差相等 σ
2
1
= σ
2
2
= σ
2
3
= σ
2
4
= σ
2
，那如何检验各组的
方差是否相等呢？下面，我们来介绍用于方差是否相等的检验-Bartlett 检验，也就是 SPSS 中的
ANOVA 分析的方差齐性检验的首选检验方法。
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.3 方差分析
仍然考虑安财 4 个学院的男生身高数据，见图 (1.8)，定义两个量 MSE 和 GM SE。
MSE =
SSE
n − k
=
k
i=1
n
i
j=1
(x
ij
− ¯x
i
)
2
n − k
=
k
i=1
Q
i
n − k
=
k
i=1
s
2
i
(n
i
− 1)
n − k
=
k
i=1
n
i
− 1
n − k
s
2
i
其中：Q
i
是第 i 组的离差平方和，s
2
i
是第 i 组的样本方差，有 s
2
i
=
Q
i
n
i
−1
。从 M SE 的计算公
式可以看出，MSE 是各总体方差 s
2
i
的加权算术平均。定义 GM SE 为
GMSE =
n−k
k
i=1
(s
2
i
)
(n
i
−1)
从 GM SE 的计算公式来看，GMSE 是各总体方差 s
2
i
的几何平均。对于算术平均数和几何平均
数，我们有
x
1
+ x
2
+ ··· + x
n
n
⩾
n
√
x
1
x
2
. . . x
n
即算术平均数大于几何平均数，只有在 x
1
= x
2
= x
n
时，等号才成立。根据算术平均数和几何
平均数的这个性质，我们有
GMSE ⩽ MSE
只有当各总体方差 s
2
i
相等时，等号才成立，各组方差差异越大，GMSE 越小于 MSE。此外，
我们还可以尝试用熵来度量方差的离散程度。我们可以构建检验统计量
T =
MSE
GMSE
将
T
和
1
比较大小。同样可以构建
T = ln
MSE
GMSE
将 T 和 0 比较大小。
Bartlett 检验的原假设为 H0 : σ
2
1
= σ
2
2
= σ
2
3
= σ
2
4
= σ
2
，备择假设为 H1 : ∃i, j ∈ K, σ
2
i
= σ
2
j
。
检验统计量为
T ≜ B =
n − k
C
(ln MSE − ln GM SE)
·
∼ χ
2
(k − 1)
其中：C = 1 +
1
3(k+1)
k
i=1
1
n
i
−1
−
1
n−k
为校正数。
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.3 方差分析 第一章 统计基础
修正 Bartlett 检验：针对样本量小于 5 不能用 Bartlett 检验的缺点，Box 修正了 Bartlett
检验，修正后的检验统计量为
B
′
=
f
1
BC
f
1
(A − B − C)
·
∼ F (f
1
, f
2
)
其中：
f
1
=
k
−
1
，
f
2
=
k+1
(C+1)
2
，A =
f
1
2−C+2/f
2
，
B, C
同前面介绍。
Levene
检验是对模型残差进
行组间齐性检验 (残差方差齐性检验)。
1.3.4 MATLAB 应用实例
1). 单因素一元方差分析的 matlab 命令为
[p,table,stats] =anova1(x,group,displayopt)
其中：p 为 p 值；table 为元胞数组数据类型的方差分析表，包含组内方差、组间方差和总方差
等以及他们对应的自由度、均方等等；stats 为结构体，可作为后续的多重比较的输入参数；x 可
以为一列数据或者多列分组数据 (矩阵)，和前面介绍的一致；group 为分组，当 x 为一列数据时，
group 的长度和 x 相等，用来指定 x 中各元素所在的组，当 x 为矩阵时，可以忽略 group 这个
参数，或者用它来指定组名。displayopt 设置为’o’ 时，不显示方差分析表 table 和箱线图。
单因素方差分析的多重比较命令为
[c,m,h,gnames] = multcompare(stats,param,. . . )
其中：c 两两比较结果矩阵，是一个多行 5 列的矩阵，每一行对应某两两比较结果，例如某一行
为：2, 5, 1.9442, 8.2206, 14.4971，表示第 2 组和第 5 组进行比较，两组的均值差为 8.2206，均值
差的 95% 置信区间为 [1.9442, 14.4971]，这个区间不包含 0，说明在 0.05 显著水平下，两组之间
的均值是显著差异的；m 为多行两列矩阵，第 1 列为每组均值，第 2 列为相应的标准误差；h 为
交互式多重比较图形句柄；games 是一个元胞数组，每一行对应一个组名；stats 为 anova 返回
的结构体；param 为参数，其可选参数和参数值如表 (1.2) 所示。
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.3 方差分析
表 1.2: 多重比较参数表
参数名 参数值 说明
’alpha’ (0,1) 内的标量 用来指出输出矩阵 c 中的置信区间和
交互式图形上的置信区间的置信水平：
100(1-alpha)%.’alpha’ 的默认值是 0.05
’display’ ’on’, ’o’ 用 来 指 定 是 否 显 示 交 互 式图 形，若
为’on’(默 认 情 况)， 则 显 示 图 形； 若
为’o’，则不显示图形
’ctype’ ’hsd’ 或’tukey-kramer’, ’lsd’
或’bonferroni’, ’dunn-sidak’, ’schee’
用来指定多重比较中临界值的类型，实
际上就是指定多重比较的方法。可用
的方法有 Tukey-Kramer 法 (默认情
形)、最小显著差数法 (LSD 法)、Bon-
ferroni t 检验法、Dunn-Sidak 检验法
和 Schee 法
’dimension’ 正整数向量 对于多因素方差分析的比较检验，用来
指定要比较的因素 (分组变量) 的序号，
默认值为 1，表示第 1 个分组变量。仅
适用于 stats 是 anovan 函数的输出的
情形
’estimate’ 依赖于生成结构体变量 stats 所用的函
数
指定要比较的估计，其可能取值依赖于
生成结构体变量 stats 所用的函数。对
于 anoval、anovan(多因素方差分析)、
friedman(Friedman 秩方 差分析) 和
kruskalwallis(Kruskal-Wallis 单因素方
差分析) 函数，该参数将被忽略；对
于 anova2(双因素方差分析) 函数，’es-
timate’ 参数的可能取值为’column’(默
认) 或’row’，表示对列均值或行均值进
行比较；对于 aoctool(交互式协方差分
析) 函数，’estimate’ 参数的可能取值
为’slope’、’intercept’ 或’pmm’，分别表
示对斜率、截距或总体边缘均值进行比
较。后面会陆续介绍这些函数的用法
2). 双因素一元方差分析的 matlab 命令为
[p,table,stats] = anova2(x,reps,displayopt)
其中：x 的每一列对应因素 A 的一个水平 (因子)，每一行对应因素 B 的一个水平，anova2 检验
http://www.ma-xy.com 19 http://www.ma-xy.com
http://www.ma-xy.com
1.3 方差分析 第一章 统计基础
矩阵 x 的各列是否具有相同的均值，即因素 A 对 x 的影响是否显著，同时还检验 x 的各行是否
具有相同均值；reps 默认为 1，当 reps 的取值大于 1 时，anova2 还检验因素 A 和 B 的交互作
用是否显著；p 为检验的 p 值，如果仅考虑因素 A 和 B，则 p 是一个 2 个元素的行向量，如果
考虑 A 和 B 的交互作用，p 是一个 3 个元素的行向量。
在多重比较中，将参数’estimate’ 的值设置为’column’ 以对因素 A 进行多重比较，设置为’row’
以对因素 B 进行多重比较。
3). 多因素一元方差分析的 matlab 命令为
[p,table,stats,terms] = anovan(y,group,param,. . . )
其中：group 是一个元胞数组，每一个元胞对应一个因素，元胞长度和 y 相同，其元素可以是分
类数据、数值和字符串等等；anovan 支持的参数 param 和参数值 paramvalue 如表 (1.3) 所示，
多因素一元方差分析模型的类型如表 (1.4) 所示
表 1.3: anovan 参数表
参数名 参数值 说明
’alpha’ (0,1) 内的标量 用来指定置信区间的置信水平：100(1-
alpha)%。’alpha’ 的默认值是 0.05
’continuous’ 下标变量 用来指明哪些分组变量被作为连续变
量，而不是离散的分类变量
’display’ ’on’, ’o’ 用来 指定是 否显 示交互 式图 形，若
为’on’(默认 情况)， 则显 示图 形；若
为
’o’
，则不显示
’model’ 所有可能取值如表 (1.4) 所列 用来指定所用模型的类型
’nested’ 由 0 和 1 构成的矩阵 M 指定分组变量之间的嵌套关系。若第 i
个变量嵌套于第 j 个变量，则 M(i,j)=1
http://www.ma-xy.com 20 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.3 方差分析
表 1.4: 多因素一元方差分析模型表
’model’ 参数取值 模型说明
只对 N 个主效应进行检验，不考虑交互效应，这是默认情况。例如考
虑三因素 (A,B,C) 的试验，此模型对主效应 A,B,C 进行检验
’interaction’ 对 N 个主效应和 C
2
N
个两因素进行检验。例如考虑三因素 (A,B,C)
的试验，此模型对主效应 A,B,C 和两因素交互效应 AB,AC,BC 进行
检验
’full’ 对 N 个主效应和全部交互效应进行检验。例如考虑三因素 (A,B,C) 的
试验，此模型对主效应 (A,B,C)、两因素交互效应 (AB,AC,BC) 和三
因素交互效应 ABC 进行检验
正整数 k(k ≤ N ) 当 k=1 时，此模型相当于’linear’ 模型；当 k=2 时，此模型相当
于’interaction’ 模型；当 2<k<N 时，对 N 个主效应和 2 至 k 因
素交互效应进行检验；当 k=N 时，此模型相当于’full’ 模型；
由 0 和 1 构成的矩阵 用矩阵精确控制模型中的效应项。例如考虑三因素 (A,B,C) 的试验，
有以下对应关系：[100] 主效应 A；[010] 主效应 B；[001] 主效应 C；
[110] 交互效应 AB；[101] 交互效应 AC；[011] 交互效应 BC；[111] 交
互效应 ABC。若矩阵为 [010;001;011]，则表示模型中有主效应项 B,C
以及交互效应项 BC
在多重比较中，将参数’dimension’ 设置为向量，用来指定要比较的分组变量序号。
4). 单因素多元方差分析的 matlab 命令为
[d,p,stats] = manova1(x,group,alpha)
其中：d 取值 0 时，接受原假设；取值为 1 时，拒绝原假设但无法拒绝共线性假设；取值为 2 时，
拒绝原假设，拒绝共线假设，但无法拒绝共面假设。p 为 p 值向量，它的第 i 个元素对应的原假
设为：各组均值向量位于一个 i − 1 维空间。stats 结构体包含的字段如表 (1.5) 所示
http://www.ma-xy.com 21 http://www.ma-xy.com
http://www.ma-xy.com
1.4 分布检验 第一章 统计基础
表 1.5: manova1 结构体字段表
字段名 说明
W 组内平方和及交叉乘积和矩阵
B 组间平方和及交叉乘积和矩阵
T 总平方和及交叉乘积和矩阵
dfW W 的自由度
dfB B 的自由度
dfT T 的自由度
lambda Wilk’s Λ 检验统计量的观测值向量，第 i 个元素对应的原假设是：各组的均值向量
位于一个 i-1 维空间
chisq 由 Wilk’s Λ 检验统计量转换得到的近似卡方检验统计量的观测值向量
chisqdf 卡方检验统计量的自由度
eigenval 使 W
−1
B 的特征值
eigenvec 使 W
−1
B 的特征向量，它们是典型变量 C 的系数向量，经过了单位化，使得典型
变量的组内方差为 1
canon 典型变量 C，等于 XC*eigenvec ，其中 XC 是将 X 按列中心化 (每列元素减去每
列的均值) 后得到的矩阵
mdist 每个点到它的组均值的马氏距离 (Mahalanobis 距离) 向量
gmdist 各组组均值的之间马氏距离矩阵
1.4 分布检验
在前面的单总体参数估计和检验，或者方差分析中，我们都假设总体是正态的，那么，如何
检验假设是否为正态呢？或者更广泛的说，如何依据样本估计和检验总体的分布呢？下面，我们
来解决分布检验问题。我们先从正态分布检验入手，再讨论指数分布，最后扩展到任意分布的检
验。
常用的分布检验方法有：1. 卡方拟合优度检验 (1900)；2. Jarque-Bera 检验 (J-B 正态检验)；
3. Kolmogonov-Smirnov 检验；4. lillietest 检验 (1967)；5. Shaoiro-Wilk 检验，要求样本数
8 ⩽ n ⩽ 50(1965，正态)；6.Epps—Pully 检验，要求 n ⩾ 8(正态)。对于指数分布的检验，有卡
方拟合优度检验和格列坚科检验等检验方法。
1.4.1 正态分布检验方法
一提到正态性检验，我们会第一时间想到 JBtest，是的，JBtest 是可行的，但这里我们先
不去了解有哪些检验方法，我们来试着自己考虑一下。我们假设 x 服从正态分布，设置原假
设为 H0 : f (x) = N (µ, σ
2
)。但是正态分布的参数 µ, σ 未知，这使得我们根本无从下手，因
为 H0 中，我们需要一个具体的正态分布。我们用样本均值和方差来代替 µ, σ，把假设变为
http://www.ma-xy.com 22 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.4 分布检验
H0 : f(x) = N(ˆµ, ˆσ
2
)。如果 x 不服从 N(ˆµ, ˆσ
2
)，那么 x 更不可能服从其它均值方差的正态分布。
有了 ˆµ, ˆσ
2
，也就有了 N(ˆµ, ˆσ
2
)，我们将样本 x 经验密度 f
n
(x) 带入 (关于经验密度，可以使用
直方图或者其它)，可能得到像图 (1.9) 那样的情况
(a) (b) (c)
图 1.9: 正态分布检验图
如果假设 H0 成立，则会有图1.9(a) 的情况，如果 H0 不成立，则可能会有图1.9(b)(c) 的情
况。我们知道直方图 F
n
(x) 和分割窗宽 h 有关，在 h 给定后，我们可以构建如下统计量
T =
n
i=1
(f(x
i
) − f
n
(x
i
))
2
·
∼ χ
2
(n − 1)
离差平方和服从 χ
2
分布是合适的。对于超参数 (外来参数)h，我们如何确定 h 呢？我们可以让
h 逐步逼近 0，然后在每个 h 下做检验，记录每次的检验效果，选取最好的。
下面，我们介绍标准的正态分布检验方法：JBtest、Shapiro-Wilk 检验和 Epps-Pully 检验。
这里，我们不对统计量的分布进行证明。
JBtest 是基于这样的想法：如果 H0 : f = N (µ, σ
2
) 成立，正态总体的偏度为 0，峰度为 3，
那么样本的偏度应该接近 0，峰度应该接近 3。为此，构造检验统计量为
T ≜ JB =
n
6
s
2
+
(k − 3)
2
4
·
∼ χ
2
(2)
其中：s 为样本偏度 (skewness)，k 为样本峰度 (kurtosis)。对于 JB
·
∼ χ
2
(2)，其实是这样的
JB =
n
6
(s − 0)
2
+
(k − 3)
2
4
·
∼ χ
2
(2)
这样写的话，JB
·
∼ χ
2
(2) 也就合理了。
JBtest 的 matlab 命令为
[h,p,jbtest,critral] = jbtest(x,alpha,mctol)
其中：h=0 表示在 alpha 显著水平下无法拒绝原假设 H0；p 为 p 值 (统计量的右向累计概率值)；
jbtest 是 JB 统计量的值；critual 是 χ
2
α
2
的值，用于和 jbtest 比较；alpha 是规则 (显著水平)；
mctol 利用 MC(蒙特卡洛模拟) 计算 p 值的近似值，当 alpha 或者 p 值不在 [0.01�0.5] 上，需要
进行 p 值的近似，p 值的标准误为
p(1 − p)
mcreps
< mctol
其中：mcreps 是重复模拟的次数。
http://www.ma-xy.com 23 http://www.ma-xy.com
http://www.ma-xy.com
1.4 分布检验 第一章 统计基础
值得一提的是，JBtest 的性质并不稳定，当样本中有极端值时，检验结果经常发生错误，因
此，在进行 JBtest 之前，最好将样本数据整理一下，例如：去掉 1/4 分位数之外的样本。我们
姑且认为 JBtest 这种在样本异常时表现的不稳定性为鲁棒性，鲁棒性是前面提到过的所有方法
都应该考虑的问题。
Shapiro-Wilk 检验于 1965 年提出，可以参考《正态性检验》梁小筠。SWtest 对样本量有
要求，要求样本量在 8 到 50 之间 (8 ⩽ n ⩽ 50)，其原假设 H0 如前，检验统计量为
T
sw
=
n
i=1
a
i
x
(i)
2
n
i=1
(x
(i)
− ¯x)
2
=
n
i=1
(x
(i)
− ¯x)(a
i
− ¯a)
2
n
i=1
(x
(i)
− ¯x)
2
n
i=1
(a
i
− ¯a)
2
=
[n/2]
i=1
a
i
(x
(n+1−i)
− x
(i)
)
2
n
i=1
(x
(i)
− ¯x)
2
其中：a
′
= m
′
v
−1
/
√
m
′
v
−1
v
−1
m；m, v 的来源是 x
(i)
= µ + σm
i
+ a
i
。
对于上述统计量 T
sw
，我们可能并不能够得到它的具体分布，但是我们仍然可以通过模拟来
得到它的近似分布，或者用其他的方法都可以。接下来，设置决策准则 (阈值)α ，将样本统计量
T
sw
(x) 和阈值统计量 (1 − α 分位数)T
−1
sw
(1 − α) 相比较。关于 T
sw
的分位数表可以参考《数理
统计学》茆诗松附表 10。
Epps-Pully 检验要求样本量大于 8。其基本思路 (统计量构造方法) 是：样本的特征函数与
正态分布特征函数之差的模的平方的加权积分 (和) 的形式
w
g(f(x)) − g(
ˆ
f(x))
2
或者
w
g(f(x)) − g(
ˆ
f(x))
2
dx
上述两式可以简单理解为离差平方和。我们的原假设 H0 : f(x) = N(µ, σ
2
)，令样本方差为
s
2
=
1
n
n
i=1
(x
i
− ¯x)，构造检验统计量为
T
EP
= 1 +
n
√
3
+
2
n
n
k=2
k−1
j=1
exp
−
(x
j
− x
k
)
2s
2
−
√
2
n
j=1
exp
−
(x
j
− x
k
)
4s
2
这里，我们仍不能确定检验统计量 T
EP
，其 α 分位数表及构造方法可以参考《数理统计学》
茆诗松附表 1。
http://www.ma-xy.com 24 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.4 分布检验
下面我们要讨论的问题是：检验统计量 T
EP
和 T
sw
的分布都并非常见的卡方分布、t 分布
或者 F 分布，那么，我们应该如何求解其分布，或者如何构造其 α 分位数表 (分位数是自变量
轴上的一点，比如 1/2 分位数是自变量轴的中间点)？思路：我们可以利用计算机多次采样，产
生多个样本，每个样本都会有一个统计量具体的取值，当多次抽样后，我们就可以近似给出 T
EP
的密度函数的估计。当然，这是一个估计值，与抽样方法样本数量等皆有关，并且也会有相应的
标准误。
上述思路的关键问题是：我们多次采样的总体不是原假设中的 f (x) = N (µ, σ
2
)，而是一个
分布待检验的总体。有时候，我们并不能多次采样，所以，我们要想办法从样本中生成样本。
关于指数分布的检验，我们这里不做介绍，常见的检验方法有两种：1.χ
2
检验；2. 格列坚
科检验。这两种方法不仅对完全样本适用，对截尾样本也适用，详细内容可参考《数理统计学》
P295。
1.4.2 适用所有分布的检验方法
下面，我们来介绍适用于所有分布检验的两种检验方法：1.Komogrov-Smirnov 检验；2.χ
2
拟合优度检验。
Komogrov-Smirnov 检验
总结前面我们所建立的假设检验，对于假设检验，我们的步骤一般为：1、给出原假设和备
择假设；2、构建检验统计量 (检验统计量很随意，一定要灵活)；3、分析检验统计量的分布；4、
设定阈值 (置信水平 α) 或者计算 p 值，进行统计决策 (给出检验结果)。
Step1. 原假设和备择假设。我们假设 H0 : x 的分布函数为 F (x)；备择假设 H1 : x 的分布不是
F (x)。注意：这里的分布函数 F (x) 是不含参数的确定的函数 (比如参数 µ, σ 已知)。
Step2. 构建检验统计量 T 。我们从分布函数开始，F (x) 是总体的分布函数，F
n
(x) 是样本的分
布函数 (经验分布)
F
n
(x) =
n
i=1
I
i
(x)
n
其中：
I
i
(x) =
1 x
i
⩽ x
0 x
i
> x
∼ b(1, F (x)) iid
总体分布和样本分布情况如图 (1.10) 所示
http://www.ma-xy.com 25 http://www.ma-xy.com
http://www.ma-xy.com
1.4 分布检验 第一章 统计基础
图 1.10: 分布检验图
将 F
n
(x) 视为一个动态曲线，随着样本量 n 的变化而变化，这样有利于理解极限定理 (局部
相似和全局相似)。我们知道 F
n
(x) 是一个统计量，是样本的函数，如果 F
n
(x) 和 F (x) 的最大
离差都很小，则我们就没有理由拒绝原假设。构建如下统计量
T ≜ D
n
= sup
−∞<x<∞
|F
n
(x) − F (x)|
下面，我们来研究一下检验统计量 D
n
的分布情况。关于统计量 D
n
，有格里汶科定理：
P
lim
n→∞
D
n
= 0
= 1
D
n
几乎处处以概率 1 趋于 0。Komogrov 于 1933 年给出了 D
n
的精确分布
P
D
n
⩽ λ +
1
m
=
0, λ ⩽ 0
1
2n
+λ
1
2n
−λ
3
2n
+λ
3
2n
−λ
···
2n−1
2n
+λ
2n−1
2n
−λ
n!dy
1
dy
2
. . . dy
n
, 0 < λ ⩽
2n − 1
2n
1, λ >
2n − 1
2n
D
n
的渐近分布：当 n 较大时，上式就变得不易求解，在 n = 30 时，D
n
渐近分布就和精确
分布相似了。D
n
的渐近分布为
P {
√
nD
n
< λ } =
∞
j=−∞
(−1)
j
e
−2j
2
λ
2
j > 0
有了分布之后，我们就可以设置阈值 α，然后进行假设检验了。K-Stest 的改进方向：上面
提到过的 KS 检验需要有很高的要求，要求 H0 中的分布 F (x) 完全已知，不能有任何参数，如
果含有未知参数，F (x|θ), θ ∈ Θ，则量 D
n
就不再是统计量了，因为统计量只要求是样本的函数，
而不能有外来参数 θ。如果我们用样本来估计参数 θ 时，D
n
= sup
x
|F
n
(x) − F (x|
ˆ
θ)| 的分布是
未知的，这使得检验难以进行。幸运的是，lillie 于 1967 年提出用
ˆ
θ 来代替 θ，然后再进行 KS
检验的 lillie 检验，值得注意的是，lillietest 不适用与非位置尺度分布的检验。这里，我们直接给
出其 matlab 命令，对其原理不做介绍。
[h,p,kstat,crival] = lillietest(x,alpha,dister)
其中：参数 dister 用于指定分布，参数值为’norm’ 时表示正态分布，’exp’ 表示指数分布，’ev’ 表
示极值分布。其余参数和 jbtest 相似，不做说明。
http://www.ma-xy.com 26 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.5 非参数统计
卡方拟合优度检验
接下来介绍一个应用最为广泛的分布检验方法-χ
2
拟合优度检验。从检验的名字也可以看出
该检验时如何进行的，χ
2
拟合优度检验和前面我们刚开始提到的检验有极大的相似性。χ
2
检验
原本是为属性 (分类) 数据服务的，由于连续分布可以离散化，所以其对连续分布检验也同样适
用，但要注意的是：窗宽 h 会严重影响检验结果。卡方拟合优度检验的原假设和最终的决策方法
和之前介绍的方法相同，不同的是统计量的构造。对于统计量的构造
Step1. 将 x 的取值划分为几个区间段。设置 h 为窗宽 (区间段长度)，则有区间段数目为
Range(x)
h
(这
里，我们假设 Range(x) 可以被 h 整除)，Range(x) = max(x) − min(x)。对于外来参数 h 的确
定，可以参考前面介绍的方法，也可以参考后面非参数统计部分核密度估计中 h 的确定方法。
Step2. 统计落入各区间段的样本数量 n
i
，这里的 n
i
是实际频数。如果 x
i
落在区间段边界上，
我们采用“计上不计下”的原则来计数。
Step3. 计算落入各区间段的理论概率 p
i
和理论频数 np
i
。设第 i 个区间为 [a
i−1
, a
i
]，则
p
i
= P {a
i−1
⩽ x ⩽ a
i
} = F (a
i
) − F (a
i−1
)
并且，如果 F (x) 中包含参数 θ，则用样本先估计 θ，然后再计算理论概率。
Step4. 计算统计量。
T =
[Range(x)/h]
i=1
(n
i
− np
i
)
2
np
i
·
∼ χ
2
([Range(x)/h] − 1)
如果 p
i
是估计值 ˆp
i
，则
T =
[Range(x)/h]
i=1
(n
i
− nˆp
i
)
2
nˆp
i
·
∼ χ
2
([Range(x)/h] − 1)
χ
2
拟合优度检验的改进思路：卡方拟合优度检验的核心工作是样本密度
ˆ
f(x) 的计算 (总体
密度的估计)。在 Step1 中，我们硬性的将 x 划分了区间段，下面，我们换一种思路，在给定窗宽 h
后，我们不去统计落入 [a
i−1
, a
i
] 的样本个数 n
i
，而是统计落入以 x
i
为中心的 [x
i
−h/2, x
i
+ h/2]
的样本个数，然后再计算理论频数，由实际频数减去理论频数求平方后除理论频数，继续构建 χ
2
统计量，只不过这里的自由度变为 ( n − 1)。
χ
2
拟合优度检验的 matlab 命令为
[h,p,stats] = chi2gof(x,param,. . . )
其中：参数’nbins’ 是组数；’ctrs’ 是组中点；’edges’ 是组边界；’cdf’ 的参数值可以是 {@normcdf, µ, σ
2
}；
’nparm’ 分布中未知参数的个数；’emin’ 最小理论频数；’frequency’x 中个元素出现个频数；’alpha’
显著水平。
1.5 非参数统计
前面介绍的估计和检验都是基于参数的，比如单总体均值估计，我们是假设总体服从正态分
布，然后再估计正态分布的均值。这一部分，我们扔掉正态总体的前提，让数据自己说话。在前
http://www.ma-xy.com 27 http://www.ma-xy.com
http://www.ma-xy.com
1.5 非参数统计 第一章 统计基础
面的参数统计部分，我们讨论了单总体参数估计与检验、两总体均值估计与检验和多总体方差分
析等问题。在非参数统计部分，我们仍然简单的介绍一下这些问题：单总体推断 (估计) 和检验、
两总体位置尺度推断和检验以及多总体方差分析。
前面介绍的许多方法都是基于正态总体的，如果我们对总体不是很了解，不能确定其是否为
正态，或者经过分布检验后，确定总体分布不是正态，那么这些方法就失效了。现在的问题是：
如果对总体分布不确定，或者总体不是正态分布，那么我们就不能估计均值、检验均值了吗？可
以的，我们可以根据均值等统计数字特征的最原始定义和性质来进行估计和检验，比如：如果我
们要检验总体的中位数 M
0.5
是否为 q
0
，可以依据中位数的原始定义和性质，将样本数据排序后，
有一半样本在 q
0
左边，一半在 q
0
右边，那么我们就不能拒绝原假设。同样，我们可以检验总体
的分位数。
1.5.1 单总体推断与检验
符号检验
我们先来介绍用于中位数检验的符号检验方法
°
。我们要检验总体的中位数 M
0.5
是否为 q
0
。
Step1. 首先设定原假设 H0 : M
0.5
= q
0
，备择假设 H1 : M
0.5
= q
0
。
Step2. 然后再来构建检验统计量，当 H0 成立的时候，样本应该有一半落在 q
0
左边，一半落在
q
0
的右边，即 x
i
落在 q
0
左边的概率为
1
2
，既然如此，我们可以记录一下落在 q
0
左右的数据量，
记
s
+
=
I
(x
i
<q
0
)
s
−
=
I
(x
i
>q
0
)
令 s = min{s
+
, s
−
}，有 s ∼ b(n, q
0
)
L
∼ N(nq
0
, nq
0
(1 − q
0
))。我们设定检验统计量为 T = s ∼
b(n, q
0
)。
Step3. 最后我们来进行决策，我们可以根据 p 值进行决策，也可以设定显著水平 α。
Wilcoxon 符号秩检验
在介绍 Wilcoxon 符号秩检验之前，我们先介绍一下样本秩的概念。我们给样本 x
i
秩 R
i
，
比如有 5 个样本：x
1
= 0.38, x
2
= 0.23, x
3
= 0.51, x
4
= 0.40, x
5
= 0.15，那么，他们对应的秩为
R
1
= 3, R
2
= 2, R
3
= 5, R
4
= 4, R
5
= 1。可以看出，R 是一个统计量，R
i
也是一个统计量，它
标志着样本 x
i
在所有样本中的大小，其均值和方差为
E(R
i
) =
n + 1
2
V ar(R
i
) =
n
2
− 1
12
Cov(R
i
, R
j
) = −
n + 1
12
°
单总体中位数检验的 R 命令为 binom.test。
http://www.ma-xy.com 28 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.5 非参数统计
R
i
只是 x
i
的一个得分，当然，这个得分可以定义成其他形式，比如我们在 R
i
外面套一层
函数 a(R
i
)，示例：正态得分
a(r) = Z
−1
1
2
+
r
2(n + 1)
Step1. 原假设 H0 : 单峰分布 F (x) 的对称中心为 θ。
Step2. 构建统计量。前面构建的符号统计量 s 仅考虑了样本点的大小，而未考虑其绝对值大小，
但其绝对值的大小有时候是很重要的，例如对样本 −0.21, −0.2, −0.13, −0.01, 0, 15, 50, 100, 150
来说，0 是中位数，正号负号数目一样，如果只看秩而不看数据，给人的印象是一个很对称的样
本。但实际则不然，问题出来样本的绝对值大小没有考虑。假设样本为 x
1
, . . . , x
n
，样本的顺序
(统计量) 为 x
(1)
, . . . , x
(n)
，样本绝对值 (统计量) 为 |x
1
|, . . . , |x
n
|，样本绝对值的顺序统计量为
|x|
(1)
, . . . , |x|
(n)
，用 R
+
j
表示 |x
j
| 在绝对值样本中的秩，则原假设成立的情况下，有
I{x
i
> θ} · R
+
i
=
I{x
i
< θ}R
+
i
令 u
i
= I{x
i
> θ}，定义 Wilcoxon 符号秩统计量为
W
+
=
u
i
· R
+
i
它是正的样本点按绝对值所得的秩的和。W
+
不仅有样本的正负信息，而且还有样本的数值大小
信息，和符号检验相似，我们在 W
+
较小的时候拒绝原假设，认为总体的对称中心不为 θ。关于
检验统计量 W
+
的分布，可以参考《非参数统计》王星 P75。在大样本下，W
+
具有近似正态性
E(W
+
) =
1
2
n
i=1
i =
1
4
n(n + 1)
V ar(W
+
) =
1
4
n
i=1
i
2
=
1
24
n(n + 1)(2n + 1)
于是有
W
+
− E(W
+
)
√
W
+
L
−→ N (0, 1), n → ∞
上述证明可以参考《高等数理统计》P251。Wilcoxon 符号秩检验的 R 程序为 Wilcox.test。
MATLAB 用 ranksum 进行 Wilcoxon 秩和检验，用 signtest 做成对样本的符号检验，用 signrank
做成对样本的 Wilcoxon 符号秩检验。
1.5.2 两总体位置与尺度推断与检验
对独立两样本位置检验，可以使用 Broon-Mood 检验和 Wilcoxon-Mann-whithey 秩和检验，
其 R 命令为 wilcox.test(x,y,alt = ”less”)。对成对两样本位置检验，可以使用 Wilcoxon 检验，其
R 命令 wilcoxon(x,y,paired = T,alt = ”less”)。对独立两样本尺度检验，可以使用 Mood 检验方
法和 Moses 检验方法，在 R 包 stats 和 fBsics 有命令 mood.test(x,y,alt = ”greater”)。
http://www.ma-xy.com 29 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
1.5.3 非参数方差分析
单因素方差分析的 K-S 检验
Kruskal-Wallis 检验用于单因素方差分析，要求样本量足够大，其原假设为 H0 : k 个总体位
置相同。其检验统计量为
H =
12
n(n + 1)
k
j=1
R
2
:j
n
j
− 3(n + 1) ∼ χ
2
(k + 1)
其中：n 为样本总数，k 为因素水平数，R
:j
=
n
j
i=1
R
ij
为第 j 组样本的秩和，n
j
为第 j 组样本
数。
Kruskal-Wallis 检验的 R 命令为 kruskal.test，其 MATLAB 命令为
[p,table,stats] = kruskalwallkis(x,group,displayopt)
双因素方差分析的 Friedman 检验
Friedman 检验要求各总体分布的形式相同，位置参数可以不同。其原假设为 H0 : µ
i
= µ
j
，
检验统计量为
Q =
12
nk(k + 1)
k
j=1
R
2
:j
− 3n(k + 1) ∼ χ
2
(k − 1)
其中：n 为样本数，k 为组数，R
:j
为第 j 组样本的秩和。Friedman 检验的 R 命令为 fried-
man.test(y,groups,blocks)，其 MATLAB 命令为
[p,table,stats] = friedman(x,reps,displayopt)
思考：对于同一个原假设，有不同的检验方法，那么如何评价各检验方法的功效呢？
1.6 相关性分析
1.6.1 相关性简介
相关性分析是统计分析的一个重点，前面我们分析的一般是单总体 (单个变量) 的估计与检
验问题，在方差分析中分析了分类变量对连续变量的影响 (两变量或多变量的相关性)。后面的内
容我们要讨论变量相关性问题，比如：连续变量和连续变量的相关性、分类变量和连续变量的相
关性、分类变量和分类变量的相关性、有序分类变量的相关性以及多变量之间的相关性。当然，
相关性分析是具有一定假设和前提的，比如常见的 Pearson 相关性分析就要求变量是正态的，并
且这种分析仅能检验两连续变量之间是否具有线性关系，如果变量之间是其它关系，Pearson 就
变得无能为力了。
这里顺便说一下数据分析的一般流程。如果数据已经规整好了 (结构化数据)，我们的分析
步骤一般是：1、设置变量属性和变量值标签，缺失值、异常值的处理。常见的变量属性为：连
续变量、二分类变量、多分类变量以及有序变量，其实就是 SPSS 中的变量属性设置。变量值标
http://www.ma-xy.com 30 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
签即为分类变量的变量值注释。2、变量的描述统计 (均值、方差以及极大极小值等)。3、单总体
(单变量) 估计与检验。4、变量相关性检验。分类变量对连续变量的 ANOVA，连续变量对连续
变量的 Pearson 相关性以及 Kendall 秩相关分析 (非参数)，分类变量对分类变量的列联表分析，
有序变量对有序变量的列联表分析和 spearman 秩、kendall 秩相关分析。5、建模。设置目标变
量，即 y，y 可以是一个变量也可以是多个变量，可以是连续变量也可以是分类变量，挑选自变
量 x
1
, x
2
. . . ，建立变量关系模型 y = f(x)。6、模型的检验与比较。关于定题和结论这里不介绍。
下面，我们将介绍一些变量相关性检验的方法，在这一部分后面的章节中，我们将介绍变量关系
模型：线性回归、SVM、RVM 和神经网络 ANN 等。
1.6.2 连续变量对连续变量的相关性
首先，我们先来介绍用于连续变量相关性检验的 Pearson 积矩线性相关系数 ρ、Kendall 秩
相关系数 τ 和 Spearman 秩相关系数 ρ
s
。设我们现在要考虑变量 x 和 y 的关系，有样本数据
x
1
, x
2
, . . . , x
n
和 y
1
, y
2
, . . . , y
n
。
Pearson 相关系数 在变量 x, y 为正态总体的假设前提下，我们有 x, y 的 Pearson 线性相关系
数
ˆρ =
n
i=1
(x
i
− ¯x)(y
i
− ¯y)
n
i=1
(x
i
− ¯x)
2
n
i=1
(y
i
− ¯y)
2
其中：¯x =
n
i=1
x
i
, ¯y =
n
i=1
y
i
是极大似然估计，由极大似然估计得不变性，所以 ˆρ 也是总体
相关性的极大似然估计。关于极大似然估计的不变性，可以参考《数理统计学》。
Kendall 秩相关系数 是 Maurice Kendall 于 1938 年提出的方法，常用希腊字母 τ 表示。设有
两个随机变量 x, y，x, y 可以是 (有序) 分类变量也可以是连续变量，设样本量为 n。定义 (x
i
, y
i
)
为一个样本对，样本对的一致性定义如下：
定义 (样本对的一致性) 当 x
i
> x
j
并且 y
i
> y
j
，或者 x
i
< x
j
并且 y
i
< y
j
，即 (x
i
−x
j
)(y
i
−
y
j
) > 0，则称 x
i
, y
i
一致；当 (x
i
−x
j
)(y
i
−y
j
) < 0，则称 x
i
, y
i
不一致；当 (x
i
−x
j
)(y
i
−y
j
) = 0，
称 x
i
, y
i
既不是一致的也不是不一致。
记 G 为样本中样本对一致的样本对数量，H 为样本中样本对不一致的样本对数量，则 Kendall
秩相关系数的计算公式为
τ =
G − H
1
2
n(n − 1)
注意：上述公式仅适用于样本 x, y 中不存在相同元素的情况。对于样本中存在相同元素的情
况，比如：x 的样本值为 1, 2, 3, 4, 3, 3, 2，我们就需要对 Kendall 计算公式进行修正，修正计算公
http://www.ma-xy.com 31 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
式为
τ
′
=
G − H
(n
3
− n
1
)(n
3
− n
2
)
其中：
n
3
=
1
2
n(n + 1)
n
1
=
s
i=1
1
2
u
i
(u
i
− 1)
n
2
=
t
i=1
1
2
v
i
(v
i
− 1)
对于上面的 s, t 的计算，以 s 为例，将 x 中的相同元素分别组合成小集合，s 表示集合 x 中拥有
的小集合数，例如 x 包含元素 1, 2, 3, 4, 3, 3, 2，那么这里得到的 s 则为 2，因为只有 2 和 3 有相
同元素，u
i
表示第 i 个小集合所包含的元素个数。
Spearman 秩相关系数 属于非参数方法，从其名称中也可以看出它是一种基于秩的方法。其
定义为
ˆρ
s
=
n
i=1
(R
i
−
¯
R)(Q
i
−
¯
Q)
(R
i
−
¯
R)
2
(Q
i
−
¯
Q)
2
其中：
¯
R =
1
n
n
i=1
R
i
¯
Q =
1
n
n
i=1
Q
i
由于
n
i=1
R
i
=
n
i=1
Q
i
=
n(n + 1)
2
n
i=1
R
2
i
=
n
i=1
Q
2
i
=
n(n + 1)(2n + 1)
6
所以，Spearman 秩相关系数也可以写为
ˆρ
s
= 1 −
6
n(n
2
− 1)
n
i=1
(R
i
− Q
i
)
2
可以发现，Spearman 秩方法和 Pearson 方法很相似。上面介绍的三种相关性度量方法中，后两
种属于非参数统计方法，第一种是参数统计方法。MATLAB 中的命令为
coe = corr(x,y,’type’,’Kendall’)
其中：参数’type’ 指定了相关系数的类型，其参数值可以是’Pearson’、’Spearman’ 和’Kendall’。
http://www.ma-xy.com 32 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
互信息 是信息论中的概念，这个概念我们在后面还会用到。Shannon 于 1948 年首次定义了互
信息 (mutual information)，用于度量两个变量间的相互依赖程度，其定义为
I
(
x, y
) =
p
(
x, y
)
log
p(x, y)
p(x)p(y)
d
x
d
y
其中：p(x, y) 为联合概率密度，p(x), p(y) 为边际密度，log 的底数可以取 2、10、e 或者其它。互
信息的样本估计为
I(x, y) =
i
j
p(x
i
, y
j
) log
p(x
i
, y
j
)
p(x
i
)p(y
j
)
由于总体概率密度未知，所以我们一般才用样本的概率密度来替代总体的密度，相应的样本
密度计算方法可以采用非参数统计中的核密度估计和 K 邻近距离 KNN。
距离相关性 是 Szkely 于 2007 年提出的一种相关性方法，设 X 是 p 维随机向量，Y 是 q 维随
机向量，X, Y 具有有限 1 阶矩。定义 X, Y 总体距离协方差为
V (X, Y ) =
V
2
(X, Y ) =
||f
X,Y
(t, s) − f
x
(t)f
Y
(s)||
2
=
1
C
p
C
q
|f
X,Y
(t, s) − f
x
(t)f
Y
(s)|
2
|t|
1+p
p
|s|
1+q
q
dtds
其中：f
X
(t) 是 X 的特征函数，f
XY
(t, s) 是 X, Y 的联合特征函数，
C
d
=
π
1+
d
2
Γ
1+d
2
, d = p or q
总体距离方差为
V (X, X) =
V
2
(X, X) =
||f
X,X
(t, s) − f
x
(t)f
X
(s)||
2
总体距离相关性为
R(X, Y ) =
R
2
(X, Y ) =
V
2
(X, Y )
V
2
(X)V
2
(Y )
V
2
(X)V
2
(Y ) > 0
0 V
2
(X)V
2
(Y ) < 0
MIC 相关性 是哈佛大学 Broad 研究院的 Reshef 在 2011 年提出的一种数据相关性挖掘方法，
MIC 是 maximal information Coecient 的缩写，其 R 命令在 R 包 minerva 中，命令为
mine(x,y,master,alpha,C,n.cores,varthr)
其中：x,y 指定了数据集，alpha 为网格分割的大小 B(n) = n
α
，C 为起点，n.cores 是并行计算，
varthr 是新的样本最小方差。
考虑 x, y 两个变量，样本容量为 n。MIC 算法主要由以下 2 个因素决定：1、网格划分数，
即在给定的数据集形成的散点图上，在 x 轴和 y 轴上分别进行多少次的划分；2、网格划分的位
http://www.ma-xy.com 33 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
置，即如果在 x 轴上划分 n
x
次，那么这 n
x
个划分点是等距放置还是以某种其他方式放置在 x
轴上？若给定划分数和划分位置，则给定了一种划分，计算该划分下的互信息值
I(x, y) =
p(x, y) log
p(x, y)
p(x)p(y)
dxdy
其中：p(x, y) 用落入格子中的样本频率来估计，p(x), p(y) 用落入 (k, k + 1) 和 (l, l + 1) 区间的
样本频率估计 (0 ⩽ k ⩽ n
x
− 1, 0 ⩽ l ⩽ n
y
− 1)，n
x
, n
y
为 x, y 轴划分段数。
如果固定网格划分数为 n
x
, n
y
，则通过改变网格划分位置，会得到不同的互信息值，记其中
最大的互信息值为 I
n
x
×n
y
(x, y)。进一步，为了方便在不同维度下进行比较，将其标准化，使其
取值在 [0, 1]
M
n
x
×n
y
(x, y) =
I
n
x
×n
y
(x, y)
log(min{x, y})
设定网格划分数目的上限 B(n)，则最大互信息 MIC 定义为
MIC(x, y) = max
n
x
×n
y
<B(n)
{M
n
x
×n
y
(x, y)}
1.6.3 分类变量对分类变量的相关性
上面介绍了几种用于处理连续变量相关性的方法，下面，我们来处理分类变量相关性问题。
四格表
四格表即 2 × 2 列联表，用于处理变量 A, B 皆是二分类变量的情况。四格表一般形式如表
(1.6) 所示
表 1.6: 四格表
B
¯
B
A n
11
n
12
n
1+
¯
A n
21
n
22
n
2+
n
+1
n
+2
n
例如：A 表示是否吸烟，B 表示是否得肺癌。表中的 A 表示吸烟，
¯
A 表示不吸烟，A 行 B 列的
格 n
11
表示吸烟得肺癌的人数，n
1+
表示吸烟的人数，n 为总样本数，n = n
1+
+n
2+
= n
+1
+n
+2
。
我们一般默认把原因变量“吸烟”作为行，把结果变量 (被影响)“肺癌”作为列。
下面，我们主要说明完全随机四格表，即 n
11
, n
12
, n
21
, n
22
都是随机变量。首先，要说明的
是，对四格表而言，不相关和独立是等价的，其证明可参考《属性数据分析》王静龙 P45。假设
n
ij
服从泊松分布 n
ij
∼ p(λ
ij
)
定理 若 n
ij
∼ p(λ
ij
)，并且相互独立，则在 n = n
11
+n
12
+n
21
+n
22
给定后，(n
11
, n
12
, n
21
, n
22
)
的条件分布为多项分布
(n
11
, n
12
, n
21
, n
22
) = M (n, p
11
, p
12
, p
21
, p
22
)
http://www.ma-xy.com 34 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
其中：
p
ij
=
λ
ij
ij
λ
ij
证明 由泊松分布的可加性，有
n ∼ p(λ
11
+ λ
12
+ λ
21
+ λ
22
)
所以，当 n 给定后，(n
11
, n
12
, n
21
, n
22
) 的条件分布为
λ
n
11
11
n
11
e
−λ
11
λ
n
12
12
n
12
e
−λ
12
λ
n
21
21
n
21
e
−λ
21
λ
n
22
22
n
22
e
−λ
22
(
∑
ij
λ
ij
)
n
n!
e
−
∑
ij
λ
ij
=
n!
n
11
!n
12
!n
21
!n
22
!
p
n
11
11
p
n
12
12
p
n
21
21
p
n
22
22
由于四格表的独立性和不相关等价，所以四格表的独立性检验 (A, B 之间不相关) 的原假设
为
H0 : p
ij
= p
i+
p
+j
∀ij ∈ {1, 2}
H1 : p
ij
= p
i+
p
+j
∃ij ∈ {1, 2}
接下来的问题就是构建检验统计量，下面我们介绍两种检验统计量：1、卡方检验统计量 2、
似然比检验统计量。
(1) 构建卡方统计量
T =
2
i=1
2
j=1
(n
ij
− np
ij
)
2
np
ij
其中：n
ij
为实际频数，np
ij
为期望频数。由于当然，总体的 p
ij
未知，所以用样本估计 ˆp
ij
，于
是统计量变为
T =
2
i=1
2
j=1
(n
ij
− nˆp
ij
)
2
nˆp
ij
下面讨论总体概率 p
ij
的估计 ˆp
ij
。在 H0 中，我们假设 p
ij
= p
i+
p
+j
，而总体概率 p
1+
, p
2+
, p
+1
, p
+2
的极大似然估计 (频率估计) 为
ˆp
1+
=
n
1+
n
, ˆp
2+
=
n
2+
n
, ˆp
+1
=
n
+1
n
, ˆp
+2
=
n
+2
n
所以，ˆp
ij
为
ˆp
ij
= ˆp
i+
ˆp
+j
=
n
i+
n
+j
n
2
http://www.ma-xy.com 35 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
于是检验统计量变为
T =
2
i=1
2
j=1
(n
ij
− nˆp
ij
)
2
nˆp
ij
=
2
i=1
2
j=1
n
ij
−
n
i+
n
+j
n
2
n
i+
n
+j
n
=
2
i=1
2
j=1
n
2
ij
n
i+
n
+j
n
− n
=
n(n
11
n
22
− n
12
n
21
)
n
1+
n
2+
n
+1
n
+2
和前面分布估计中的卡方拟合优度检验一样，检验统计量 T 服从卡方分布，并且，由于 p
1+
+p
2+
=
p
+1
+ p
+2
= 1 ，所以未知参数只有 2 个，卡方分布的自由度为 df = 4 − 1 − 2，即
T
·
∼ χ
2
(1)
经过连续性修正的检验统计量为
T
′
=
n
|n
11
n
22
− n
12
n
21
| −
n
2
n
1+
n
2+
n
+1
n
+2
(2) 构建似然比检验统计量
−2 ln(Λ) = −2
2
i=1
2
j=1
n
ij
ln
ˆp
ij
n
ij
/n
由
ˆp
ij
=
n
i+
n
+j
n
2
于是似然比检验统计量写为
−2 ln(Λ) = −2
2
i=1
2
j=1
n
ij
ln
n
i+
n
+j
n
ij
n
·
∼ χ
2
(1)
二维列联表
前面介绍的四格表是用于处理分类变量 A, B 是二分类的情况，仿照其表格结构，我们可以
处理 A 是二分类变量 B 是多分类的情况，同样也可以处理 A, B 皆是多分类变量的情况。设有
多分类变量 A, B，A 有 r 个水平，B 有 c 个水平，其形成的二维列联表如表1.7所示
http://www.ma-xy.com 36 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
表 1.7: 二维列联表
B
1
B
2
··· B
c
A
1
— — — — —
A
2
— n
22
— — —
.
.
. — — — — n
i+
A
r
— — — — —
— — n
+j
— —
表中 A
i
表示变量 A 的第 i 个水平，n
ij
表示是 A
i
和 B
j
的频数，n
i+
是 A
i
的样本数量，n
为总的样本数量，p
ij
是 A
i
和 B
j
的概率 (理论)，ˆp
ij
是 p
ij
的样本估计。
现在，我们要检验变量 A, B 是否相关。我们设置原假设为 H0 : ∀i ∈ r, j ∈ c, p
ij
= p
i+
p
+j
，
即变量 A, B 相互独立，P {A
i
∧ B
j
} = P (A
i
)P (B
j
)。下面我们来构建检验统计量，我们有实际
频数 n
ij
，其期望频数为 np
ij
，可以构建卡方统计量
T =
r
i=1
c
j=1
(n
ij
− np
ij
)
2
np
ij
但是，一般情况下 p
ij
是未知的，所以我们要用样本来估计之。设其样本估计量为 ˆp
ij
，于是期望
频数为 nˆp
ij
，卡方统计量 T 变为
T =
r
i=1
c
j=1
(n
ij
− nˆp
ij
)
2
nˆp
ij
接下来讨论 ˆp
ij
的计算。由原假设 H0 : p
ij
= p
i+
p
+j
，所以，ˆp
ij
= ˆp
i+
ˆp
+j
，而 ˆp
i+
=
n
i+
n
, ˆp
+j
=
n
+j
n
，所以，期望频数为
np
ij
≈ nˆp
ij
= nˆp
i+
ˆp
+j
= n
n
i
+
n
n
+
j
n
=
n
i+
n
+j
n
卡方统计量 T 变为
T =
r
i=1
c
j=1
(n
ij
− nˆp
ij
)
2
nˆp
ij
=
r
i=1
c
j=1
n
2
ij
n
i+
n
+j
/n
− n
∼ χ
2
((r − 1)(c − 1))
http://www.ma-xy.com 37 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
我们还可以构建如下似然比统计量
−2 ln Λ = −2
r
i=1
c
j=1
n
ij
ln
ˆp
ij
n
ij
/n
= −2
r
i=1
c
j=1
n
ij
ln
n
i+
n
+j
n
ij
n
上面介绍的四格表或者二维列联表仅用于检验两个变量 x, y 是否具有相关性，并没有给出
x, y 是正相关还是负相关，也没有给出相关性强度。对于 x, y 的正负相关性及相关性强度的检
验，我们可以参考两连续随机变量相关性的 Pearson 相关性、Spearman 相关性和 Kendall 相关
性。在分类变量中，我们考虑如下相关性 (相合性) 问题：随着抑郁程度的增加，自杀倾向是否增
加？即自杀倾向是否和抑郁程度成正比？设自杀倾向为 y，抑郁程度为 x，我们将自杀倾向和抑
郁程度分为 3 个等级：1, 2, 3，标签为低中高。设样本数为 n，抽取样本，可构建如表 (1.8) 形式
的列联表
表 1.8: 顺序变量列联表
1 2 3
1 195 93 34
2 20 27 27
3 26 37 39
对于 x, y 的相关性，我们可以采用前面介绍的 Kendall 相关性
τ =
2
n(n + 1)
Z
其中：Z = G − H，G 为样本中中一致序列的样本对数量，H 为不一致的数量
Z =
1⩽i<j⩽n
sign((x
i
− x
j
)(y
i
− y
j
)) = G − H
不过，在列联表中，我们需要进行一些修改。Kendall 相关系数 τ 修改为
τ =
Z
[n(n − 1)/2 − T
A
][n(n + 1)/2 − T
B
]
其中：T
A
, T
B
为
T
A
=
r
i=1
n
i+
2
=
r
i=1
n
i+
(n
i+
− 1)
2
T
B
=
c
j=1
n
+j
2
=
c
j=1
n
+j
(n
+j
− 1)
2
http://www.ma-xy.com 38 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
Z = G − H，列联表中 G, H 可以写为
G =
r−1
i=1
c−1
j=1
n
ij
r
k=i+1
c
t=j+1
n
kt
=
i<k
j<t
n
ij
n
kt
H =
r−1
i=1
c−1
j=1
n
ij
r
k=i+1
j−1
t=j+1
n
kt
=
i<k
j>t
n
ij
n
kt
我们还可以用 Gamma 系数 γ 和 Somers 系数 d 来衡量 x, y 相关性强弱。Gamma 系数的
计算公式为
γ =
G − H
G + H
Somers 系数只使用与 2 × c 列联表，其计算公式为
d
y |x
=
G − H
n(n − 1)/2 − T
A
d
x|y
=
G − H
n(n − 1)/2 − T
B
上面介绍的相关性方法都是用来衡量 x, y 之间正负相关性强弱的 (估计)，并没有涉及到假
设检验技术。对于相关性检验问题，可以参考《属性数据分析》王静龙 P92。书中说的是：先计
算 τ ，或者 γ, d
x|y
, d
y |x
以及它们的标准误 se(*)，然后在原假设 (H0 : x, y 相互独立) 成立时，检
验统计量 T = τ/se(τ ) 的渐近服从标准正态分布。上面介绍的相关性度量和检验都是针对两个分
类变量而言，可以尝试将其推广到多个分类变量。还可以考虑有缺失数据的情况，或者小样本的
情况。
MATLAB 中用 crosstab 来实现列联表分析，其命令为
[table,chi2,p,labels] = crosstab(x1,...,xn)
其中：table 即为列联表，chi2 为卡方统计量值，p 为 p 值，labels 为标签和标签值，是一个元
胞数组。R 中的列联表有多种实现方式 (毕竟 R 是专业统计工具嘛)。
1. table(var1,var2,. . . ,varN)：使用 N 个分类变量 (因子) 创建一个 N 维列联表。
2. xtabs(formula,data)：根据一个公式和一个矩阵或数据框创建一个 N 维列联表。
3. prop.table(table,margins)：依 margins 定义的边际列表将表中条目表示为分数形式。
4. margin.table(table,margins)：依 margins 定义的边际列表计算表中条目的和。
5. addmargins(table,margins)：将概述边 margins(默认是求和结果) 放入表中。
6. ftable(table)：创建一个紧凑的“平铺”式列联表
table 函数生成的结果对象类型为“table”，可以直接作为 chisq.test 等函数的参数输入，进
行检验等。在计算列联表的时候，一般都希望可以同时获得行、列和总和范围内的百分比，而默
认安装的 R 没有内置的给出这个功能的函数，不过还是提供了 margin.table 和 prop.table 分别
解决这个问题。代码示例
http://www.ma-xy.com 39 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
1 #ta b le 创 建 列 联 表
2 mytable2 <− tabl e (A, B) #A是 行 变 量 ，B是 列 变 量
3 #xtab 创 建 列 联 表
4 xtab (A ~ B, data=mydata , chi s q = TRUE)
5 #或 者
6 mytable3 <− xtabs ( ~A + B, data=mydata , chis q = TRUE) #mydata为 一 矩 阵 ， 需 要 交 叉 分 类 的
变 量 放 在 ”~” 右 侧 ， 以”+”作 为 分 隔 符
7 #margin . ta b l e ( ) 和 prop . ta b le () 函 数 分 别 生 成 边 际 频 数 和 比 例 。
8 margin . table ( mytable2 ,1) #1 指 的 是 第1个 变 量
9 prop . tabl e ( mytable3 , , 2 ) #2 指 的 是 第 2 个 变 量
10 prop . tabl e ( mytable3 ) #各 单 元 格 所 占 的 比 例
11 #addmargins ( ) 函 数 为 这 些 表 格 添 加 边 际 和
12 addmargins ( mytable2 ) #行 与 列 分 别 相 加
13 addmargins ( prop . tabl e ( mytable2) ) #行 与 列 分 别 相 加
14 #使 用 gmodels 包 中 的 CrossTable () 函 数 是 创 建 出 类 似 于SAS的 表 格 ：
15 l i b ra r y (gmodels )
16 CrossTable ( A r t hr i t i s $Treatment , A r t h r i ti s $Improved)
17 #高 维 列 联 表
18 mytable4 <− xtabs ( ~Treatment + Sex + Improved , data = A r t h ri t i s )
19 f t a b l e (mytable4 ) #将 mytable4 三 维 列 联 表 转 化 为 扁 平 形 式
20 margin . table ( mytable4 ,1) #对 第 1 个 变 量 ( Treatment ) 求 和
21 margin . table ( mytable4 , c ( 1 ,3 ) ) #对 第 1 个 变 量 ( Treatment ) 与 第 3 个 变 量 ( Improved ) 作 边 际 频 数
22 f t a b l e (addmargins ( prop . tab l e ( mytable4 , c (1 ,2) ) , 3) )
23
1.6.4 Copula
Copula 简介
从随机变量的独立性入手。如果两个随机变量 x, y 独立，则其联合概率密度 f (x, y) 等
于各自边缘概率密度 f(x), f(y) 的乘积，即 f (x, y) = f(x)f(y)。现在考虑，如果我们有了
f(x, y), f(x), f(y)，那么，在 f(x, y) 中去掉 f(x), f (y)，剩下的量是否就可以用来表示 x, y 之间
的相关性？这种思想可以追溯到 1959 年 Sklar 提出的 Copula 函数。Sklar 的理论是：联合概率
密度可以分解为边缘密度函数和一个 Copula 函数，并且 Copula 函数唯一。
我们先来统一符号：f 表示密度函数；F 表示分布函数；θ 是参数；C 是 Copula 函数；x, y
是随机变量，x 是随机向量。
Copula 在法语中是链接交换的意思。Copula 模型用来研究多变量相关性和求解联合密度
(分布) 的模型。多变量联合密度的求解一直是统计学中的难题，它不像一元变量有正态分布、均
匀分布、指数分布以及各种分布族。我们后面讨论的 Copula 函数隶属于参数统计，但其多变的
参数估计方法又不局限于参数方法。下面，我们先来简单介绍一下 Copula，然后讨论 Copula 建
模步骤。
注：现代金融分析中，组合投资和资产定价都离不开资产相关性的度量。在后面的建模章节
中，我们简单的讨论了关于股票的组合投资策略，其关键点就是收益率序列的建模及 Copula-VaR
相关性度量。
http://www.ma-xy.com 40 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
定义 (Copula 函数的定义) Copula 是一个多元分布函数，其边缘分布是定义在 [0, 1] 上的
均匀分布，即 U (0, 1)，常用 C 表示，并且满足一下三个条件 1). C : [0, 1]
n
→ [0, 1]；2). C 是
有 grounded 的递增函数；3). C 的所有边缘分布函数 C
i
满足：C
i
= C(1, . . . , 1, u, 1, . . . , 1) =
u, u ∈ [0, 1], i = 1, . . . , n。
设 F
1
, . . . , F
n
是随机变量 x
1
, . . . , x
n
的分布函数，则 C[F
1
(x
1
), . . . , F
n
(x
n
)] 是表示已多变量
的分布函数，其边缘分布函数就为 F
1
, . . . , F
n
。由此可见 Copula 吧一些边缘分布链接为一个联
合分布函数。
定理 (Sklar 定理) 设有 n 个随机变量 x
1
, . . . , x
n
。若 F (x
1
, . . . , x
n
) 是 n 个随机变量的联
合分布函数，其边缘分布为 F
1
, . . . , F
n
，那么一定存在一个 Copula 函数 C，使得
F (x
1
, . . . , x
n
) = C[F
1
(x
1
), . . . , F
n
(x
n
)]
并且，如果边缘分布是连续的，那么 Copula 函数形式唯一。
利用上面的定理，我们可以将一个多维分布函数拆分成边缘分布函数组合的形式
f(x
1
, . . . , x
n
) =
∂F (x
1
, . . . , F
n
)
∂x
1
. . . , ∂x
n
=
∂C[F
1
(x
1
), . . . , F
n
(x
n
)]
∂x
1
. . . , ∂x
n
=
∂C(u
1
, . . . , u
n
)
∂u
1
. . . , ∂u
n
×
i
∂F
i
(x
i
)
∂x
i
= c(u
1
, . . . , u
n
) ×
i
f
i
(x
i
)
= c(˜u) ×
i
f
i
(x
i
)
其中：f (x
1
, . . . , x
n
) 为联合概率密度函数；u
i
= F
i
(x
i
), i = 1 . . . , n；˜u = (u
1
, . . . , u
n
)；c(˜u) 为
Copula
的密度函数。
也就是说，我们可以将联合概率密度 f(·) 拆分为两部分，前部分 c(˜u) 为 Copula 的密度
函数，表示变量 x
1
, . . . , x
n
的相关性；后部分为边缘概率密度的乘积。当各变量相互独立时，
C(u
1
, . . . , u
n
) = u
1
× ··· × u
n
。
上面给出了 Copula 函数的定义，下面给出它的性质。Copula 函数有许多优良的性质，比如：
1、在利用 Copula 构建多维随机变量的联合分布 F 时，不需要对边缘分布 F
i
(x
i
) 作任何假设。
2、Copula 函数导出的相关性度量在严格单调增的变换下都不改变，有如下不变性定理
定理 (Copula 不变性) (Nelsen.2006.) 就两变量而言，对随机变量 x, y 做严格单调增变换
x
∗
= h
1
(x), y
∗
= h
2
(y)，相应的 Copula 函数不变，即
C[F
1
(x), F
2
(y)] = C[F
1
(x
∗
), F
2
(y
∗
)]
简写为 C
xy
= C
x
∗
y
∗
。
http://www.ma-xy.com 41 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
前面介绍的 Spearman 秩相关系数 ρ
s
和 Kendall 秩相关系数 τ 都可以用 Copula 函数表示
ρ
s
(x, y) = 12
[0,1]
2
u
1
u
2
dC(u
1
, u
2
) − 3
τ(x, y) = 4
[0,1]
2
C(u
1
, u
2
)dC(u
1
, u
2
) − 1
尾部相关性的 Copula 表示。在金融风险分析中，随机变量的尾部相关性是相当重要的。设
X, Y 为两个随机变量，条件概率 P {X
1
> x
1
|X
2
> x
2
} 反应了股票市场中一种股票涨高时，引
起另一种股票涨高的概率。
定义 (尾部相关性) 设 X, Y 的分布函数为 F, G，X, Y 的上端尾部相关性和下端尾部相关性
定义为
λ
u
= lim
α→1−
P {Y > G
−1
(α)|X > F
−1
(α)}
λ
l
= lim
α→0+
P {Y ⩽ G
−1
(α)|X ⩽ F
−1
(α)}
其中：F
−1
(α) = inf(x|X ⩾ α)，G
−1
(α) = inf(y|Y ⩾ α)。
如果 λ
u
(或 λ
l
) 存在且 λ
u
∈ [0, 1]，那么 X, Y 具有渐近上 (下) 端尾部相关性，如果 λ
u
= 0 ，
则 X, Y 在上 (下) 端尾部独立。尾部相关系数的 Copula 表示形式为
λ
u
= lim
α→1−
¯
C(α, α)
1 − α
其中：
¯
C(u, u) = 1 − 2u + C(u, u)。推导如下：
λ
u
= lim
α→1−
P {Y > G
−1
(α)|X > F
−1
(α)}
= lim
α→1−
P {Y > G
−1
(α), X > F
−1
(α)}
P {X > F
−1
(α)}
= lim
α→1−
1 − P {Y ⩽ G
−1
(α)} − P {X ⩽ F
−1
(α)} + P {Y ⩽ G
−1
(α), X ⩽ F
−1
(α)}
1 − P {X ⩽ F
−1
(α)}
= lim
α→1−
1 − 2α + C(α, α)
1 − α
不同类型的阿基米德 Copula 族对尾部相关性的刻画是不一样的，Gumbel Copula 和 Clayton
Copula 具有非对称性，可以用来描述变量间非对称的相关模式，不同的是 Gumbel Copula 强调
随机变量间具有更高的上端尾部相关性，而 Clayton Copula 则强调随机变量间具有更高的下端
尾部相关性，Frank Copula 则更高强调分布尾部相关性的变化，在其分布的上尾和下尾，变量间
的相关性是对称增长的。
一些 Copula 函数
(1) 二元正态 Copula 函数
C(u, v|ρ) =
Φ
−1
(u)
−∞
Φ
−1
(v )
−∞
1
2π
1 − ρ
2
exp
−
s
2
− 2ρst + t
2
2(1 − ρ
2
)
dsdt
http://www.ma-xy.com 42 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
其中：ρ 为线性相关系数，Φ
−1
为标准正态分布逆函数，对应 x, y 的一个值。
(2) 二元 t-Copula 函数
C(u, v|ρ, k) =
t
−1
k
(u)
−∞
t
−1
k
(v )
−∞
1
2π
1 − ρ
2
1 +
s
2
− 2ρst + t
2
k(1 − ρ
2
)
−(k+2)/2
dsdt
其中：ρ 为线性相关系数，t
−1
k
是自由度为 k 的 t 分布的逆函数。
(3) 二元 Gumbl Copula 函数
C(u, v|θ) = exp
−
(−log u)
1
θ
+ (−log v)
1
θ
θ
θ ∈ [0, ∞]
(4) 二元 Clayton Copula 函数
C(u, v|θ) = max
[u
−θ
+ v
−θ
− 1]
1
θ
,
0
θ
∈
[
−
1
,
∞
]/
{
0
}
(5) 二元 Frank Copula 函数
C(u, v|θ) =
1
θ
log
1 +
(e
−θu
− 1)(e
−θv
− 1)
e
−θ
− 1
θ ∈ (−∞, ∞)/{0}
上述的 (1)(2) 属于椭圆形 Copula 函数组，(3)(4)(5) 属于二元阿基米德 Copula 函数族。
定义 (阿基米德 Copula 函数族) 阿基米德 Copula 函数族是通过算子 φ(一个完全单调函
数) 构造而成的，其表示形式为：
C(u
1
, . . . , u
n
) = φ
−1
(φ(u
1
) + ··· + φ(u
n
))
其中：φ
−1
是 φ 的逆函数。
更多的阿基米德 Copula 参考谢中华4P190 或者 2007. 王红莲1附录 A。阿基米德 Copula 函
数有许多优良的性质，设 C 是一个具有算子 φ 的阿基米德 Copula 函数：
1). C 是对称的。C(u, v) = C(v , u), u, v ∈ [0, 1]；
2). C 满足结合律。C(C(u, v), w) = C(u, C(v, w)), u, v, w ∈ [0, 1]；
3). 对任意正数 k,kφ 也是 C 的算子；
4). C(u, 1) = u, C(1, v) = v，u, v ∈ [0, 1]
生成 Copula 随机数
下面介绍如何生成服从 Copula 函数的随机数。
(1) 正态 Copula 随机数的模拟
Step1. 如果多元正态分布的相关系数矩阵 Σ 是正定的，对其进行 Cholesky 分解，有
Σ = AA
T
Step2. 生成 n 维独立标准正态随机变量 z = (z
1
, . . . , z
n
)
T
。
Step3. 生成变量组 x，x = Az。
http://www.ma-xy.com 43 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
Step4. 确定分量 u
i
，u
i
= Φ(x
i
), i = 1, 2, . . . , n。
Step5. 如果要生成 X, Y 的随机数，只要 x = F
−1
(u), y = G
−1
(v)
(2) 阿基米德 Copula 函数的随机数
Step1. 生成两个独立随机变量 s, q ∼ U (0, 1)。
Step2. 设 t = K
−1
C
(q)。其中：
K
C
(x) = P {C(U, V ) < x} = x −
φ(x)
φ
′
(x)
Step3. 令 u = φ
−1
(sφ(t))，v = φ
−1
((1 − s)φ(t))s。
(3) 多变量模拟的递归算法。除了以上介绍的针对不同类型 Copula 的数据模拟方法外，还
有一个通用的方法：单变量条件分布的递归模拟。首先，定义
C
i
(u
1
, . . . , u
i
) = C(u
1
, . . . , u
i
, 1, . . . , 1), i = 2, 3, . . . n − 1
因此，C
1
(u
1
) = u
1
，C
n
(u
1
, . . . , u
n
) = C(u
1
, . . . , u
n
)。假设 (U
1
, . . . , U
N
)
t
∼ C，那么给定 (U
1
, . . . , U
N
)
t
的前 i − 1 个分量时，U
i
的条件分布函数可以用 i 维边缘密度函数以及导数表示：
C
i
(u
i
|u
1
, . . . , u
i−1
) = P {U
i
⩽ u
i
|U
1
⩽ u
1
, . . . , U
i−1
⩽ u
i−1
}
=
∂
i−1
C
i
(u
1
, . . . , u
i
)
∂u
1
. . . ∂u
i−1
∂
i−1
C
i−1
(u
1
, . . . , u
i
)
∂u
1
. . . ∂u
i−1
假设上式中的分子和分母都存在，由此，我们可以得到基于条件分布的递归算法：
Step1. 生成数据 u
1
∼ U(0, 1)；
Step2. 通过 C
2
(u
2
|u
1
) 计算生成 u
2
；
Step3. 通过 C
3
(u
3
|u
1
, u
2
) 生成 u
3
；
Step4. 重复上面步骤，通过 C
n
(u
n
|u
1
, . . . , u
n−1
) 生成 u
n
，从而得到最终数据 u
1
, . . . , u
n
)
T
∼ C。
上述算法是从 u
1
∼ U(0, 1) 开始的，之后每一步都要计算条件分布的逆函数 C
i
(u
i
|u
1
, . . . , u
i−1
)，
从而得到结果。
Copula 模型的参数估计
在实际应用中，我们关心的是 Copula 函数的选择以及其参数估计。下面，我们介绍一些
Copula 函数的参数估计方法：极大似然估计、2 步估计法和半参数估计法等。
极大似然估计 设有随机变量 X, Y ，样本为 x
i
, y
i
(i = 1, . . . , n)，n 为样本数。X 的边缘分布为
F (X|θ
1
)，边缘密度为 f (X|θ
1
)，Y 的边缘分布为 G(Y |θ
2
)，边缘密度为 g(Y |θ
2
)，θ
1
, θ
2
为参数。
设 X, Y 的联合分布为 H(X, Y )，联合密度为 h(X, Y )。连接函数 Copula 选择 C(u, v|α ) 为分布
函数，c(u, v|α) 为密度函数，有
c(u, v|α) =
∂
2
C
∂u∂v
http://www.ma-xy.com 44 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
α 为函数中的待求参数。由 Sklar 定理，有
H(X, Y |θ
1
, θ
2
, α) = C[F (X|θ
1
), G(Y |θ
2
)|α]
h(x, y|θ
1
, θ
2
, α) =
∂H
∂x∂y
= c[F (X|θ
1
), G(Y |θ
2
)|α] × f(x|θ
1
)g(y|θ
2
)
我们要求参数 θ
1
, θ
2
, α，极大似然估计的思想是样本出现的概率最大。我们有样本 x
i
, y
i
出
现的概率为
h(x
i
, y
i
|θ
1
, θ
2
α) = c[F (x
i
|θ
1
), G(y
i
|θ
2
)|α] × f(x
i
|θ
1
)g(y
i
|θ
2
)
于是样本 {x
i
, y
i
}
n
i=1
出现的联合概率为
L(θ
1
, θ
2
, α|x
i
, y
i
) =
n
i=1
h
i
(
x
i
, y
i
|
θ
1
, θ
2
, α
)
对上式取对数，有
ln L(θ
1
, θ
2
, α|x
i
, y
i
) =
n
i=1
ln c[F, G|α] +
n
i=1
ln f(x
i
|θ
1
) +
n
i=1
ln g(y
i
|θ
2
) (1.2)
最终，我们的目标是在 θ
1
, θ
2
, α 的可行空间 Θ 中找到最优的 θ
∗
1
, θ
∗
2
, α
∗
来使得 ln L 最大，即
arg max
θ
1
,θ
2
,α∈Θ
ln L(θ
1
, θ
2
, α|x
i
, y
i
)
由极大似然估计的性质可知，(θ
∗
1
, θ
∗
2
, α
∗
)
ML
是 (θ
1
, θ
2
, α) 的相合估计，且具有渐近正态性。
两步法 在上面的 ML 估计中，我们直接在参数 θ
1
, θ
2
, α 的可行空间 Θ 中找到最优的 θ
∗
1
, θ
∗
2
, α
∗
，
回顾目标函数 (1.2)，可以看到目标函数有 3 部分组成，如果我们先让后面两部分先达到最大值，
可以求解边缘分布参数
ˆ
θ
1
,
ˆ
θ
2
，然后再将其带入第一部分，可以继续求解 Copula 参数 α，进而使
目标近似最大
ln L =
ln c(θ
1
, θ
2
, α) +
ln f(θ
1
) +
ln g(θ
2
)
我们记这样得到的参数估计为 (θ
∗
1
, θ
∗
2
, α
∗
)
SML
。为什么说上面的目标函数是近似最大呢？是
因为
ˆ
θ
1
,
ˆ
θ
2
使目标的后两部分最优，并不一定使第一部分最优，因而可能是次优解。其实用一个
二维函数的寻优来解释即可，我们先固定 y，在 x 轴上找最优，然后再在 y 轴上找最优，最后的
到的并不一定是全局最优解。在实际中，SML 和 ML 得到的估计结果基本相似，并且 SML 也具
有相合性和渐近正态性，只不过渐近正态性有所减小。
半参数估计 (Canoical Maxinnum likelihood Method)。在上面的求解中，我们在利用 Copula
来计算联合分布函数时，直接运行
H
=
C
[
F, G
|
α
]
http://www.ma-xy.com 45 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
但在大多数情况下，F, G 是未知的，或者有未知参数 θ
1
, θ
2
。对于此问题，随机变量 X, Y 的边
缘分布 F, G 可以通过非参数统计中的核密度、最邻近估计以及小波估计等进行计算，然后以不
存在参数的形式和 Copula 结合，这样，未知参数就只有 α。用 ML 估计 α，有
α
∗
= arg max
α
n
i=1
ln c(u
i
, v
i
|α)
其中：u
i
=
ˆ
F (x
i
), v
i
=
ˆ
G(y
i
)。
2014. 张连增5表明：当我们选择错误的边缘分布时，CML 估计方法较优，TSML 的偏差和
均方误差都很大；在 AIC 准则下，当边缘分布指定错误时，TSML 方法会错误的选择 Copula 函
数。
经验估计法 对于阿基米德 Copula 函数族的参数 α 有如下经验估计法
f(α) = τ = 1 + 4
1
0
φ
α
(t)
φ
′
α
(t)
dt
其中：φ
α
(t) 是阿基米德 Copula 函数的生成元。
Copula 模型的选择
在给定两个随机变量 X, Y 的样本 x
i
, y
i
, i = 1, . . . , n 后，我们可以通过前面介绍的 Copula
参数估计方法得到 X, Y 的联合概率密度，即 C[F,G]。在计算时，可供我们选择的 Copula 函数
C 和边缘分布函数 F, G 多种多样，那么，我们应该如何选择 Copula 函数呢？下面介绍一些准
则，可以用来评价不同的 Copula 函数。
AIC(BIC) 准则 在许多模型选择问题中都会看到 AIC 的身影，比如时间序列建模 ARMA 等，
但它的计算形式使它局限于极大似然估计模型。其计算形式为
AIC = −2 log(M IE) + 2k
BIC = −2 log(MIE) + k log n
其中：M IE 是样本的极大似然值 (在参数 θ
1
, θ
2
, α 估计出来之后，极大似然函数的函数值也就
有了)；n 为样本数；k 为模型中参数的个数。AIC 和 BIC 的值越小，说明模型越好，但是该方
法不能证明 Copula 明显优于另一个 Copula。
L
2
距离最小法 L
2
距离最小法如其名，类似于离差平方和最小。考虑如果我们有了分布函数
H = C[F, G] 函数以及它的估计
ˆ
C，我们可以证明 H 和
ˆ
C 之间的离差平方和最小。但问题是：
如何得到分布 H = C[F, G] 的估计呢？可以利用前面提到过的经验分布来进行，定义 C 的经验
分布为：
http://www.ma-xy.com 46 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
定义 (经验分布) 设 x
i
, y
i
为样本，
ˆ
F (X) 为变量 X 的经验分布，
ˆ
G(Y ) 为变量 Y 的经验分
布，则有联合分布 H = C 的经验分布
ˆ
C
ˆ
C
n
(u, v) =
1
n
n
i=1
I
{F
n
(x
i
)⩽u}
I
{G
n
(y
i
)⩽v }
u, v ∈ [0, 1]
其中：I 为特征函数。
于是我们可以定义 L
2
范数上的距离
d =
n
i=1
ˆ
C
n
(u, v) − C(u, v)
最后，选择使 d 最小的 Copula 函数即可。
假设检验法 这里的参数 α 是分布中的参数，我们可以使用前面介绍的 χ
2
拟合优度检验和 K-S
检验来检验分布 H = C[F , G] 是否合理。在许多统计软件中，我们都会看到，在模型的参数估计
结果后面伴随着参数的检验，有人会问，模型中的参数都估计出来了，还检验什么？其实，对于
一个数据集，即便我们设置了一个非常不合理的模型，也能估计出其参数值，那么，这些参数有
意义吗？或者说这些参数有必要存在吗？这就是为什么结果中会自带假设检验。
下面，我们也用假设检验的思想来检验我们选择的 Copula 模型的好坏。回忆前面的分布拟
合部分，χ
2
拟合优度检验是样本经验分布和假设分布的离差平方和的检验，于是，对于 C 和
ˆ
C
n
，
我们有检验统计量
T =
n
i=1
n
j=1
n
ˆ
C
n
(x
i
, y
j
) − nC[F (x
i
), G(y
j
)]
2
nC[F (x
i
), G(y
j
)]
∼ χ
2
((n − 1)
2
)
其中；n
ˆ
C
n
(x
i
, y
j
) 表示实际频数，nC[F (x
i
), G(y
j
)] 表示期望频数。χ
2
拟合优度检验需要一定的
样本量。用 K-S 检验，有
T = D
n
= max{|
ˆ
C
n
− C|}
自助法 Bootstrap Bootstrap 方法最早由 Efron 提出
±
。Bootstrap 算法通常分为 3 类：¬非参
数 Bootstrap，也是最简单最常用的 Bootstrap 方法，就是将样本进行有放回的随机抽样，以获
取 Bootstrap 样本；如果对总体已经有所推断，可以从推断的总体重复抽样；®残差 Bootstrap，
对回归模型，首先通过估计得到残差，然后对残差进行再抽样。
用 Bootstrap 方法检验 Copula 函数选择是否正确的基本做法是：利用给定的 Copula 生成
X, Y 的“样本”，然后从中进行重复抽样，对每一次抽样的样本计算对数最大似然值 ln L，最后，
在利用对数似然函数值和再抽样模拟值之间的距离构造拟合优度统计量。
Step1. 对数据集采用多个 Copula 函数进行拟合 (估计)，得到模型 C
1
, C
2
, . . . , C
m
，存储每个
Copula 函数的参数以及对数极大似然函数值。
±
下面这个方法有些问题要捋一捋
http://www.ma-xy.com 47 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
Step2. 从每一个 Copula 模型的分布中抽取 (生成) 样本，重新拟合 Copula 模型，存储其对数
似然函数值。
Step3. 重复 Step2 多次。
Step4. 对比基于原始样本数据的对数似然函数值和基于模拟数据的对数似然函数值。记第 k 个
Copula 模型的原始样本数据的对数似然函数值为 λ
k
(记原始样本数据的对数似然函数值向量为
λ = (λ
1
, λ
2
, . . . , λ
m
))，记第 k 个 Copula 模型模拟得到的对数似然函数值向量为 λ
k
。计算第 k
个 Copula 函数样本对数似然和模拟对数似然之间的马氏距离
D
2
k
= ( λ
k
− λ
k
)S
−1
(λ
k
− λ
k
)
T
k = 1 , 2, . . . , m
其中：S 为模拟对数似然的协方差矩阵。由中心极限定理，对数似然函数的联合分布为渐近正态
分布，我们假设 H0 : 第 i 个 Copula 分布为正态分布，在 H0 成立时，有如下统计量
T = mD
2
k
∼ F (m, B − 1)
其中：B 为模拟抽样次数。
Step5. 对每次模拟计算 p 值。
关于 Copula，我们这里介绍的仅是九牛一毛。有两个重要的方向这里不得不提一下：1. 是
极值 Copula 在变量极端值关系的应用，即 x
1
突变与 x
2
突变的关系。2. 是 Copula 在函数优化
方面的应用，因为 Copula 能够刻画 x
1
和 x
2
之间的关系，这使得它能够在最优化方面进行应用，
这一方面的内容可以参考 2011. 王丽芳2。更详实的 Copula 理论可以参考 2010. 吴娟3。
MATLAB 中的 Copula 函数
MATLAB 提供了 copulat、copulastat、copulaparam、copulapdf、copulacdf 和 copularnd 六
个与 Copula 有关的函数，它们可用于多元正态 Copula 函数、多元 t-Copula 函数、二元 Gumbel、
二元 Clayton 和二元 Frank Copula 函数。
(1)copulat 函数根据样本观测数据估计 Copula 函数中的未知参数。其命令格式为
rhohat = copulat(’Gaussian’,u,param)
[rhohat,nuhat,nuci] = copulat(’t’,u,param)
[paramhat,paramci] = copulat(family,u,param)
其中：’Gaussian’、’t’ 和 family 指定了 Copula 函数的类型，family 的阿基米德函数类型可
选’Clayton’、’Frank’ 和’Gumbel’；u 是边缘分布函数值构成的 n × p 的矩阵，p 是变量个数
u
1
, u
2
, . . . , u
p
，n 为样本数，取值范围在 [0, 1]，如果是阿基米德 Copula，则只能是二元 Copula，
即 u 是 n ×2 矩阵。rhohat 是线性相关系数 ρ 的估计，是 p ×p 矩阵。muhat 是 t 分布的自由度
k。nuci 是 k 的置信下限和置信上限。paramhat 是二元阿基米德 Copula 的参数 α。paramci 是
置信上限和置信下限。param 是辅助参数，Example: ’Alpha’,0.01,’Method’,’ApproximateML’，
Method for tting t copula, specied as the comma-separated pair consisting of ’Method’ and
either ’ML’ or ’ApproximateML’.If you specify ’ApproximateML’, then copulat ts a t copula
for large samples by maximizing an objective function that approximates the prole log likelihood
for the degrees of freedom parameter. This method can be signicantly faster than maximum
http://www.ma-xy.com 48 http://www.ma-xy.com
http://www.ma-xy.com
第一章 统计基础 1.6 相关性分析
likelihood (’ML’), but the estimates and condence limits may not be accurate for small to
moderate sample sizes.
(2)copulastat 函数根据 Copula 函数计算相关量 ρ, τ, ρ
s
。
r = copulastat(’Gaussian’,rho,param)
r = copulastat(’t’,rho,nu,param)
r = copulastat(family,alpha,param)
其中：r is the Kendall’s rank correlation. 可以用 param 来指定返回的相关量，可设置为’type’,’Spearman’。
(3)copulaparam 函数利用 Kendall 秩相关或者 Spearman 秩相关计算 Copula 中的参数 ρ 或
α。
rho = copulaparam(’Gaussian’,r)
rho = copulaparam(’t’,r,nu)
alpha = copulaparam(family,r)
其中：r 是 kendall 相关系数或者 Spearman 相关系数，可以用过’type’,’Spearman’ 来指定。nu
是 t 分布的自由度 k。rho 和 alpha 是 Copula 函数的参数。
(4)copulapdf 函数用来计算 Copula 的密度函数值。
y = copulapdf(’Gaussian’,u,rho)
y = copulapdf(’t’,u,rho,nu)
y = copulapdf(family,u,alpha)
其中：y 是密度函数值。示例：
1 u = linspace (0 ,1 , 1 0) ;
2 [ u1 , u2 ] = meshgrid (u , u) ;
3 y = copulapdf ( ’ Clayton ’ , [ u1 ( : ) , u2 ( : ) ] , 1 ) ;
4 su r f (u1 , u2 , reshape ( y , 1 0 ,1 0 ) )
5 xlabe l ( ’ u1 ’ ) , yl a b e l ( ’ u2 ’ )
6
(5)copulacdf 函数用于计算 Copula 的分布函数值。
y = copulacdf(’Gaussian’,u,rho)
y = copulacdf(’t’,u,rho,nu)
y = copulacdf(family,u,alpha)
(6)copularnd 函数用于生成 Copula 随机数。
u = copularnd(’Gaussian’,rho,n)
u = copularnd(’t’,rho,nu,n)
u = copularnd(family,alpha,n)
其中：n 为随机数的个数。u 为 u
1
, u
2
, . . . , u
p
。示例：
1 rng d e f a u l t % For r e pr o d u c i b i l i t y
2 tau = −0.5;
3 rho = copulaparam( ’ Gaussian ’ , tau )
4 u = copularnd ( ’ g aussian ’ , rho ,1 0 0) ;
5 f i g u r e
6 s c a t t e r h i s t (u ( : , 1 ) ,u ( : , 2 ) )
7
http://www.ma-xy.com 49 http://www.ma-xy.com
http://www.ma-xy.com
1.6 相关性分析 第一章 统计基础
http://www.ma-xy.com 50 http://www.ma-xy.com
http://www.ma-xy.com
参考文献
[1] 王红莲. 链接函数 (Coupla) 理论及其在金融中的应用. 上海财经大学. 硕士学位论文.
[2] 2011. 王丽芳. 基于 copula 理论的分布估计算法研究. 兰州理工大学. 博士学位论文.
[3] 2010. 吴娟. 基于 copula 理论的分布估计算法研究. 华中科技大学. 博士学位论文.
[4] 谢中华.MATLAB 统计分析与应用 40 个案例分析.
[5] 2014. 张连增.Copula 的参数与半参数估计方法的比较.
51