所有由左辰发布的文章

关于左辰

左辰, University of Wisconsin Madison, PhD in Statistics, 学术兴趣:数理统计,概率论

强大数定律与康托三分集

首先从博雷尔正轨数定律(Borel’s Normal Number Theorem)说起。众所周知,(0,1]区间上的每一个实数$\omega$都与一列唯一的无穷的二进制展开序列$\{X_k(\omega)\}$一一对应,其中$X_k (\omega)$表示二进制展开的第k位,对应关系为:

$
\omega=\sum_{k=1}^n\frac{X_k(\omega)}{2^k}
$

如果用(0,1]上的勒贝格测度$\lambda$作为概率测度来计算$X_k(\omega)$的分布,我们容易发现

$
\{\omega:X_1(\omega)=1\}=(\frac{1}{2},1]\\
\{\omega:X_2(\omega)=1\}=(\frac{1}{4},\frac{1}{2}]\cup(\frac{3}{4},1]\\
\vdots\\
\{\omega:X_k(\omega)=1\}=(\frac{1}{2^k},\frac{2}{2^k}]\cup(\frac{3}{2^k},\frac{4}{2^k}]\cup\cdots\cup(\frac{2^k-1}{2^k},1]
$

而其中等式右边每一个集合的勒贝格测度都为1/2。由此,$X_k(\omega)$是事件集$\Omega=(0,1]$导出的概率为1/2的0-1分布的随机变量。

事实上,$X_k(\omega)$不仅是同分布的,而且是相互独立的。强大数定律的结论告诉我们:

$
\lambda\{\omega:\frac{1}{n}\sum_{k=1}^n X_k(\omega)\rightarrow \frac{1}{2}\}=1
$

即:如果考虑区间(0,1]上这样一个子集$A$,这个子集的构成元素是所有与均值收敛于 1/2的二进制展开序列相对应的实数,那么这个子集的“长度”,即勒贝格测度为1($\lambda(A)=1$)。这就是博雷尔正轨数定律。

以上所介绍的(0,1]上实数与其二进制展开序列的对应关系,我认为是解释概率空间基本概念的最经典、最直观的模型:它不仅近乎完美地解释了事件集与随机变量的对应关系,进而直观揭示分布与测度的关系,而且证明了如何在同一个概率空间中生成无穷维相互独立的随机变量。

如果对以上模型稍作推广,可以得到一些更有趣的结论。对于任意的$0-\frac{1}{2}$之间的常数$p$,我们考虑$\omega$的“非对称”的无穷二进制展开序列$X_k(\omega)$,其中每一个$X_k(.)$都是如下集合$A_k$的指示函数:

$
A_1=(1-p,1]\\
A_2=(p-p^2,p]\cup(1-2p+2p^2,1-p]\cup(1-p^2,1]\\
\vdots\\
$

由对称性可知,$p>1/2$的情形可以类似得到。表达式虽然复杂,但是构造的原理很直观:首先将区间分为长度为p:(1-2p):p的三段,然后将每一小段安比例p:(1-2p):p继续分割,如此不断进行下去,将第k步得到的3k个区间中每隔两个选取一个合并,如此就可以得到$A_k$。容易得到,这种无穷序列和(0,1]区间的实数也是一一对应的,$X_k$服从相互独立的0-1分布,$P\{X_k=1\}=p$。

这里采用了的分割方法和二进制分割略有区别;当$p=\frac{1}{2}$时,每次分割的中间那段区间是退化的,因此得到的结果与二进制的分割本质相同。在这种形式下,与集合$A$对应的集合$B$包含所有那些与收敛于$p$的“非对称”二进制展开序列对应的实数集合,即:

$
B=\{\omega:\frac{1}{n}\sum_{k=1}^nX_k(\omega)\rightarrow p\}
$

强大数定律告诉我们:

$
\lambda(B)=1
$

当p不为1/2时,我们显然有

$
\lambda\{\frac{1}{n}\sum_{k=1}^nX_k(\omega)\rightarrow \frac{1}{2}\}=0
$

因为大括号内的集合属于集合$B$的补集,所以其测度为0。我们把括号内的集合记为$A’$。注意:虽然$A’$和$A$的表达式相同,但是两个集合对应的二进制展开形式是不同的,因此它们是两个截然不同的集合,所以它们虽然在相同的测度空间内,但有着截然不同的测度。

最后,我们按照如下次序构造集合$C$:

$
C_1=(0,p]\cup(1-p,1]\\
C_2=(0,p^2]\cup(p-p^2,p]\cup(1-p,1-p+p^2]\cup(1-p^2,1]\\
C_3=(0,p^3]\cup(p^2-p^3,p^2]\cup\cdots\cup(1-p^3,1]\\
\vdots\\
C_k=(0,p^k]\cup(p^{k-1}-p^k,p^{k-1}]\cup\cdots\cup(1-p^k,1]\\
\vdots\\
C=\cap_{k=1}^\infty C_k
$

我们首先挖去(0,1]上中央长度为$(1-2p)$的区间;在剩下的两个区间内,我们继续挖去每一个区间中央长度为$(1-2p)p$的区间,如此不断进行下去,最后得到的集合就是$C$。

给定$\omega \in C$,那么该条件下的$X_k$如何分布呢?

首先考虑$X_1$。容易知道,$P(X_1=0,\omega \in C)=\lambda(C\cap(0,p])$,而$\lambda(C\cap(1-p,1])=P(X_1=1,\omega \in C)$。由集合C的对称性,$\lambda(C\cap(0,p])=\lambda(C\cap(1-p,1])$,因此,$P(X_1=1,\omega \in C)=P(X_1=0,\omega \in C)$。即:在集合C上,$X_1$服从等概率的0-1分布。以此类推,所有的$X_k$在集合$C$上都服从等概率的0-1分布。

强大数定律告诉我们:几乎所有$C$中的实数,与之对应的二进制展开序列均收敛于1/2。而这样的实数构成的子集,属于零测集$A’$的子集,因此也是零测集。这就说明$\lambda(C)=0$。

如果令$p=1/3$,我们就解释了在强大数定律意义下康托三分集测度的含义。

类似地,我们可以进一步解释广义康托集与强大数定律的关系。有兴趣的读者可以作出相应推广。

注:本文已经由COS编辑部整理为PDF(LaTeX)版本,读者可以下载:[button type=”icon” link=”http://cos.name/wp-content/uploads/2010/10/强大数定律与康托三分集.pdf”]强大数定律与康托三分集(PDF,105K)[/button]

比率估计为什么精确

一、比率的方差估计式

比率估计量是抽样技术理论里一大重要估计量,其定义为两个总体总量或总体均值之比。借助适当的辅助变量,比率估计也可以得到主要变量的参数估计

由于通过辅助变量实质上引入了更多的信息,因此有理由猜测比率估计量可能更加精确。但是比率估计的方差和简单估计相比所谓的改进是否确切的存在,即使存在,改进的程度又有多大呢?

记总体大小为$N$,抽样大小为$n$,抽样比例为$f=\frac{n}{N}$,辅助变量的总体值为$ X_1,X_2,…,X_N$,样本值为$x_1,x_2,…,x_n$:主要变量的总体值为$Y_1,Y_2,…,Y_N$,样本值为$y_1,y_2,…,y_n$。教材上常见的一个估计式是:

$ Var(\hat R) \approx \frac{1-f}{\overline{X} n (N-1)}\sum(Y_i – R X_i)^2$
据此,可以给出主要变量相应参数的估计方差。以总体总值为例:

$Var(\hat{Y}) =Var(\hat{R}N \overline{X}) \approx \frac{1-f}{n}\frac{1}{N-1} \sum (Y_i-RX_i)^2 $

注意到上式使用了$\approx$而不是“=”;也就是说是一个近似值。更确切地说,上式估计的只是一个方差下界,因为上式右端实质上是$ Var(\frac{\overline{y}}{\overline{X}}) $;而$\hat{R}=\frac{\overline{y}}{\overline{x}}$。可以看到,比率估计$\hat{R}$方差包括分子、分母两部分波动因素,而估计式中忽略了分母部分的波动,因此得到的方差估计是偏小的。

要使等号严格成立的条件是:

$\overline{x} = \overline{X},a.s.$

在有限总体的情况下,表示辅助变量恒为定值。注意:此时辅助变量已经没有意义了,因为它不能带来更多的信息,比率估计量与简单估计量的精度是完全相同的。

实际应用的时候,为了使方差估计式成立,我们也必须保证:

$\overline{x} \approx \overline{X}$

即样本均值$\overline{x}$总在$\overline{X}$附近波动,且波动范围很小。在这种情况下,辅助变量的意义也很小.

这就是矛盾的所在:比率估计量的方差估计严格成立的场合,也是比率估计量失去应用价值的时候。

二、一个模拟的例子

在样本均值$\overline{x}$波动比较大的时候,比率估计的方差究竟有多大的改进呢?对于这个问题,可以用统计模拟来实现。

我的例子如下:数据来源是人民大学版的《抽样技术》例题4.3,估计33个乡的粮食总产量,抽样得到10个乡粮食产量Y,耕地面积X,村的数量M。Y= (22, 22.8, 30.2, 21.7, 24.3, 31.2, 26, 20.5, 33.8, 23.6),X= (800, 780, 1000, 700, 880, 1100, 850, 800, 1200, 830),M= (15, 18, 26, 14, 20, 28, 21, 19, 31, 17)。

我们可以比较三种方法估计的理论方差:简单估计,以耕地面积作辅助变量的比率估计,以村数量作辅助变量的比率估计。因为总体数据未知,我首先以有放回的抽样模拟一个样本量为33的数据;然后枚举所有可能抽样组合,计算三种估计量。另一方面,对于每种抽样结果,我也采用方差估计式求方差估计值。最后可以将不同方差进行比较。考虑到计算量的问题,仅模拟了样本量为5的情形.

考虑到数据量大,在生成全组合时,采用了字典排序的算法,(可参见http://www.blogjava.net/stme/archive/2007/10/23/94361.html)

#放回抽样,生成总体数据
INDEX = sample(1:10, 33, rep = T)
M = c(15, 18, 26, 14, 20, 28, 21, 19, 31, 17)[INDEX]
Y = c(22, 22.8, 30.2, 21.7, 24.3, 31.2, 26, 20.5,
    33.8, 23.6)[INDEX]
X = c(800, 780, 1000, 700, 880, 1100, 850, 800, 1200,
    830)[INDEX]
#总体总值计算
M0 = sum(M)
X0 = sum(X)
y.simple <- y.m <- y.x <- var.m <- var.x <- NULL

index = c(rep(1, 5), rep(0, 28))

y.simple = c(y.simple, 33 * sum(Y_M * index)/5)
R.m = sum(Y * index)/sum(M * index)
R.x = sum(Y * index)/sum(X * index)
y.m = c(y.m, M0 * R.m)
y.x = c(y.x, X0 * R.x)
var.m = c(var.m, sum((Y[index == 1] - R.m * M[index ==
    1])^2))
var.x = c(var.x, sum((Y[index == 1] - R.x * X[index ==
    1])^2))

i = 1
j = 0
while (prod(index[29:33]) != 1) {
    while (i < 33) {
        if (index[i] == 1 & index[i + 1] == 0) {
            index[i] = 0
            index[i + 1] = 1
            k = sum(index[1:(i - 1)])
            if (k > 0) {
                index[1:k] = 1
                index[(k + 1):i] = 0
            }
            y.simple = c(y.simple, 33 * sum(Y * index)/5)
            R.m = sum(Y * index)/sum(M * index)
            R.x = sum(Y * index)/sum(X * index)
            y.m = c(y.m, M0 * R.m)
            y.x = c(y.x, X0 * R.x)
            var.m = c(var.m, sum((Y[index == 1] - R.m * M[index ==
                1])^2))
            var.x = c(var.x, sum((Y[index == 1] - R.x * X[index ==
                1])^2))
            i = 1
            j = j + 1
            print(j)
        }
        else {
            i = i + 1
        }
    }
}
var.m = var.m/4 * (1/5 - 1/33) * 33 * 33
var.x = var.x/4 * (1/5 - 1/33) * 33 * 33  # simple sampling
> mean(y.simple)
[1] 844.8968
> var(y.simple)
[1] 2678.197
 #ratio with respect to M
> mean(y.m)
[1] 847.828
> var(y.m)
[1] 1156.886
> mean(var.m)
[1] 1111.191
 #ratio with respect to X
> mean(y.x)
[1] 844.9335
> var(y.x)
[1] 221.964
> mean(var.x)
[1] 220.7296

模拟的均值估计结果为:三种方法均值估计为:844.90, 847.83, 844.93;方差为2678.20, 1156.89, 221.96;方差估计的期望为2678.20, 1111.19, 220.73。

这个结果有些出人意料:虽然采用方差估计式得到了低估的结果,但是低估的程度很低,甚至可以忽略不计。也就是说,即使在样本均值波动比较大的场合,比率方差估计的偏误并不大。

这就启示我们对方差估计式的含义重新思考。

三、方差估计式的另一种解释

比率估计量的偏误为:

$E(\frac{\overline y -R \overline x}{\overline x})^2=E( \frac{\overline y -R \overline x}{\overline X})^2 (\frac{\overline X}{\overline x})^2 $

如果假设每次抽样的残差$\overline y -R \overline x$都是一个与 $ \overline x$独立的随机变量,则有:

$E(\frac{\overline y -R \overline x}{\overline x})=E(\frac{\overline y-R \overline x}{\overline X})^2E(\frac{\overline X}{\overline x})^2$

由Jensen不等式,得到

$E(\frac{\overline X}{\overline x})^2 \geq \frac{\overline X ^2}{E ^2 \overline x }=1$

这解释了方差确实存在低估的,而且低估的比例为$ E(\overline X^2 / \overline x^2)$。

采用之前模拟的例子计算这个比例,得到利用耕地面积作辅助变量的抽样方差为121356,但是方差的低估比例仅为1.0035。用此比例修正方差估计,结果为221.51,和真实值221.96几乎相同。

由此可见,即使在辅助变量波动较大,样本两较小,辅助变量抽样均值$\overline x$方差较大的情形,方差低估的比例也可能是很低的,所以采用方差估计式依然可以得到较好的结果。

四、题外话

这个问题给我们的启示:统计学归根结底离不开数学,定量的分析才能给予问题严格的解决。

关于定性和定量的话题,让我想到关于正态分布均值的T检验问题,有的统计学教材上刻意强调了这样一句话:当样本量无限增大的时候,检验结果总是趋向于拒绝。果真如此吗?

上述论断的依据是:随着样本两n的增大,样本均值方差$s^2/n \rightarrow 0$,所以非拒绝域 $(\overline x -t_{\alpha/2}s/\sqrt n,\overline x+t_{\alpha/2}s/\sqrt n)$收缩为一点,因此该拒绝域包括均值$\mu$的可能性为零。

事实上,该论断的谬误是显而易见的:虽然样本方差趋向于零,但是拒绝域包括均值的概率是恒定不变的,这是由拒绝域的构造得到的:

$Pr\{\mu \in (\overline x -t_{\alpha/2}s/\sqrt n,\overline x+t_{\alpha/2}s/ \sqrt n) \} \equiv 1-\alpha$

即使在大样本情形,虽然均值方差趋近于0,非拒绝域的区间非常短,但是只要样本服从原假设下的正态分布,样本均值偏离真实值的可能性也会很小。大数定律告 诉我们,在大样本的情形,样本均值哪怕偏离一丝一毫的概率都为0。因此,哪怕样本均值只有很小的偏离,拒绝零假设也是没有任何问题的。

这就再次说明定量分析的重要。

Hilbert空间视角下的时间序列模型

Hilbert空间说起来和我国古代数学有着一定的渊源。《九章算术》里记载:“勾股术曰:勾股各自乘,并,而开方除之,即弦”。这条著名的勾股定理实质上蕴含了Hilbert空间中对于距离和正交的核心性质。

Hilbert空间特点是通过定义内积$ < \cdot , \cdot > $来导出范数,进而导出距离函数。由内积可以定义正交关系,即若$ <x,y>=0 $则定义$ x , y $在空间中正交。一列彼此正交的元素$ e_j $称为一组正交基,内积的概念给出了任意元素$ x $相对于这组正交基的坐标,即元素$ x $在各个正交基上的投影$ < x , e_j > $.

最常见的Hilbert空间是Euclidean空间$ R^n $,即是线性代数研究的范畴。平面几何研究的是$ R^2 $,“勾股定理”反映的性质就是属于这个空间内的。“勾股定理”用Hilbert空间中的概念表达即为:
给定$ R^2 $中的一组标准正交基${e_1 , e_2} $,则任意点$ x $关于0元素的距离$ \parallel x \parallel $满足:$ \parallel x \parallel^2=\mid<x , e_1>\mid^2 + \mid <x , e_2>\mid^2 $.

这条定理的一般形式正是Hilbert空间中著名的Parseval公式,前提是正交系$ \{e_j\} $是完备的。在此前提下,Parseval公式成立:

$ \parallel x \parallel^2 =\sum_j \mid < x , e_j >\mid^2 $

Parseval公式揭示了Hilbert空间的优良特点,即任一点的位置可以相对于一组完备正交基完全定义,而且这种定义完全保留了点与点之间的距离关系。这就给计算带来极大的便利。

时序模型研究的是一类特殊的Hilbert空间:$ \mathscr{L}^2\{\Omega,\mathscr{F},P\} $,即一个定义在概率空间$ \{\Omega,\mathscr{F},P\} $上的平方可积函数集合,其中的内积定义为:

$ <X,Y>=E(XY)=\int_\Omega xydP(x,y) $

在提出这个定义以前,时序中经常提到的“均方收敛”的概念难以理解:一个随机变量“均方收敛”的结果是另一随机变量,也就是说这个序列收敛的结果取值仍然是不确定的;这和数学分析中研究的收敛范畴大相径庭。实际上,“收敛”的确切含义是指点与点之间距离趋向于0,但是对于这种“距离”的定义在不同的空间中可以是不同的。数学分析中的距离由绝对值定义,而这里的“距离”则由内积定义。均方收敛实际上与上述距离的定义完全一致,即:

$ \parallel X_m-X_n \parallel =(E(X_m-X_n)^2)^{\frac{1}{2}} $

再来看$ \mathscr{L}^2\{\Omega,\mathscr{F},P\} $的正交基的形式。其实最简单的正交基就是白噪声序列$ \{a_t\} \sim i.i.d(0,\sigma^2) $,因为

$ <a_i,a_j> = E a_ia_j=\delta_{ij}\sigma^2 $

为了应用之前的Parseval公式,严格的说需要证明$ \{a_t\} $是完备的。但如果我们只研究这组基生成的线性子空间$ \mathscr{H}:=\overline{sp}{a_j} $,那么显然在$ \mathscr{H} $上这组正交基是完备的。这个空间对线性时序模型已经完全足够了。

对于平稳的ARMA序列,$ X_t $的信息是由t时刻以前的白噪声序列决定的,写成数学形式即

$ X_t \in\mathscr{H}_t:=\overline{sp}{a_j,j \leq t} $

这反映了序列的鞅性质。更一般地,$ X_t $在Hilbert子空间$ \overline{sp}{a_j,j \leq t} $有固定的表达,这就是模型的$ MA(\infty) $表示。但是在实际情况中对$ \{X_j,j < t\} $的观测更为直接,因此我们也关心$ X_t $ 在$ \overline{sp}{X_j,j < t} $上的投影表达,即模型的$ AR(\infty) $表达形式。注意到$ \overline{sp}{X_j,j<t}=\overline{sp}{a_j,j<t} $,因此这两种表达形式本质是一致的,只是前者子空间的生成元不是正交基。

线性时序模型ARMA的预测问题可以归结为求解$ X_t $在空间$ \mathscr{H}_{t-1} $的投影,即在Hilbert子空间$ \mathscr{H}_{t-1} $中的最佳逼近元。最佳逼近元的泛函概念是在$ \mathscr{H}_{t-1} $中寻找使得距离$ \parallel X-\hat{X} \parallel $最小的估计值$ \hat{X} $,从范数的定义很容易看出,这个概念和统计上“均方误差最小”的概念是一致的。注意到在线性回归模型中,我们求解的目标也是$ \parallel y-\hat{y} \parallel $最小,但是线性回归模型中范数$ \parallel \cdot \parallel $是定义在Euclidean空间$ R^n $上的。最小二乘模型中的$ R^n $也是一种Hilbert空间,而$ R^n $和$ \mathscr{L}^2\{\Omega,\mathscr{F},P\} $的区别也说明了ARMA回归和线性回归的差别。线性回归模型的最小二乘法完全源于$ R^n $的性质,但是对于ARMA 模型来说,$ \parallel X-\hat{X} \parallel $的计算是几乎不可能的,因此首先都是通过构造一组由$ \{a_j,1 \leq j \leq t-1\} $作为正交基生成的Euclidean空间$ R^{t-1} $去对$ \mathscr{H}_{t-1} $作近似,用$ R^{t-1} $中的范数近似$ \mathscr{H}_{t-1} $的范数。

这种做法背后的假设是,随机白噪声序列的观测值$ \{a_j,1 \leq j \leq t-1\} $在距离上典型的“代表意义”,即以观测值为基生成的欧式空间$ R^{t-1} $可以反映$ \mathscr{H}_{t-1} $空间中的距离信息,这个思想和极大似然估计是相似的,即把观测值都是具有典型性代表性的。有了这个空间的变换,就把原先非欧空间中的优化问题转化为欧式空间中的问题,解决起来也就容易多了。从算法上看,回归模型和ARMA模型都是在欧式空间中进行,方法十分类似;但是二者本质上是存在差别的,回归模型本身就是欧式空间中的问题,而把ARMA模型放到欧式空间中进行求解只是一种简化。当对统计量极限分布进行研究时,ARMA模型必须重新考虑在原有空间$ \mathscr{L}^2\{\Omega,\mathscr{F},P\} $中距离的定义,这个过程就远比线性回归模型要复杂了。

时间序列的谱分析实质上是在另一种Hilbert空间$ \mathscr{L}^2(F) $的视角下进行研究。Fourier变换定义了一个映射$ T : \mathscr{L}^2\{\Omega,\mathscr{F},P\}\rightarrow \mathscr{L}^2(F) $,该映射满足

$ TX_t=e^{it} $

这两个空间是通过一个$ \mathscr{L}^2\{\Omega,\mathscr{F},P\} $上的正交增量过程$ Z(v) $发生联系的,其中

$ X_t=\int_{(-\pi,\pi]}e^{itv}dZ(v) $

定义$ F(v):=var(Z(v)-Z(-\pi)) $,则在空间$ \mathscr{L}^2(F) $上,利用Ito积分的有关性质得到:

$ <e^{i(t+h)},e^{it}>=\int e^{i(t+h)v-itv}dF(v)=\int e^{ihv}dF(v)=\gamma(h) $

因此在空间$ \mathscr{L}^2(F) $上可以直观地对自相关函数进行描述。函数$ F(v) $的导数$ f(v) $就是谱函数。对有限序列$ \{X_t,1 \leq t \leq n \} $作Fourier变换得到的序列$ f_n(\omega_j) $,就是$ f(v) $在各个频率上的估计。可以证明对固定的时序模型,该映射$ T $对应的正交增量过程$ Z(v) $以概率1恒定;所以相应的谱函数$ f(v) $是唯一的。

以上提到的两种Hilbert空间为时序模型提供了时域和频域两种不同的视角,也在泛函领域奠定了两种分析方法的的理论基础。在对时序模型统计量的收敛性进行分析时,以上两种视角是必不可少的。

注:本文已经由COS编辑部整理为PDF(LaTeX)版本,读者可以下载:
http://cos.name/wp-content/uploads/2009/03/Hilbert空间视角下的时间序列模型.pdf