Back
Featured image of post 学习笔记|数理统计

学习笔记|数理统计

数理统计、贝叶斯统计、高维统计相关知识

课程介绍

数理统计部分大多数是统计研究院陳鄰安老师数理统计,陈老师年纪比较大了,但是是很有激情的老师,他的课值得一听。

参考:

  1. 网站:统计学1
  2. 网站:统计学2
  3. youtube课程:交大數理統計學
  4. youtube课程:陳鄰安老師 - 高等統計學
  5. 教材:Introduction to Mathematical Statistics ,Hogg
  6. 教材:统计推断,Casella
  7. 教材:应用随机过程,Ross
  8. 网页:统计计算
  9. 讲义:Markov chain Monte Carlo methods

陈老师的高等统计学没有讲义,而且讲的过于简单,本人使用了Casella的统计推断进行学习,并且配合一些其他的统计课程进行的补充。

Fundamentals of probability theory

Random Variable

Define 1(random variable). Let $(\Omega, \mathcal{F}, P)$ be a probability space and $(E, \mathcal{E})$ a measurable space. Then an $(E, \mathcal{E})$-valued random variable is a measurable function $X: \Omega \rightarrow E$.
Remark.

  1. 随机变量是从样本空间到$\mathbb{R}$的映射,概率测度是从$\sigma$-域到$[0,1]$的映射,并且满足整个样本空间的测度等于1。在实际中我们一般不会考虑原来的样本空间,因为一般考虑由随机变量诱导出来的概率测度,它是由博雷尔集($\mathcal{B}$)到$\mathbb{R}$的映射。
  2. $\forall B\in \mathcal{B},X^{-1}(B)\in \mathcal{F},\text{where} X^{-1}(B) = \{\omega \in \Omega\mid X(\omega)\in B\}$,which means that, for every subset $B \in \mathcal{E}$, its preimage is $\mathcal{F}$-measurable; $X^{-1}(B) \in \mathcal{F}$, where $X^{-1}(B)=\{\omega: X(\omega) \in B\}$. This definition enables us to measure any subset $B \in \mathcal{E}$ in the target space by looking at its preimage, which by assumption is measurable.
  3. 由随机变量很容易定义累积分布函数(CDF)
$$ F(x) = P(X(\omega) \leq x) = P(\{\omega\in\Omega:X(\omega)\in (-\infty,x]\}) $$

Def 2(m.g.f.). Let $X$ be a random variable with CDF $F_X$. The moment generating function (mgf) of $X$ (or $F_X$ ), denoted by $M_X(t)$, is $$ M_X(t)=\mathrm{E}\left[e^{t X}\right] $$ provided this expectation exists for $t$ in some open neighborhood of 0 . That is, there is an $h>0$ such that for all $t$ in $-h<t<h, \mathrm{E}\left[e^{t X}\right]$ exists. If the expectation does not exist in an open neighborhood of 0 , we say that the moment generating function does not exist.
In other words, the moment-generating function of $X$ is the expectation of the random variable $e^{t X}$. More generally, when $\mathbf{X}=\left(X_1, \ldots, X_n\right)^{\mathrm{T}}$, an $n$-dimensional random vector, and $\mathbf{t}$ is a fixed vector, one uses $\mathbf{t} \cdot \mathbf{X}=\mathbf{t}^{\mathrm{T}} \mathbf{X}$ instead of $t X$ : $$ M_{\mathbf{X}}(\mathbf{t}):=\mathrm{E}\left(e^{\mathbf{t}^{\mathrm{T}} \mathbf{X}}\right) $$ $M_X(0)$ always exists and is equal to 1 .
Note. 几种常见分布的 m.g.f: $$ \begin{cases}X \sim N\left(\mu, \sigma^2\right) & , M_x(t)=e^{\mu t+\frac{\sigma^2}{2} t^2}, \forall t \in R \\ X \sim \operatorname{Gamma}(\alpha, \beta) & , M_x(t)=(1-\beta t)^{-\alpha}, t<\frac{1}{\beta} \\ X \sim b(n, p) & , M_x(t)=\left(1-p+p e^t\right)^n, \forall t \in R \\ X \sim \operatorname{Poisson}(\lambda) & , M_x(t)=e^{\lambda\left(e^t-1\right)}, \forall t \in R\end{cases} $$

Lemma 3. 若$X$为r.v.,令g为$\mathbb{R}\rightarrow\mathbb{R}$的映射,则$Y=g(X)$同样为r.v.,且可以确定其分布。
Remark.
两种求得$g(X)$的pdf的方法:

  • Distribution method:

    Suppose that X is a continuous r.v.. Let $Y=g(X)$ The d.f(distribution function) of Y is $$ G(y)=P(Y \leq y)=P(g(X) \leq y) $$ If G is differentiable then the p.d.f. of $Y=g(X)$ is $g(y)=G^{\prime}(y)$.

  • mgf method :(moment generating function) $$ \begin{aligned} E\left[e^{t x}\right]= \begin{cases}\sum e^{t x} f(x) & \text { (discrete) } \\ \int_{-\infty}^{\infty} e^{t x} f(x) d x & \text { (continuous) }\end{cases} \end{aligned} $$

Thm 4. m.g.f. $M_x(t)$ and its distribution (p.d.f. or d.f.) forms a $1-1$ functions.

Thm 5. $\quad M_X^{(k)}(0)=E\left[X^k\right], k=1,2, \cdots$

Probability distribution

Def 1. Gamma function $\Gamma(\alpha)=\int_0^{\infty} x^{\alpha-1} e^{-x} d x$
Properties:

  1. $\Gamma(\alpha)=(\alpha-1) \Gamma(\alpha-1)$, if $\alpha>1$
  2. $\Gamma(1)=\int_0^{\infty} e^{-x} d x=1$
  3. $\Gamma(n)=(n-1) \Gamma(n-1)=\cdots=(n-1)(n-2) \cdots 1 \cdot \Gamma(1)=(n-1)$ !
  4. $\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi}$

Def 2. We say that $X$ has a Gamma distribution if it has p.d.f $$ f(x)=\frac{1}{\Gamma(\alpha) \beta^\alpha} x^{\alpha-1} e^{-\frac{x}{\beta}}, x>0, \text { for some } \alpha>0, \beta>0 $$ We denote by $X \sim \operatorname{Gamma}(\alpha, \beta)$.
Note. If $X$ has Gamma distribution with $\beta=2$ and $\alpha=\frac{r}{2}$, we say that $X$ has a chi-square distribution with degrees of freedom r. The p.d.f is $$ f(x)=\frac{1}{\Gamma\left(\frac{r}{2}\right) 2^{\frac{r}{2}}} x^{\frac{r}{2}-1} e^{-\frac{x}{2}}, x>0 $$ We denote by $X \sim \chi^2(r)$.

Ex 1. If $X \sim \operatorname{Gamma}(\alpha, \beta)$, m.g.f of $X$ is

$$ \begin{aligned} M_X(t) & =E\left[e^{t X}\right]=\int_0^{\infty} e^{t x} \frac{1}{\Gamma(\alpha) \beta^\alpha} x^{\alpha-1} e^{-\frac{x}{\beta}} d x=\int_0^{\infty} \frac{1}{\Gamma(\alpha) \beta^\alpha} x^{\alpha-1} e^{-\frac{(1-\beta t) x}{\beta}} d x \\ & =(1-\beta t)^{-\alpha} \int_0^{\infty} \frac{1}{\Gamma(\alpha)\left(\frac{\beta}{1-\beta t}\right)^2} x^{\alpha-1} e^{-\frac{x}{1-\beta t}} d x \\ & =(1-\beta t)^{-\alpha}, t<\frac{1}{\beta} \\ & \because \frac{\beta}{1-\beta t}>0 \Rightarrow 1-\beta t>0 \Rightarrow t<\frac{1}{\beta} \\ M_X^{\prime}(t) & =\alpha(1-\beta t)^{-\alpha-1} \beta \\ & \Rightarrow \operatorname{Mean} \mu=E[X]=M_X^{\prime}(0)=\alpha \beta \\ M_X^{\prime \prime}(t) & =\alpha(\alpha+1)(1-\beta t)^{-\alpha-2} \beta^2 \\ & \Rightarrow E\left[X^2\right]=M_X^{\prime \prime}(0)=\alpha(\alpha+1) \beta^2 \\ & \Rightarrow \text { Variance } \sigma^2=M_X^{\prime \prime}(0)-\left(M_X^{\prime}(0)\right)^2=\alpha(\alpha+1) \beta^2-(\alpha \beta)^2=\alpha \beta^2 \end{aligned} $$

Random Voctor

Definition 1. Suppose that the random variables $X$ and $Y$ are defined to assume values in $\displaystyle I\subseteq \mathbb {R}$. We say that $n$ random variables $X_1, \ldots, X_n$ are i.i.d(Independent and identically distributed random variables). If and only if

  • $F_{X_1}(x)=F_{X_k}(x)$ , $\forall k \in{1, \ldots, n}$ and $\forall x \in I$
  • $F_{X_1, \ldots, X_n}\left(x_1, \ldots, x_n\right)=F_{X_1}\left(x_1\right) \cdot \ldots \cdot F_{X_n}\left(x_n\right)$ , $\forall x_1, \ldots, x_n \in I$

where $F_{X_1, \ldots, X_n}\left(x_1, \ldots, x_n\right)=\mathrm{P}\left(X_1 \leq x_1 \wedge \ldots \wedge X_n \leq x_n\right)$ denotes the joint cumulative distribution function of $X_1, \ldots, X_n$

Thm 2(Theory for change variables). $$ P\left(\left(\begin{array}{c} x_1 \\ \vdots \\ x_n \end{array}\right) \in A\right)=\int \cdots \int f_{X_1, \ldots, X_n}\left(x_1, \ldots, x_n\right) d x_1 \cdots d x_n $$ Let $y_1=g_1\left(x_1, \ldots, x_n\right), \cdots, y_n=g_n\left(x_1, \ldots, x_n\right)$ be a $1-1$ function with inverse $x_1=w_1\left(y_1, \ldots, y_n\right), \cdots, x_n=w_n\left(y_1, \ldots, y_n\right)$ and Jacobian $$ J=\left|\begin{array}{ccr} \frac{\partial x_1}{\partial y_1} & \cdots & \frac{\partial x_1}{\partial y_n} \\ \vdots & & \vdots \\ \frac{\partial x_n}{\partial y_1} & \cdots & \frac{\partial x_n}{\partial y_n} \end{array}\right| $$ Then $$ \begin{aligned} & \int \cdots \int f_{X_1, \ldots, X_n}\left(x_1, \ldots, x_n\right) d x_1 \cdots d x_n \\ & =\int \cdots \int f_{X_1, \ldots, X_n}\left(w_1\left(y_1, \ldots, y_n\right), \ldots, w_n\left(y_1, \ldots, y_n\right)\right)|J| d y_1 \cdots d y_n \end{aligned} $$ Hence, joint p.d.f. of $Y_1, \cdots, Y_n$ is $$ f_{Y_1, \ldots, Y_n}\left(y_1, \ldots, y_n\right)=f_{X_1, \ldots, X_n}\left(w_1, \ldots, w_n\right)|J| $$
Remark. 可以用该方法构造从$n$个随机变量的联合分布到另外$n$个随机变量联合分布的映射(注意映射到的值域), 然后再通过积分的方式获得边缘联合分布。

Statistical basis

Def 1(随机样本). If a sequence of r.v.’s $X_1,\cdots,X_n$ are i.i.d.,then they are called a random sample.

Def 2(统计量和参数). Any function $g(X_1,\cdots,X_n)$ of a random sample $X_1,\cdots,X_n$ which is not dependent on a parameter $\theta$ is called a statistic.
Note. If $X$ is a random sample with p.d.f. $f(x;\theta)$, where $\theta$ is an unknown constant, then $\theta$ is called parameter.统计量是一个可观察的随机变量,这是它与统计参数的直接区别 ,统计参数包含一个不可观察的随机变量从而可以描述统计总体。 统计参数必须要在整体数据都可被观察的时候才能计算,例如,一个完美的人口普查。

为了介绍 Thm 8,我们介绍以下定义和定理作为铺垫,这些定理同样是很实用的定理。

Lemma 2(独立的等价判定). $X_1$ and $X_2$ are independent if and only if $$ M_{X_1, X_2}\left(t_1, t_2\right)=M_{X_1}\left(t_1\right) M_{X_2}\left(t_2\right), \forall t_1, t_2 . $$

Thm 3. 关于独立的两个定理:

  1. If $\left(X_1, \ldots, X_n\right)$ and $\left(Y_1, \ldots, Y_m\right)$ are independent(that is:$f(X_1, \ldots, X_n,Y_1, \ldots, Y_m) = f(X_1, \ldots, X_n)f(Y_1, \ldots, Y_m)$), then $g\left(X_1, \ldots, X_n\right)$ and $h\left(Y_1, \ldots, Y_m\right)$ are also independent.
  2. If $X, Y$ are independent, then $$ \mathrm{E}[g(X) h(Y)]=\mathrm{E}[g(X)] \mathrm{E}[h(Y)] $$

Thm 4. Joint m.g.f. $M_{X,Y}(0,t) = M_Y(t)$.

Thm 5. $X \sim \chi^2(r_1)$,$Y \sim \chi^2(r_2)$,if $X\perp Y$,then $X+Y\sim \chi^2(r_1+r_2)$.

Thm 6. $Z\sim N(0,1)$,then $Z^2\sim \chi^2(1)$.

Def 7(样本均值,样本方差). 一组随机样本 $X_1,\ldots,X_n$,样本均值为 $\bar{X} = \frac{1}{n}\sum_{i=1}^{n} X_i$,样本方差为 $S^2 = \frac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar{X})^2$.二者都是统计量。

Thm 8. If $\left(X_1, \ldots, X_n\right)$ is a random sample from $N\left(\mu, \sigma^2\right)$, then
(a) $\bar{X} \sim N\left(\mu, \frac{\sigma^2}{n}\right)$
(b) $\bar{X}$ and $S^2$ are independent.
(c) $\frac{(n-1) S^2}{\sigma^2} \sim \chi^2(n-1)$
Remark.

  1. 样本均值不等于均值,样本均值的期望才等于均值,注意区分以下两者的关系。 $$ \begin{cases} \sum_i \left( \frac{X_i-\mu}{\sigma}\right)^2\sim\chi^2(n) \\ \frac{(n-1) S^2}{\sigma^2} = \sum_i \left( \frac{X_i-\bar{X}}{\sigma}\right)^2 \sim \chi^2(n-1) \end{cases} $$
  2. Thm 8 的(b)是十分难得的,因为二者都可以看作 $X_i$ 的函数,而存在独立性, 归根结底是因为 $\bar{X}$ 和 $\left(X_1-\bar{X},\ldots,X_n-\bar{X}\right)$ 独立, 即 $$M_{\bar{X}}(t)M_{X_1-\bar{X},\ldots,X_n-\bar{X}}(t_1,\ldots,t_n) = M_{\bar{X},X_1-\bar{X},\ldots,X_n-\bar{X}}(t,t_1,\ldots,t_n)$$.

Statistical Inference - Point Estimation

Problem in statistics:
A random variables $X$ with p.d.f. of the form $f(x, \theta)$ where function $f$ is known but parameter $\theta$ is unknown. We want to gain knowledge about $\theta$.
What we have for inference:
There is a random sample $X_1, \ldots, X_n$ from $f(x, \theta)$. $$ \text { Statistical inferences }\left\{\begin{array}{l} \text { Estimation }\left\{\begin{array}{l} \text { Point estimation: } \hat{\theta}=\hat{\theta}\left(X_1, \ldots, X_n\right) \\ \text { Interval estimation: } \\ \text { Find statistics } T_1=t_1\left(X_1, \ldots, X_n\right), T_2=t_2\left(X_1, \ldots\right) \\ \text { such that } 1-\alpha=P\left(T_1 \leq \theta \leq T_2\right) \end{array}\right. \\ \text { Hypothesis testing: } H_0: \theta=\theta_0 \text { or } H_0: \theta \geq \theta_0 . \\ \text { Want to find a rule to decide if we accept or reject } H_0 . \end{array}\right. $$

Unbiased and consistent

Def 1(估计量). We call a statistic $\hat{\theta}=\hat{\theta}\left(X_1, \ldots, X_n\right)$ an estimator of parameter $\theta$ if it is used to estimate $\theta$. If $X_1=x_1, \ldots, X_n=x_n$ are observed, then $\hat{\theta}=\hat{\theta}\left(x_1, \ldots, x_n\right)$ is called an estimate of $\theta$.

Def 2(无偏估计). We call an estimator $\theta$ unbiased for $\theta$ if it satisfies $$ \begin{gathered} E_\theta\left(\hat{\theta}\left(X_1, \ldots, X_n\right)\right)=\theta, \forall \theta \\ \mathrm{E}_\theta\left(\hat{\theta}\left(X_1, \ldots, X_n\right)\right)=\left\{\begin{array}{l} \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} \hat{\theta}\left(x_1, \ldots, x_n\right) f\left(x_1, \ldots, x_n, \theta\right) d x_1 \cdots d x_n \\ \int_{-\infty}^{\infty} \theta^* f_{\hat{\theta}}\left(\theta^*\right) d \theta^* \text { where } \hat{\theta}=\hat{\theta}\left(X_1, \ldots, X_n\right) \text { is a r.v. with pdf } f_{\hat{\theta}}\left(\theta^*\right) \end{array}\right. \end{gathered} $$
Remark. $\bar{X}$ 和 $S^2$ 都是无偏估计量。

既然是点估计,我们有必要选择一种度量来衡量点估计是否准确(不可以选择概率作为度量,因为对于连续分布来说,统计量等于任意一点的概率为0,没有意义)。

Def 3(依概率收敛). We say that $X_n$ converges to $X$,(a r.v. or a constant) in probability if for $\epsilon>0$, $$ P\left(\left|X_n-X\right|>\epsilon\right) \longrightarrow 0, \text { as } n \longrightarrow \infty $$ In this case, we denote $X_n \xrightarrow{P} X$.
Remark. 该定理的意思是,$\{\omega\in\Omega:|X_n(\omega)-X(\omega)|<\epsilon\}$这个事件的概率测度收敛到1.不加证明的说明一个定理:对于连续函数 $g$,若$X_n \xrightarrow{P} X$,则$g(X_n) \xrightarrow{P} g(X)$.

下面引入两种容易判断随机变量依概率收敛的定理,在引入定理之前,回顾一下集中不等式。

Thm 4(马尔可夫不等式). If $X$ is a nonnegative random variable and $a>0$, then the probability that $X$ is at least $a$ is at most the expectation of $X$ divided by a: $$ \mathrm{P}(X \geq a) \leq \frac{\mathrm{E}(X)}{a} $$ When $\mathrm{E}(X)>0$, we can take $a=\tilde{a} \cdot \mathrm{E}(X)$ for $\tilde{a}>0$ to rewrite the previous inequality as $$ \mathrm{P}(X \geq \tilde{a} \cdot \mathrm{E}(X)) \leq \frac{1}{\tilde{a}} $$

Thm 5(判断依概率收敛). If $E\left(X_n\right)=a$ or $E\left(X_n\right) \longrightarrow a$ and $\operatorname{Var}\left(X_n\right) \longrightarrow 0$ when $n\rightarrow\infty$, then $X_n \xrightarrow{P} a$.
Proof. $$ \begin{aligned} \mathrm{E}\left[\left(X_n-a\right)^2\right] & =\mathrm{E}\left[\left(X_n-\mathrm{E}\left(X_n\right)+\mathrm{E}\left(X_n\right)-a\right)^2\right] \\ & =\mathrm{E}\left[\left(X_n-\mathrm{E}\left(X_n\right)\right)^2\right]+\mathrm{E}\left[\left(\mathrm{E}\left(X_n\right)-a\right)^2\right]+2 \mathrm{E}\left[\left(X_n-\mathrm{E}\left(X_n\right)\right)\left(\mathrm{E}\left(X_n\right)-a\right)\right] \\ & =\operatorname{Var}\left(X_n\right)+\mathrm{E}\left(\left(X_n\right)-a\right)^2 \end{aligned} $$ By Chebyshev’s Inequality(可以由马尔可夫不等式得到) : $$ P\left(\left|X_n-X\right| \geq \epsilon\right) \leq \frac{\mathrm{E}\left(X_n-X\right)^2}{\epsilon^2} \\ \text { or } P\left(\left|X_n-\mu\right| \geq k \sigma\right) \leq \frac{1}{k^2} $$ For $\epsilon>0$, $$ \begin{aligned} & 0 \leq P\left(\left|X_n-a\right|>\epsilon\right)=P\left(\left(X_n-a\right)^2>\epsilon^2\right) \\ & \leq \frac{\mathrm{E}\left(X_n-a\right)^2}{\epsilon^2}=\frac{\operatorname{Var}\left(X_n\right)+\left(\mathrm{E}\left(X_n\right)-a\right)^2}{\epsilon^2} \longrightarrow 0 \text { as } n \longrightarrow \infty . (\text{条件})\\ & \Rightarrow P\left(\left|X_n-a\right|>\epsilon\right) \longrightarrow 0, \text { as } n \longrightarrow \infty . \Rightarrow X_n \xrightarrow{P} a . \end{aligned} $$ Q.E.D.
Remark. 上述定理中 $X_n$ 是统计量,我们同样可以用 $\hat{\theta}$ 将其表达,通过该定理,可以容易的判断统计量是否依概率收敛于参数。

Thm 6(Weak Law of Large Numbers,WLLN). If $X_1, \ldots, X_n$ is a random sample with mean $\mu$ and finite variance $\sigma^2$, then $\bar{X} \xrightarrow{P} \mu$.
Proof. $$ \mathrm{E}(\bar{X})=\mu, \operatorname{Var}(\bar{X})=\frac{\sigma^2}{n} \longrightarrow 0 \text { as } n \longrightarrow \infty . \Rightarrow \bar{X} \xrightarrow{P} \mu $$ Q.E.D.
Remark. 上述定理中的 $X_1, \ldots, X_n$同样可以写为 $X_1^k, \ldots, X_n^k$,则可以得到,若 $Var(X^k)$ 有限,则 $\bar{X^k} \xrightarrow{P} \mathbb{E}X^k$.

对于满足依概率收敛的统计量,有以下定义:

Def 7(一致估计量). We sat that $\hat{\theta}$ is a consistent estimator of $\theta$ if $\hat{\theta} \xrightarrow{P} \theta$ .
Remark. $\bar{X}和S^2$分别是 $\mu$ 和 $\sigma^2$一致估计量.

Def 8(Moments): Let $X$ be a random variable having a p.d.f. $f(x, \theta)$, the population $k_{t h}$ moment is defined by $$ \mathrm{E}_\theta\left(X^k\right)= \begin{cases}\sum_{\text {all } x} x^k f(x, \theta) & , \text { discrete } \\ \int_{-\infty}^{\infty} x^k f(x, \theta) d x & , \text { continuous }\end{cases} $$ The sample $k_{t h}$ moment is defined by $\frac{1}{n} \sum_{i=1}^n X_i{ }^k$.
Remark. 假设方差有限,样本原点矩是原点矩的无偏估计量和一致性估计量 (WLLN)。

Point estimation method

Method 1:矩估计,Method of Moments, MoM
使用样本矩估计原点矩,然后反解出参数 $\theta$ 作为 $\hat{\theta}$。
If $\theta=\left(\theta_1, \ldots, \theta_k\right)$ is $k$-variate, the method of moment estimator $\left(\hat{\theta_1}, \ldots, \hat{\theta_k}\right)$ solves $\theta_1, \ldots, \theta_k$ for $$ \frac{1}{n} \sum_{i=1}^n X_i{ }^j=E_{\theta_1, \ldots, \theta_k}\left(X^j\right), j=1, \ldots, k $$ Remark. 矩估计方法可以保证估计量相对于k阶原点矩是无偏且一致的,但无法保证中心矩的无偏性和一致性。


Method 2:,最大似然估计,Maximum Likelihood Estimator
Let $X_1, \ldots, X_n$ be a random sample with p.d.f. $f(x, \theta)$. The joint p.d.f. of $X_1, \ldots, X_n$ is $$ f\left(x_1, \ldots, x_n, \theta\right)=\prod_{i=1}^n f\left(x_i, \theta\right), x_i \in R, i=1, \ldots, n $$ Let $\Theta$ be the space of possible values of $\theta$. We call $\Theta$ the parameter space.

Def 1(似然函数). The likelihood function of a random sample is defined as its joint p.d.f. as $$ L(\theta)=L\left(\theta, x_1, \ldots, x_n\right)=f\left(x_1, \ldots, x_n, \theta\right), \theta \in \Theta $$ which is considered as a function of $\theta$. Remark. 这里穿插一点贝叶斯学派和频率学派不同的一点,频率学派认为 $\theta$ 是一个确定的数值,尽管难以得到,但是这个数值是确定的,而分布一般写成 $F(X_1,\cdots,X_N,\theta)$,这是为了强调 $\theta$ 作为函数的参数这一点,但是按照严格的概率论写法,括号内应为随机变量;贝叶斯学派认为 $\theta$ 本身就是服从一定分布的随机变量,一般将似然函数写为 $L(\theta) = f(X_1,X_2,\cdots,X_n\mid \theta)$,写法的不同反应了对参数空间认知的根本差异。

Def 2(似然估计量). Let $\hat{\theta}=\hat{\theta}\left(X_1, \ldots, X_n\right)$ be any value of $\theta$ that maximizes $L\left(\theta, X_1, \ldots, X_n\right)$. Then we call $\hat{\theta}=\hat{\theta}\left(X_1, \ldots, X_n\right)$ the maximum likelihood estimator (m.l.e) of $\theta$. When $X_1=x_1, \ldots, X_n=x_n$ is observed, we call $\hat{\theta}=\hat{\theta}\left(x_1, \ldots, x_n\right)$ the maximum likelihood estimate of $\theta$.
Remark. 即似然估计量为 $$ \hat{\theta} = \argmax_{\theta\in\Theta}L(\theta) $$ 理解最大似然估计有很多种方法,一种最直观的方法就是使得当前采样事件发生概率最大的参数就是最大似然估计量需要估计的参数,可以证明,最大似然估计量可以等价表示为(或者具有如下性质) $$ \hat{\theta} =\argmin_{\theta\in\Theta} D_{k L}\left(f_\theta(x) \| f_{\hat{\theta}}(x)\right) $$ 其中,$f_\theta(x) $ 和 $f_{\hat{\theta}}(x)$ 分别表示真实分布和估计分布。
Remark 2. 凡所发生,皆应发生。最大似然估计就是找到使得当前事件发生概率最大的估计量。

在实际的分析中,注意 $ln L(\theta)$ 并不会影响最大似然的解,使用对数函数可以极大的简化运算。另外,这里引入顺序统计量的概念和一个定理,可能会在求mle的过程中遇到:

Def 3(Order statistics). Let $\left(X_1, \ldots, X_n\right)$ be a random sample with d.f. F and p.d.f. f. Let $\left(Y_1, \ldots, Y_n\right)$ be a permutation $\left(X_1, \ldots, X_n\right)$ such that $Y_1 \leq Y_2 \leq \cdots Y_n$. Then we call $\left(Y_1, \ldots, Y_n\right)$ the order statistic of $\left(X_1, \ldots, X_n\right)$ where $Y_1$ is the first (smallest) order statistic, $Y_2$ is the second order statistic,…, $Y_n$ is the largest order statistic.
Remark. $Y_1 = \min\{X_1, \ldots, X_n\}$,$Y_n = \max\{X_1, \ldots, X_n\}$;顺序统计量一般不相互独立,也不同分布。

Thm 4(两个顺序统计量的pdf). Let $\left(X_1, \ldots, X_n\right)$ be a random sample from a “continuous distribution” with p.d.f. $f(x)$ and d.f $F(x)$. Then the p.d.f. of $Y_n=\max \left\{X_1, \ldots, X_n\right\}$ is $$ g_n(y)=n(F(y))^{n-1} f(y) $$ and the p.d.f. of $Y_1=\min \left\{X_1, \ldots, X_n\right\}$ is $$ g_1(y)=n(1-F(y))^{n-1} f(y) $$

Thm 5(Invariance property of a mle). If $h(\cdot)$ is a $1-1$ function,and $\hat{\theta}$ is a mle,then $$ \hat{h(\theta)} = h(\hat{\theta}) $$

UMVUE

现在我们有了两种点估计的方法,矩估计方法和最大似然估计方法,我们可以得到一些估计量,并且判断它们的无偏性和一致性,但是如何进一步评判一个统计量的好坏仍然是现在不知道的,接下来我们要定义一种对统计量好坏的度量,先启发的想,方差代表了统计量测度集中的程度,对于点估计而言,我们自然是希望方差越小越好:

$$ \begin{aligned} MSE(\hat{\theta})& =\mathrm{E}_\theta\left[(\hat{\theta}-\theta)^2\right] \\\ & =\mathrm{E}_\theta\left[\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]+\mathrm{E}_\theta[\hat{\theta}]-\theta\right)^2\right] \\\ & =\mathrm{E}_\theta\left[\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)^2+2\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)+\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)^2\right] \\\ & =\mathrm{E}_\theta\left[\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)^2\right]+\mathrm{E}_\theta\left[2\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)\right]+\mathrm{E}_\theta\left[\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)^2\right] \\\ & =\mathrm{E}_\theta\left[\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)^2\right]+2\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right) \mathrm{E}_\theta\left[\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right]+\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)^2 \\\ & =\mathrm{E}_\theta\left[\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)^2\right]+2\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)\left(\mathrm{E}_\theta[\hat{\theta}]-\mathrm{E}_\theta[\hat{\theta}]\right)+\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)^2 \\\ & =\mathrm{E}_\theta\left[\left(\hat{\theta}-\mathrm{E}_\theta[\hat{\theta}]\right)^2\right]+\left(\mathrm{E}_\theta[\hat{\theta}]-\theta\right)^2 \\\ & =\operatorname{Var}_\theta(\hat{\theta})+\operatorname{Bias}_\theta(\hat{\theta}, \theta)^2 \end{aligned} $$

于是有以下定义:

Def 1(UMVUE). An unbiased estimator $\hat{\theta}=\hat{\theta}\left(X_1, \ldots, X_n\right)$ is called a uniformly minimum variance unbiased estimator (UMVUE) or best estimator if for any unbiased estimator $\hat{\theta^{*}}$, we have $$ Var_\theta \hat{\theta} \leq Var_\theta \hat{\theta}^{*}, \text { for } \theta \in \Theta $$ ( $\hat{\theta}$ is uniformly better than $\hat{\theta^{*}}$ in variance.)

Lemma 2 (Uniqueness of UMVUE). Assume that a UMVUE for $\theta$ exists and is denoted by $T^{*}$. Then $T^{*}$ is unique.

对于如何寻找这个统计量,我们将会引入两种方式,这里先介绍第一种:Cramer-Rao lower bound for variance of unbiased estimator :

Some regularity conditions:
(a) Parameter space $\Theta$ is an open interval. $(a, \infty),(a, b),(b, \infty)$, a,b are constants not depending on $\theta$.
(b) Set ${x: f(x, \theta)=0}$ is independent of $\theta$.
(c) $\int \frac{\partial f(x, \theta)}{\partial \theta} d x=\frac{\partial}{\partial \theta} \int f(x, \theta) d x=0$
(d) If $T=t\left(x_1, \ldots, x_n\right)$ is an unbiased estimator, then $$ \int t \frac{\partial f(x, \theta)}{\partial \theta} d x=\frac{\partial}{\partial \theta} \int t f(x, \theta) d x $$

Thm 3(Cramer-Rao (C-R)). Suppose that the regularity conditions hold. If $\hat{\tau}(\theta)=t\left(X_1, \ldots, X_n\right)$ is unbiased for $\tau(\theta)$, then

$$ \operatorname{Var}_\theta \hat{\tau}(\theta) \geq \frac{\left(\tau^{\prime}(\theta)\right)^2}{n E_\theta\left[\left(\frac{\partial \ln f(x, \theta)}{\partial \theta}\right)^2\right]}=\frac{\left(\tau^{\prime}(\theta)\right)^2}{-n E_\theta\left[\left(\frac{\partial^2 \ln f(x, \theta)}{\partial \theta^2}\right)\right]} \text { for } \theta \in \Theta $$

Remark. C-R bound同样可以表示为 $\operatorname{Var}_\theta \hat{\tau}(\theta) \geq \frac{\left(\tau^{\prime}(\theta)\right)^2}{nI(\theta)}$,其中 $I(\theta)$ 为 Fisher information。若我们可以找到一个无偏统计量,使得该无偏统计量的方差为 C-R bound,则该统计量就是 UMVUE,并且我们定义这个统计量为 efficient estimator(有效估计量),是所有无偏估计量中最有效的统计量,并且定义,C-R bound 比上真实估计量的方差为 efficiency。

Sufficient statistics and Coditional expectation

我们在上一节中介绍了 UMVUE,但是使用 C-R bound来寻找UMVUE 不是一种很聪明的方法, 因为我们只是去验证而不是去寻找,本章节我们将引入充分统计量(Sufficient statistics) 的概念,并发掘一套寻找UMVUE的方法。

Def 1(Sufficient statistics). Let $X_1;\cdots ;X_n$ be a random sample from a distribution with p.d.f $f(x; \theta); \theta \in \Theta$. We call a statistic $U = u(X_1;\cdots ;X_n)$ a Sufficient statistics if, for any value $U = u$, the conditional p.d.f $f(X_1;\cdots ;X_n;\theta\mid u)$ and its domain all not depend on parameter $\theta$.
Remark. 可以认为充分统计量蕴含的随机样本中关于原始分布的所有信息。

实际上,随机样本全体以及其顺序统计量全体都是充分统计量,但是我们却很少会直接用它们当作统计量, 因为这样的统计量维度太高,难以应用。另外,充分统计量的函数也是充分统计量(该函数不可引入 $\theta$)。

Def 2(minimal sufficient statistics). A statistic $T$ is minimal sufficient if $T$ is sufficient, and for every sufficient statistic $\tilde{T}$ there exists a function $f$ such that $T=f(\tilde{T})$ (a.e. $\mathcal{P}$ ). Here (a.e. $\mathcal{P}$ ) means that the set where equality fails is a null set for every $P \in \mathcal{P}$.

判断充分统计量并不是一种很简单的事情,有些条件分布的 p.d.f. 难以计算,下面的分解定理将会 极大的简化充分统计量的判断:

Thm 3(Factorization Theorem). Let $X_1, \ldots, X_n$ be a random sample from a distribution with p.d.f $f(x, \theta)$. A statistic $U=u\left(X_1, \ldots, X_n\right)$ is sufficient for $\theta$ $\iff$ there exists functions $K_1, K_2 \geq 0$ such that the joint p.d.f of $X_1, \ldots, X_n$ may be formulated as $$f\left(x_1, \ldots, x_n, \theta\right)=K_1\left(u\left(X_1, \ldots, X_n\right), \theta\right) K_2\left(x_1, \ldots, x_n\right)$$ where $K_2$ is not a function of $\theta$.
Remark. 使用该方法判断统计量的充分性,甚至连条件概率都不用计算,只需计算联合分布是否可以拆分 即可,该定理同样适应高维的情况,可以判断高维充分统计量;另外,充分统计量可能以多种形式存在,不唯一。

Def 4(Conditional expectation,variance). Conditional expectation of $Y$ given $X=x$ is $$ \mathrm{E}(Y \mid x)=\int_{-\infty}^{\infty} y f(y \mid x) d y $$ The random conditional expectation $\mathrm{E}(Y \mid X)$ is function $\mathrm{E}(Y \mid x)$ with x replaced by X. Conditional variance of Y given $X=x$ is $$ \operatorname{Var}(Y \mid x)=\mathrm{E}\left[(Y-\mathrm{E}(Y \mid x))^2 \mid x\right]=\mathrm{E}\left(Y^2 \mid x\right)-(\mathrm{E}(Y \mid x))^2 $$ The conditional variance $\operatorname{Var}(Y \mid X)$ is $\operatorname{Var}(Y \mid x)$ replacing x by X .
Remark. 给定 $X = x$ 的条件下,条件期望和方差都是 $x$ 的函数,若给定条件换为随机变量, 条件期望和方差为随机变量,并且可以看作条件的函数;在严格的测度论定义中,取代given条件的位置是随机 变量生成的 $\sigma-$代数,可以理解为一种信息流。

Thm 5. Let $Y$ and $X$ be two r.v.’s.
(a) $E[E(Y \mid x)]=E(Y)$
(b) $\operatorname{Var}(Y)=E(\operatorname{Var}(Y \mid X))+\operatorname{Var}(E(Y \mid X))$
Reamrk.(1)可以看作tower性质的一个特例; (2)可以得出一个不等式 $\operatorname{Var}(Y)\geq\operatorname{Var}(E(Y \mid X))$, 可以直观的理解为,依据已经发生事件做出的判断相比于直接对事件本身做出判断的方差更小,更可靠。

Completeness

Thm 1(Rao-Blackwell). If $\hat{\tau}\left(X_1, \ldots, X_n\right)$ is unbiased for $\tau(\theta)$ and $U$ is a sufficient statistic, then
(a) $E_\theta\left[\hat{\tau}\left(X_1, \ldots, X_n\right) \mid U\right]$ is a statistic.
(b) $E_\theta\left[\hat{\tau}\left(X_1, \ldots, X_n\right) \mid U\right]$ is unbiased for $\tau(\theta)$.
(c) $Var_\theta\left(E\left[\hat{\tau}\left(X_1, \ldots, X_n\right) \mid U\right]\right) \leq Var_\theta\left(\hat{\tau}\left(X_1, \ldots, X_n\right)\right), \forall \theta \in \Theta$.
Remark. 若 $U$ 不是充分统计量,则 $E_\theta\left[\hat{\tau}\left(X_1, \ldots, X_n\right) \mid U\right]$ 未必是统计量,因为可能与 $\theta$ 相关。

If $\hat{\tau}(\theta)$ is an unbiased estimator for $\tau(\theta)$ and $U_1, U_2, \ldots$ are sufficient statistics, then we can improve $\hat{\tau}(\theta)$ with the following fact:

$$ \begin{aligned} \operatorname{Var}_\theta\left(\mathrm{E}\left[\hat{\tau}(\theta) \mid U_1\right]\right) & \leq \operatorname{Var}_\theta \hat{\tau}(\theta) \\ \operatorname{Var}_\theta \mathrm{E}\left(\mathrm{E}\left(\hat{\tau}(\theta) \mid U_1\right) \mid U_2\right) & \leq \operatorname{Var}_\theta \mathrm{E}\left(\hat{\tau}(\theta) \mid U_1\right) \\ \operatorname{Var}_\theta \mathrm{E}\left[\mathrm{E}\left(\mathrm{E}\left(\hat{\tau}(\theta) \mid U_1\right) \mid U_2\right) \mid U_3\right] & \leq \operatorname{Var}_\theta \mathrm{E}\left(\mathrm{E}\left(\hat{\tau}(\theta) \mid U_1\right) \mid U_2\right) \end{aligned} $$

Will this process ends with Cramer-Rao lower bound ? This can be solved with “complete statistic”.

事实上,由Rao-Blackwell保证,无偏估计量只需given一个充分统计量,就可以保证防擦好不大于任意无偏估计量的方差,即上述过程中条件期望的方差都是相等的,在后文中我们还将说明,若given的充分统计量满足完备性的条件,则上述过程中的条件期望在概率1的意义下都是同一个r.v.。

为了引入完备统计量的概念,我们先引入完备分布族:

Def 2(完备分布族). 设 $X$ 所在的分布族为 $\mathcal{F}=\{p(x ; \theta): \theta \in \Theta\}$ ,如果对任何实值函数 $\varphi(x)$ ,由期望 $$ \mathrm{E}_\theta \varphi(X)=0, \quad \forall \theta \in \Theta $$ 为 0 可以推出随机变量为 0 ,即

$$ \mathrm{P}_\theta\{\varphi(X)=0\}=1, \quad \forall \theta \in \Theta $$ 那么这个分布族称为完备分布族.
Remark. 上述定义可以由积分诱导的内积来表述,即 $<\phi,f_X> = 0\Rightarrow \phi = 0(a.s.),\forall \theta$,任何一个与分布 $f_X$ 正交的函数必定等于0(a.s.),可见分布 $f_X$ 已经张成的整个函数空间,这也是完备性的来源,是由正交函数论中的完备演化来的。

例1 二项分布族 $\{B(n, p): 0<p<1\}$(n已知)是完备的.
解: 依照定义,假设 $\varphi(X)$ 期望为 0 ,则

$$ \mathrm{E}_p \varphi(X)=\sum_{x=0}^n \varphi(x)\binom{n}{x} p^x(1-p)^{n-x}=0, \quad \forall p \in(0,1) $$

作相同变换,约去非零常数 $$ \sum_{x=0}^n\binom{n}{x} \varphi(x) \theta^x=0, \quad \forall \theta>0 $$ 这是关于 $\theta$ 的多项式,而 $\binom{n}{x}$ 不为 0 ,也就只能 $\varphi(x)$ 为 0 了,所以依照定义就有 n 已知的二项分布族是完备的.

例2 正态分布族 $\left\{N\left(0, \sigma^2\right): \sigma>0\right\}$ 不是完备的.
解: 只需找一个反例,取 $\varphi(X)=X$ ,它的期望显然为 0 ,但是随机变量取 0 的概率为 $$ \mathrm{P}_\theta{\varphi(X)=0}=0 \neq 1 $$ 因此该分布族不是完备的.

这两个例子中的分布族都是由一个随机变量诱导出的,自然而然的,我们想以统计量诱导出的分布族的完备性来定义统计量的完备性:如果一个统计量,它诱导出的分布族是完备的,那么它是完备统计量

对比讲义上的解释:
Def 3(complete statistic). $X_1, \ldots, X_n$ is random sample from $f(x, \theta)$. $A$ statistic $U=u\left(X_1, \ldots, X_n\right)$ is a complete statistic if for any function $h(U)$ such that $$ E_\theta(h(U))=0, \forall \theta \in\Theta \text{ ,then } P_\theta(h(U)=0)=1, \text{for } \theta \in \Theta $$

两种定义方式是等价的,但是由诱导出的分布族来定义更方便我们理解完备性的由来

那么完备性有没有什么对于统计的特殊意义呢?事实上,如果充分统计量表示统计量包含了样本中关于 $\theta$的所有信息,完备统计量就是不包含无用于推断 $\theta$ 的随机信息,有下面两个定理来帮助理解:

Thm 4. If $T$ is complete and sufficient, then $T$ is minimal sufficient.
Remark. 将该定理反过来说是不对的,即最小充分统计量不一定是完备的。

Def 5(Ancillary statistic). 设 $P_\theta$ 是一概率模型,其中 $\theta$ 是参数。若对于来自样本的数据 $\mathbf{Y}$ ,统计量 $T(\mathbf{Y})$ 的分布不依赖于 $\theta$ ,则称 $T(\mathbf{Y})$ 是关于 $\theta$ 的辅助统计量。这即是说,对于任何博雷尔集 $A$ ,有 $P_\theta(T(\mathbf{Y}) \in A)=\mu(A)$ ,其中 $\mu(\cdot)$ 是不依赖于 $\theta$ 的概率测度。

Theorem 5(Basu). If $T$ is complete and sufficient for $\mathcal{P}=\left\{P_\theta: \theta \in\Omega\right\}$, and if $V$ is ancillary, then $T$ and $V$ are independent under $P_\theta$ for any $\theta \in \Omega$.
Remark. Basu 定理表示,完备充分统计量与辅助统计量独立,说明该统计量既包含了样本中所有关于 $\theta$ 的信息,又不包含 $\theta$ 以外的信息,该信息量是相当 “紧致” 的,这也与定理4中完备充分统计量是最小充分统计量的说法相互映照。

完备充分统计量具有很好的性质,充分性可以由分解定理快速的验证,但是完备性却难以验证,但是若是分布满足指数族分布,就有下面的充分性定理帮助验证:

Thm 6. Let $X_1, \ldots, X_n$ be a random sample from $f(x, \theta)$ which belongs to an exponential family as $$ f(x, \theta)=e^{a(x) b(\theta)+c(\theta)+d(x)}, l<x<q $$ Then $\sum^n a\left(X_i\right)$ is a complete statistic.
Remark. 对上述等式再运用分解定理,不难发现, $\sum^n a\left(X_i\right)$ 其实是完备充分统计量.

铺垫了那么多之后,我们引出寻找UMVUE的主定理:

Thm 7(Lehmann-Scheffe). Let $X_1, \ldots, X_n$ be a random sample from $f(x, \theta)$. Suppose that $U=u\left(X_1, \ldots, X_n\right)$ is a complete and sufficient statistic. If $\hat{\tau}=t(U)$ is unbiased for $\tau(\theta)$, then $\hat{\tau}$ is the unique function of $U$ unbiased for $\tau(\theta)$ and is a UMVUE of $\tau(\theta)$. (Unbiased function of complete and sufficient statistic is UMVUE.)
Remark 1. 这里的 unique 指的是若有两个 UMVUE $X$ 和 $Y$,则 $P(X=Y) = 1$.
Remark 2. 在该定理的证明中,完备性保证了唯一,充分性保证了统计量的存在,Rao-Blackwell 定理保证了方差最小,该定理将前面的知识都联系在了一起,我们所定义的无偏性、充分性、完备性都是有用的,值得注意的是,UMVUE无法保证一致性。


Method: Two ways in constructing UMVUE based on a complete and sufficient statistic $U$ :
(a) If $T$ is unbiased for $\tau(\theta)$, then $\mathrm{E}(T \mid U)$ is the UMVUE of $\tau(\theta)$. This is easy to define but difficult to transform it in a simple form.
(b) If there is a constant such that $\mathrm{E}(U)=f \circ \theta$, then $T=f^{-1}\circ U$ is the UMVUE of $\theta$.

Example :
(a)Let $X_1, \ldots, X_n$ be a random sample from $U(0, \theta)$.

Want UMVUE of $\theta$.
sol: $Y_n=\max \left\{X_1, \ldots, X_n\right\}$ is a complete and sufficient statistic .
The p.d.f of $Y_n$ is $$ \begin{aligned} f_{Y_n}(y, \theta) & =n\left(\frac{y}{\theta}\right)^{n-1} \frac{1}{\theta}=n \frac{y^{n-1}}{\theta^n}, 0<y<\theta \\ \mathrm{E}\left(Y_n\right) & =\int_0^\theta y n \frac{y^{n-1}}{\theta^n} d y=\frac{n}{n+1} \theta . \end{aligned} $$ We then have $\mathrm{E}\left(\frac{n+1}{n} Y_n\right)=\frac{n+1}{n} \mathrm{E}\left(Y_n\right)=\theta$.
So, $\frac{n+1}{n} Y_n$ is the UMVUE of $\theta$.

(b) $X_1, \ldots, X_n \stackrel{i i d}{\sim} \operatorname{Possion}(\lambda)$. Want UMVUE of $e^{-\lambda}$.
sol: The p.d.f of X is $$ f(x, \lambda)=\frac{1}{x!} \lambda^x e^{-\lambda}=e^{x \ln \lambda-\lambda-\ln x!} $$ $\Rightarrow \sum_{i=1}^n X_i$ is complete and sufficient.
$\mathrm{E}\left(I\left(X_1=0\right)\right)=P\left(X_1=0\right)=f(0, \lambda)=e^{-\lambda}$ where $I\left(X_1=0\right)$ is an indicator function.
$\Rightarrow I\left(X_1=0\right)$ is unbiased for $e^{-\lambda}$
$\Rightarrow \mathrm{E}\left(I\left(X_1=0\right) \mid \sum_{i=1}^n X_i\right)$ is UMVUE of $e^{-\lambda}$.


上述章节为统计学(一)的所有内容,下面的章节为统计学(二)的内容。

Confidence Interval

在谈置信区间之间,我们先补充一些必要的概率论知识。

Order statistic, t-distribution

首先再谈变数转换,对于 $(X_1,\cdots,X_n)\rightarrow (Y_1,\cdots,Y_n)$ 若存在双射关系,则可以直接通过求Jacobian行列式的方式进行 p.d.f. 的转换,但若是序列之间并不存在双射关系,该如何做变数转换呢,下面给出一种通用的算法:

设 $(X_1,\cdots,X_n)\in \mathscr{A},$ 映射到的 $(Y_1,\cdots,Y_n)\in \mathscr{B}$ ,we can represent $\mathscr{A}$ as the union of a finite number, say $k$, of mutually disjoint sets $A_1, A_2, \ldots, A_k$ so that $$ y_1=u_1\left(x_1, x_2, \ldots, x_n\right), \ldots, \quad y_n=u_n\left(x_1, x_2, \ldots, x_n\right) $$ define a one-to-one transformation of each $A_i$ onto $\mathscr{B}$. Thus, to each point in there will correspond exactly one point in each of $A_1, A_2, \ldots, A_k$. Let $$ \begin{aligned} & x_1=w_{1 i}\left(y_1, y_2, \ldots, y_n\right), \\ & x_2=w_{2 i}\left(y_1, y_2, \ldots, y_n\right), \\ & \vdots \\ & x_n=w_{n i}\left(y_1, y_2, \ldots, y_n\right), \end{aligned} $$ denote the $k$ groups of $n$ inverse functions, one group for each of these $k$ transformations. Let the first partial derivatives be continuous and let each $$ J_i=\left|\begin{array}{cccc} \frac{\partial w_{1 i}}{\partial y_1} & \frac{\partial w_{1 i}}{\partial y_2} & \cdots & \frac{\partial w_{1 i}}{\partial y_n} \\ \frac{\partial w_{2 i}}{\partial y_1} & \frac{\partial w_{2 i}}{\partial y_2} & \cdots & \frac{\partial w_{2 i}}{\partial y_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial w_{n i}}{\partial y_1} & \frac{\partial w_{n i}}{\partial y_2} & \cdots & \frac{\partial w_{n i}}{\partial y_n} \end{array}\right|, i=1,2, \ldots, k $$ be not identically equal to zero in $\mathscr{B}$. From a consideration of the probability of the union of $k$ mutually exclusive events and by applying the change of variable technique to the probability of each of these events, it can be seen that the joint p.d.f. of $Y_1=u_1\left(X_1, X_2, \ldots, X_n\right)$, $$ \begin{aligned} & Y_2=u_2\left(X_1, X_2, \ldots, X_n\right), \ldots, Y_n=u_n\left(X_1, X_2, \ldots, X_n\right), \text { is given by } \\ & g\left(y_1, y_2, \ldots, y_n\right)=\sum_{i=1}^k\left|J_i\right| h\left[w_1\left(y_1, \ldots, y_n\right), \ldots, w_{n i}\left(y_1, \ldots, y_n\right)\right] \end{aligned} $$ provided that $\left(y_1, y_2, \ldots, y_n\right) \in \mathscr{B}$, and equals zero elsewhere.
Remark. 即对 $\mathscr{B}$ 来说,每一个元素都会同时对应 $\mathscr{A}$ 中的 $k$ 个元素,将这 $k$ 个元素所在的区域划分,算出来的概率再进行累加,就是我们需要获得的 p.d.f.

引入该变数转换的原因是,顺序统计量与原来的随机样本就不是 $1-1$ 的函数转换,因为随机样本的大小我们是未知的,若有 $n$ 个样本则可以产生 $n!$ 中排列组合的逆映射方式,通过上述方法对顺序统计量的联合p.d.f.求解,有 $$ \begin{aligned} f\left(y_1, y_2, \ldots, y_n\right)&= (n!) f\left(y_1\right) f\left(y_2\right) \cdots f\left(y_n\right) , \quad a<y_1<y_2<\cdots<y_n<b, \\ & =0 \quad \text { elsewhere. } \end{aligned} $$ 重述之前的定理:
Thm 1(两个顺序统计量的pdf). Let $\left(X_1, \ldots, X_n\right)$ be a random sample from a “continuous distribution” with p.d.f. $f(x)$ and d.f $F(x)$. Then the p.d.f. of $Y_n=\max \left\{X_1, \ldots, X_n\right\}$ is $$ g_n(y)=n(F(y))^{n-1} f(y) $$ and the p.d.f. of $Y_1=\min \left\{X_1, \ldots, X_n\right\}$ is $$ g_1(y)=n(1-F(y))^{n-1} f(y) $$

有了 $p.d.f.$ 之后,我们可以通过与之前完全不同的方式求得边缘分布,直接简单粗暴的积分即可,但是对于顺序统计量的积分,积分的顺序会极大的简化运算:大的顺序统计量从大向小积分,小的顺序统计量从小向大积分,这样会使得积分的上下限变得明了,给出一个具体的例子:

EX. it is easy to express the marginal p.d.f. of any order statistic, say $Y_k$, in terms of $F(x)$ and $f(x)$. This is done by evaluating the integral $$ g_k\left(y_k\right)=\int_a^{y_k} \cdots \int_a^{y_2} \int_{y_k}^b \cdots \int_{y_{n-1}}^b n!f\left(y_1\right) f\left(y_2\right) \cdots f\left(y_n\right) d y_n \cdots $$ $$ d y_{k+1} d y_1 \cdots d y_{k-1} $$ The result is $$ \begin{aligned} g_k\left(y_k\right) & =\frac{n!}{(k-1)!(n-k)!}\left[F\left(y_k\right)\right]^{k-1}\left[1-F\left(y_k\right)\right]^{n-k} f\left(y_k\right),\quad a<y_k<b, \\ & =0 \text { elsewhere. } \end{aligned} $$

注意,顺序统计量之间未必是独立的,所以计算某两个顺序统计量的联合分布时,依然需要计算积分,这时的积分技巧与上述相同,这里直接给出答案(详细解答可看hogg的书籍): $$ \begin{aligned} f_{i j}\left(y_i, y_j\right)=\frac{n!}{(i-1)!(i-i-1)!(n-j)!} & \left[F\left(y_i\right)\right]^{i-1}\left[F\left(y_j\right)-F\left(y_i\right)\right]^{j-i-1}\left[1-F\left(y_j\right)\right]^{n-j} f\left(y_i\right) f\left(y_j\right) \end{aligned} $$ for $a<y_i<y_j<b$, and zero elsewhere.

现在来引入t分布,t分布在以后的问题中会经常遇到:

Def 2(t分布). If $Z \sim N(0,1)$ and $\chi^2(r)$ are independent, we call the distribution of the r.v. $$ T=\frac{Z}{\sqrt{\frac{\chi^2(r)}{r}}} $$ a t-distribution with $r$ degrees of freedom.
Remark. t分布是偶对称分布,有 $f_T(-x) = f_T(x)$,另外,t分布只与其自由度 $r$ 相关,由 $r$ 所唯一确定,考虑到我们之前曾经提到 $\frac{(n-1)S^2}{\sigma^2}\sim\chi^2$,可以将任意一个正态分布标准化为标准正态分布,再比上 $\sqrt{\frac{S^2}{n-1}}=\frac{S}{\sigma}$ 来得到 t分布。

Confidence Interval

之前介绍的点估计是一种经典的统计推断方法,在一致性的保证下,估计量可以依概率收敛到真实参数,但是点估计的缺陷也是明显的,对于连续分布而言,我们的估计量是服从某一分布的,那么就必然有 $P(\hat{\theta} = \theta) = 0$,一个点对于连续分布是测度为零的,或许这在提示我们,想获得概率上的保证,估计一个区间也是可行的方法:

Def 1(C.I.). Suppose that we have a random sample from $f(x, \theta)$. For $0<\alpha<1$, if there exists two statistics $T_1=t_1\left(X_1, \ldots, X_n\right)$ and $T_2=t_2\left(X_1, \ldots, X_n\right)$ satisfying $$ 1-\alpha=P\left(T_1 \leq \theta \leq T_2\right) $$ We call the random interval $\left(T_1, T_2\right)$ a $100(1-\alpha) \%$ confidence interval(C.I.) of parameter $\theta$.
Remark. $100(1-\alpha) \%$ 又叫做 confidence coefficient,这个概率的含义是一个有意思的事情,它并不是说观测出某一组随机样本,它们的 $\theta$ 在置信区间内的概率是$100(1-\alpha) %$ ,因为对于已经观测出来的样本,他存在的概率只有 $0$ 或 $1$,这个概率的含义是:如果我们可以 Re-sampling 很多样本,那么在很多样本中,$\theta$ 属于置信区间的概率(频率比总实验次数) 为 $100(1-\alpha) \%$ .

Constructing C.I. by pivotal quantity(枢轴量,关键量):

Def 2(pivotal quantity). A function of random sample and parameter, $Q=q\left(X_1, \ldots, X_n, \theta\right)$, is called a pivotal quantity if its distribution is independent of $\theta$(与 $\theta$ 无关系).

满足定义的枢轴量会有很多,哪些枢轴量对寻找置信区间有帮助呢?下面的一段话回答了这个问题:

With a pivotal quantity $q\left(X_1, \ldots, X_n, \theta\right)$, there exists $a, b$ such that $$ 1-\alpha=P\left(a \leq q\left(X_1, \ldots, X_n, \theta\right) \leq b\right), \forall \theta \in \Theta $$ The interest of pivotal quantity is that there exists statistics $T_1=t_1\left(X_1, \ldots, X_n\right)$ and $T_2=t_2\left(X_1, \ldots, X_n\right)$ with the following 1-1 transformation $$ a \leq q\left(X_1, \ldots, X_n, \theta\right) \leq b \iff T_1 \leq \theta \leq T_2 $$ Then we have $1-\alpha=P\left(T_1 \leq \theta \leq T_2\right)$ and $\left(T_1, T_2\right)$ is a $100(1-\alpha) %$ C.I. for $\theta$

Remark. 现在来思考一下枢轴量的必要性,我们需要找到一个随机区间,使得参数属于该区间的概率为一个常数,该条件需要对所有的参数空间成立,如果一开始的枢轴量是依赖于 $\theta$ 的,那么该概率常数必然与 $\theta$ 相关,就违反了我们的条件;我们感兴趣的枢轴量是可以一对一的转化成 $T_1 \leq \theta \leq T_2$的形式的,举个例子说明这一点:

EX. $X_1, \ldots, X_n \stackrel{i i d}{\sim} N\left(\mu, \sigma^2\right)$. $\mu,\sigma$ 为参数,求 $\mu$ 的置信区间。 $$ \begin{gathered} \left\{\begin{array} { l } { \overline { X } \sim N ( \mu , \frac { \sigma ^ { 2 } } { n } ) } \\ { \frac { ( n - 1 ) s ^ { 2 } } { \sigma ^ { 2 } } \sim \chi ^ { 2 } ( n - 1 ) } \end{array} \text { indep. } \Rightarrow \left\{\begin{array}{l} \frac{\overline{X}-\mu}{\sigma / \sqrt{n}} \sim N(0,1) \\ \frac{(n-1) s^2}{\sigma^2} \sim \chi^2(n-1) \end{array}\right.\right. \text { indep. } \\ T=\frac{\frac{\bar{X}-\mu}{\sigma / \sqrt{n}}}{\sqrt{\frac{(n-1) s^2}{\sigma^2(n-1)}}}=\frac{\bar{X}-\mu}{s / \sqrt{n}} \sim t(n-1) \end{gathered} $$ Let $t_{\frac{\alpha}{2}}$ satisfies

$$ \begin{aligned} 1-\alpha & =P\left(-t_{\frac{\alpha}{2}} \leq \frac{\bar{X}-\mu}{s / \sqrt{n}} \leq t_{\frac{\alpha}{2}}\right) \\ & =P\left(-t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}} \leq \bar{X}-\mu \leq t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}}\right) \\ & =P\left(\bar{X}-t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}} \leq \mu \leq \bar{X}+t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}}\right) \end{aligned} $$

$\Rightarrow\left(\bar{X}-t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}}, \bar{X}+t_{\frac{\alpha}{2}} \frac{s}{\sqrt{n}}\right)$ is a $100(1-\alpha) %$ C.I. for $\mu$.
Remark. 可见,由于 $T$ 可以将 $\mu$ 剥离出来,而保证 $T_1,T_2$ 都是统计量,这才使得 $T$ 成为了我们感兴趣的枢轴量,若我们只是将 $\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}$ 当作枢轴量,是找不到置信区间的。

EX. 同样是上述例子,找 $\sigma^2$ 的C.I.,需要利用 $\frac{(n-1)S^2}{\sigma^2}\sim \chi^2(n-1)$ 来当枢轴量。

另外,有时候还需找到样本均值的差的C.I.,即 $\mu_1-\mu_2$ 的 C.I.,还是那一套方法,构造枢轴量,然后化简。


这里总结一下常见的枢轴量选取:

  1. 推断 $\mu$, 方差已知,构造标准正态分布作为枢轴量;方差未知,先构造标准正态分布,再构造t分布消去 $\sigma^2$作为枢轴量。
  2. 推断 $\sigma^2$,构造 $\chi^2$分布作为枢轴量。

近似置信区间

在有些时候我们并不知道母体分布,需要使用大样本理论去推断参数的近似置信区间。

下面介绍一些简单的大样本理论,首先说明一点,之前我们提到了mgf唯一对应随机变量的分布,其实这种对应是分布函数意义上的确定,也就是在分布函数的意义上相等,如果用mgf进行极限分布的推断,只可以得到依分布收敛的结果。

定理 1 (连续性定理). 令随机变量 $X_n$ 和 $X$ 的 CDF $F_n$ 和 $F$ 为对应它们的矩母函数为 $\phi_n$ 和 $\phi$ .那么, $\phi_n(t) \rightarrow \phi(t) \forall t\iff X_n \xrightarrow{\mathcal{D}} X$ .

定理 2(经典CLT). 令 $\left\{X_n\right\}$ 为均值为 $\mu$ ,方差为 $\sigma^2<\infty$ 的i.i.d随机变量序列.那么, $$ \frac{\overline{ X } -\mu}{\sigma \sqrt{n}} \xrightarrow{\mathcal{D}} Z $$ 这里 $Z \sim N(0,1)$ .

定理3 (连续映射定理). 令 $g: \mathbb{R}^{\mathrm{p}} \mapsto \mathbb{R}^{\mathrm{m}}$ 为在集合 C 中几乎处处连续的映射。如果 $\mathrm{X}_{\mathrm{n}}$ 依概率/几乎处处收玫 $/$ 依分布收敛到 X ,则 $g\left(\mathrm{X}_{\mathrm{n}}\right)$ 依概率/几乎处处收玫/依分布收敛到 $g(\mathrm{X})$

定理 4(Slutsky’s theorem). 设 $\left(X_n\right),\left(Y_n\right)$ 为随机标量,向量或矩阵序列。若 $\left(X_n\right)$ 依分布收敛至随机元素 $X$ ,且 $\left(Y_n\right)$ 依概率收敛至常数 $c$ ,则

  • $\left(X_n+Y_n\right) \xrightarrow{d} X+c ;$
  • $\left(X_n Y_n\right) \xrightarrow{d} X c$ ;
  • $\left(X_n / Y_n\right) \xrightarrow{d} X / c$, 若 $c$ 可逆。

现在讨论近似区间估计的问题,这类问题实际上很简单,即我们对样本的分布一无所知,但可以由CLT保证当样本量趋于无穷时,$\frac{\overline{ X } -\mu}{\sigma \sqrt{n}} \xrightarrow{\mathcal{D}} Z$,我们唯一需要做的就是构造 $\frac{\overline{ X } -\mu}{\sigma \sqrt{n}}$,但是这种枢轴量是无意义的,我们可以进一步构造 $\frac{\sigma}{S}\frac{\overline{ X } -\mu}{\sigma \sqrt{n}}$,再由Slutsky定理($S^2$依概率收敛到 $\sigma^2$,另外一半式子依分布收敛到标准正态分布)、连续映射定理($S$依概率收敛到 $\sigma$)以及CLT得到 $\frac{\overline{ X } -\mu}{S \sqrt{n}}$依分布收敛到标准正态分布,则可以反解出 $\mu$ 的近似置信区间。

一般而言,对于 $\mu$ 的近似区间估计,最后都会将分母的 $\sigma$ 换成 $S$,再由Slutsky定理保证收敛。


还有一些区间估计的内容,我们在假设检验的章节进行描述。

假设检验

我们所说的统计模型可以分为估计和假设检验两类,假设检验是不可以用估计的方法去做的。

在统计推断中, 参数的估计量的导 出是基于一些原则(例如无偏性、不变性、充分性、替代原理、似然原理、Bayes 原理等), 而不是基于损失函数和风险函数.

基本模型

Hogg 的书中这一部分描写的有点让人困扰,接下来使用邵军老师数理统计的内容,这部分内容的写法更加偏向于描述而不是下定义:

设 $X$ 是取自总体 $P \in \mathcal{P}$ 的一个样本.统计决策是指我们观察 $X$ 之后采取的行动,例如,关于 $P$ 的结论或 $P$ 的特征.贯穿本节,我们用 $\mathbb{A}$ 表示可允许行动的集合。设 $\mathcal{F}_{\mathbf{A}}$ 是 $\mathbb{A}$ 上的 $\sigma$ 域。可测空间 $\left(\mathbb{A}, \mathcal{F}_{\mathbf{A}}\right)$ 称为行动空间。设 $\mathscr{X}$ 是 $X$ 的值域, $\mathcal{F}_{\mathscr{X}}$ 是 $\mathscr{X}$ 的一个 $\sigma$ 域。决策准则是 $\left(\mathscr{X}, \mathcal{F}_{\mathscr{X}}\right)$ 到 $\left(\mathbb{A}, \mathcal{F}_{\mathrm{A}}\right)$ 上的一个可测函数(统计量)$T$ 。如果决策准则 $T$ 已选定,则一旦观测到 $X$ 我们就可采取行动 $T(X) \in \mathbb{A}$ .

设 $\mathcal{P}$ 是一分布族, $\mathcal{P}_0 \subset \mathcal{P}, \mathcal{P}_1=\left\{P \in \mathcal{P}: P \notin \mathcal{P}_0\right\}$ .假设检验问题可描述为确定下面两种说法哪一个为真: $$ H_0: P \in \mathcal{P}_0 \quad \text { 对 } \quad H_1: P \in \mathcal{P}_1 \text {. } $$ 这里,$H_0$ 称为原假设(有时也称为零假设和虚无假设,一般是保守的,也是主要去论证的),而 $H_1$ 称为备择假设.这个问题的行动空间仅有两个元素,即, $\mathbb{A}=\{0,1\}$ ,其中 0 是接受 $H_0$ 的行动,而 1 是拒绝 $H_0$ 的行动.决策准则(Rule)称为检验(Test).因为检验 $T(X)$ 是 $\mathscr{X}$ 到 $\{0,1\}$ 的函数,它一定有形式 $I_C(X)$ ,其中 $C \in \mathcal{F}_{\mathscr{X}}$ 称为 $H_0$ 对 $H_1$ 的拒绝域或临界域

检验给出的假设 $H_0$ 对 $H_1$ ,我们只可能犯两类统计错误:当 $H_0$ 为真时拒绝 $H_0$(称为第 I 类错误)以及当 $H_0$ 为假时接受 $H_0$(称为第 II 类错误).在统计推断中,检验 $T$, 一 个 $\mathscr{X}$ 到 $\{0,1\}$ 的统计量,可由犯两类错误的概率来进行评价: $$ \alpha_T(P)=P(T(X)=1), \quad P \in \mathcal{P}_0 $$ 和 $$ 1-\alpha_T(P)=P(T(X)=0), \quad P \in \mathcal{P}_1 $$ 分别记作 $\alpha_T(\theta)$ 和 $1-\alpha_T(\theta)$(注意,这里的 $P \in \mathcal{P}_0$ 是为了表示错误发生的条件,$\alpha_T(P)$ 本身就指的是在给定 $P$ 下拒绝 $H_0$ 的概率,在有些书中也称效用函数,Power function) ,若 $P$ 是参数为 $\theta$ 的参数族。然而,即使对一个很简单的问题和很简单的检验类,最优决策准则(检验)不存在。也就是说,犯两种错误的概率不能同时最小化.而且,对固定的样本量,这两类犯错误概率不能同时以一固定的 $\alpha \in(0,1)$ 为界.

因此,寻找"最优"检验的通常办法是将犯其中一个错误的概率限定在一个小的界 $\alpha$ 内,如 $\alpha_T(P), P \in \mathcal{P}_0$ ,然后试图最小化犯另一个错误的概率 $1-\alpha_T(P), P \in \mathcal{P}_1$ ,在以下约束条件下, $$ \sup _{P \in \mathcal{P}_0} \alpha_T(P) \leqslant \alpha $$ 上界 $\alpha$ 称为显著性水平.左侧称为检验 $T$ 的检验水准(size)。显著性水平应为正数。

接下来我们介绍一种更加一般的观点:


一个假设的检验是一个取值为 $[0,1]$ 的统计量 $T(X)$ .当观测值 $X=x$ 时,我们以 $T(x)$ 的概率拒绝 $H_0$ ,以 $1-T(x)$ 的概率接受 $H_0$ .若 $T(X)=1$ 或 0 a.s. $\mathcal{P}$ ,则 $T(X)$ 是一个非随机化检验.否则 $T(X)$ 是一个随机化检验.对于给定的检验 $T(X)$ ,它的势函数(上述效用函数)定义为 $$ \beta_T(P)=E_P[T(X)], \quad P \in \mathcal{P} . $$ 若 $P \in \mathcal{P}_0$ ,它是 $T(X)$ 的犯第 I 类错误的概率;若 $P \in \mathcal{P}_1$ ,它是 $H_1$ 成立时,正确拒绝 $H_0$ 的概率,即检验的 power,$1-\beta_T(P)$ 为犯第 II 类错误的概率.(该定义可以涵盖上述的非随机检验的效用函数,也可以涵盖随机化检验,是更一般版本的power function)

当样本量固定时,我们不能同时使犯两类错误的概率达到最小.我们的办法是在所有 $P \in \mathcal{P}_1$ 和所有满足下面条件的 $T$ 上最大化势函数 $\beta_T(P)$(即最小化犯第 II 类错误的概率)(这种方法被称为MP检验,最优势,最大效力), $$ \sup _{P \in \mathcal{P}_0} \beta_T(P) \leqslant \alpha $$ 其中 $\alpha \in[0,1]$ 是给定的显著性水平.注意到之前我们定义了的左侧为 $T$ 的检验水准(size)。

定义 1(UMP). 检验水准为 $\alpha$ 的检验 $T_{*}$ 是一致最优势(UMP)检验,当且仅当对所有 $P \in \mathcal{P}_1$ 和显著水平为 $\alpha$ 的 $T, \beta_{T_{*}}(P) \geqslant \beta_T(P)$ .

似然比检验

我们先给出一个直观的检验方式,并之后给出其理论。

定义 1(似然比检验). 关于检验 $H_0: \theta \in \Theta_0$ 对 $H_1: \theta \in \Theta_0^{\mathrm{C}}$ 的似然比检验统计量是 $$ \lambda(\boldsymbol{x})=\frac{\sup_{\theta_0} L(\theta \mid \boldsymbol{x})}{\sup _\theta L(\theta \mid \boldsymbol{x})} $$ 任何一个拒绝区域的形式为 $\{\boldsymbol{x}: \lambda(\boldsymbol{x}) \leqslant c\}$ 的检验都叫做似然比检验(likelihood ratio test,LRT),$c$ 为任意 $[0,1]$ 的数。
注. 直观理解似然比检验,如果 $\Theta_0$ 的最大似然越大,那么比值越高,零假设为真的概率就越高,所以拒绝域设置为小于一个阈值。

最大似然比检验的求法也非常简单,先求出 $\lambda(x)$,再使得显著水平为 $\alpha$ 即可,一般对其显著水平的时候需要利用函数的一些性质得出等价判决域,该步骤的计算量很大,并且比较需要技巧性,也是后续需要引入p值的源头。


推断 $\mu$ :

  1. $\sigma$ 已知,构造标准正态
  2. $\sigma$ 未知,构造$t$ 分布,$\frac{\bar{X}-\mu}{S/\sqrt{n}}$

推断 $\sigma$:

  1. 不管 $\mu$ 已知还是未知,构造 $\Chi^2(n-1)$,$\frac{(n-1)S^2}{\sigma^2}$

假设检验的评价

定义 1(无偏性). 一个功效函数为 $\beta(\theta)$ 的检验是无偏的,如果对于每一个 $\theta^{\prime} \in \Theta_0^{\mathrm{C}}$ 和 $\theta^{\prime \prime} \in \Theta_0$ 有 $\beta\left(\theta^{\prime}\right) \geqslant \beta\left(\theta^{\prime \prime}\right)$ .

定义 2(UMP). 检验水准为 $\alpha$ 的检验 $T_{*}$ 是一致最优势(UMP)检验,当且仅当对所有 $P \in \mathcal{P}_1$ 和显著水平为 $\alpha$ 的 $T, \beta_{T_{*}}(P) \geqslant \beta_T(P)$ .

为了找出UMP,我们先从简单假设开始,对于简单假设,以下定理说明必定存在UMP,且UMP的形式是确定的。

定理 3(Neyman-Pearson引理). 对于简单假设$H_0:\theta = \theta_0,H_1:\theta = \theta_1$,构造判决域满足如下条件:
(1)$C=\left\{X\mid \frac{L(\theta_0\mid X)}{L(\theta_1\mid X)}\leq k\right\}$
(2)$\alpha = P(X\in C\mid \Theta_0)$
则该判决域为UMP。

使用N-P引理可以对简单假设的UMP进行求解,但是我们的最终目标是要求得一般假设,求取的手法一般是归化到简单假设,我们可以将常见的假设划分为以下四种:

  1. 简单假设,UMP存在。
  2. 单侧假设,$H_0:\theta = \theta_0,H_1:\theta\leq\theta_0$,UMP一般存在,可以归化到简单形式进行求解。
  3. 双侧假设,$H_0:\theta = \theta_0,H_1:\theta\neq\theta_0$,UMP一般不存在
  4. 复合单侧假设,$H_0:\theta \leq \theta_0,H_1:\theta\gneq\theta_0$,UMP可能存在,由下面的定理保证。

定义 4(单调似然比,MLR). 令 $f_n(x \mid \theta)$ 表示样本 $\boldsymbol{X}=\left(X_1, \cdots, X_n\right)$ 的联合密度(质量)函数;令 $T=T(\boldsymbol{X})$ 表示一个统计量。如果对于 $\theta_1 \in \Omega_1$ 和 $\theta_2 \in \Omega_2\left(\theta_1<\theta_2\right)$ ,都有似然比 $f_n\left(x \mid \theta_2\right) / f_n\left(x \mid \theta_1\right)$ 是关于 $T(\boldsymbol{X})$ 的单调(非增或非降的)函数,且该似然比中 $\boldsymbol{X}$ 的信息只包含在 $T(\boldsymbol{X})$ 中,那么称 $\boldsymbol{X}$ 的分布相对于统计量 $T(\boldsymbol{X})$ 具有单调似然比(Monotone Likelihood Ratio,MLR)。

定理 5(Karlin-Rubin) 考虑检验 $H_0: \theta \leqslant \theta_0$ 对 $H_1: \theta>\theta_0$ .设 $T$ 是一个关于 $\theta$ 的充分统计量并且 $T$ 的概率密度函数或概率质量函数的族 ${g(t \mid \theta): \theta \in \Theta}$ 关于 $\theta$ 具有 MLR.则对于任何 $t_0$ ,"当且仅当 $T>t_0$ 时拒绝 $H_0$"的检验是一个 UMP 水平为 $\alpha$ 的检验,其中 $\alpha=P_{\theta_0}\left(T>t_0\right)$ .
注. 这仅仅是一种形式的K-R定理,该定理对于复合单侧检验普遍成立,需要利用直觉去判断判决域,非常符合直观,个人觉得用另外一种方式概括该定理是合理的:对于复合单侧检验,若为单调似然比,则UMP存在,判决域可根据直觉判断。

F-分布

在此处暂停一下,引入F-分布,并总结一下统计学三大分布:$\Chi^2,t,F$.

引理1(关于 $\chi^2$ 随机变量的若干事实) 以 $\chi_p^2$ 记自由度为 $p$ 的 $\chi^2$ 随机变量.
a.如果 $Z$ 是 $\mathrm{N}(0,1)$ 随机变量,则 $Z^2 \sim \chi_1^2$ ,即标准正态随机变量的平方是 $\chi^2$ 随机变量;
b.如果 $X_1, \cdots, X_n$ 独立且 $X_i \sim \chi_{p_i}^2$ ,则 $X_1+\cdots+X_n \sim \chi_{p_1+\cdots+p_n}^2$ ,即独立的 $\chi^2$随机变量之和仍为 $\chi^2$ 随机变量,且其自由度为原随机变量自由度之和.

定义 2(F分布). 设随机样本 $X_1, \cdots, X_n$ 取自服从 $\mathrm{n}\left(\mu_X, \sigma_X^2\right)$ 分布的总体,随机样本 $Y_1, \cdots, Y_m$ 取自服从 $\mathrm{n}\left(\mu_Y, \sigma_Y^2\right)$ 分布的总体,且与 $X_1, \cdots, X_n$ 独立.则称随机变量 $$F=\frac{\left(S_X^2 / \sigma_X^2\right)}{\left(S_Y^2 / \sigma_Y^2\right)}$$

服从自由度为 $n-1$ 和 $m-1$ 的 Snedecor $F$ 分布 (Snedecor's $F$ distribution with $n-1$ and $m-1$ degrees of freedom)。

定理 3(F分布的性质). (1)$F\sim f(p,q)\Rightarrow \frac{1}{F}\sim f(q,p)$
(2)对于F分布的分位数,有 $F_\alpha(p,q) = \frac{1}{f_{1-\alpha}(q,p)}$

引入F分布的理由是因为需要对 $\sigma_x/\sigma_y$ 进行假设检验或者区间估计,使用F分布构造分位数会将问题变的十分简便,在之前的估计中,若比较两个分布,我们经常需要做 $\sigma_x = \sigma_y$ 的假设,现在也可以对该假设进行检验了。

p值

再谈区间估计

在习题中不难发现,假设检验和区间估计具有类似的方法,再加上我们之前在区间估计中定义置信度为 $1-\alpha$,不禁让人联想,二者是否是同一个东西。

定理 1(假设检验等价于区间估计). 对每一个 $\theta_0 \in \Theta$ ,设 $A\left(\theta_0\right)$ 是 $H_0: \theta=\theta_0$ 的一个水平为 $\alpha$ 的检验的接受区域.对每一个 $x \in \mathcal{X}$ ,在参数空间里定义一个集合 $C(x)$ $$ C(\boldsymbol{x})=\left\{\theta_0: \boldsymbol{x} \in A\left(\theta_0\right)\right\} $$ 则随机集合 $C(\boldsymbol{x})$ 是一个 $1-\alpha$ 置信集合.反之,设 $C(\boldsymbol{x})$ 是一个 $1-\alpha$ 置信集合.对任意的 $\theta_0 \in \Theta$ ,定义 $$ A\left(\theta_0\right)=\left\{\boldsymbol{x}: \theta_0 \in C(\boldsymbol{x})\right\} $$ 则 $A\left(\theta_0\right)$ 是 $H_0: \theta=\theta_0$ 的一个水平为 $\alpha$ 的检验的接受区域.

所以使用假设检验得到判决域再进行反转得到置信区间也是一种置信区间的求法;在假设检验中,我们定义了UMP,那么UMP反转到区间估计是什么呢,我们自然也想定义一个一致最优区间估计。

定义 2(一致最精确置信区间,UMA). 一个在一类 $1-\alpha$ 置信区间上最小化假值覆盖率的置信区间叫做一致最精确置信区间(Uniformly Most Accurate,UMA)。

UMP一般是单侧的,UMA也是如此,并且二者相互对应:

定理 3(UMP与UMA). 设 $\boldsymbol{X} \sim f(\boldsymbol{x} \mid \theta)$ ,其中 $\theta$ 是一个实值参数.对于每个 $\theta_0 \in \Theta$ ,设 $A^{*}\left(\theta_0\right)$ 是关于 $H_0: \theta=\theta_0$ 对 $H_1: \theta>\theta_0$ 的一个 UMP 水平 $\alpha$ 检验的接受区域。设 $C^*(\boldsymbol{x})$ 是通过反转上述 UMP 接受区域所建立的 $1-\alpha$ 置信集合.则对于任何其他的 $1-\alpha$ 置信集合 $C$ ,有 $P_\theta\left(\theta^{\prime} \in C^{*}(\boldsymbol{X})\right) \leqslant P_\theta\left(\theta^{\prime} \in C(\boldsymbol{X})\right)$ 对于所有的 $\theta^{\prime}<\theta$ 成立。

同样的,我们可以对置信区间定义无偏性:

定义 4(无偏性). 置信区间 $C(x)$ 是无偏的,如果 $P_\theta(\theta^{’}\in C(X))\leq1-\alpha$ 对于所有的 $\theta\neq\theta^{’}$ 成立。
注. 假值覆盖率小于真值覆盖率,这与假设检验的思想一致,实际上二者也是相互对应的。

接下来谈一个我们在区间估计中没有提到的问题,我们自然希望一个区间估计的长度最小,那么该长度最小和UMA有何关系呢,由于下面的定理,UMA(假值覆盖率最小的置信区间)也称为Neyman最短区间,通过该名字也可以看出,极小化假值覆盖率带有期望长度最优某种程度上的保证。

定理5(Pratt定理). 设 $X$ 是一个实值随机变量,$X \sim f(x \mid \theta)$ ,其中 $\theta$ 是一个实值参数.设 $C(x)=[L(x), U(x)]$ 是一个关于 $\theta$ 的置信区间.如果 $L(x)$ 和 $U(x)$都是 $x$ 的增函数,则对于任意的值 $\theta^{*}$ ,有 $$ E_{\theta^{*}}(\text { Length }[C(X)])=\int_{\theta \neq \theta^{*}} P_{\theta^{*}}(\theta \in C(X)) \mathrm{d} \theta $$ 该定理讲的是 $C(X)$ 的期望长度等于假值覆盖概率的求和(积分),积分域取遍参数的所有假值.

渐进分析

点估计

主要分析MLE的一些渐进特性。 定理 1(MLE 的一致性). 设 $X_1, X_2, \cdots$ 是 iid $f(x \mid \theta)$ 的,$L(\theta \mid \boldsymbol{x})$ $=\prod_{i=1}^n f\left(x_i \mid \theta\right)$ 是似然函数,而 $\hat{\theta}$ 表示 $\theta$ 的 MLE.设 $\tau(\theta)$ 是 $\theta$ 的一个连续函数.那么在杂录 10.6.2 的关于 $f(x \mid \theta)$ ,从而也就是对 $L(\theta \mid \boldsymbol{x})$ 的正则性条件之下,对于每个 $\epsilon>0$ 和每个 $\theta \in \Theta$ ,有 $$ \lim _{n \rightarrow \infty} P_\theta(|\tau(\hat{\theta})-\tau(\theta)| \geqslant \epsilon)=0 $$ 这就是说,此 $\tau(\hat{\theta})$ 是 $\tau(\theta)$ 的一个一致估计量。

定义 2(极限方差). 对于一个估计量 $T_n$ ,如果 $\lim _{n \rightarrow \infty} k_n \operatorname{Var} T_n=\tau^2<\infty$ ,其中 $\left\{k_n\right\}$ 是个常数序列,则 $\tau^2$ 叫做极限方差(limiting variance)或方差的极限。

极限方差并不能很好的度量渐进性能,更常用的是下面的渐进方差:

定义 3(渐进方差). 对于一个估计量 $T_n$ ,假定有依分布收敛 $k_n\left(T_n-\tau(\theta)\right) \rightarrow \mathrm{N}(0,\sigma^2)$ ,则参数 $\sigma^2$ 叫做 $T_n$ 的渐近方差或 $T_n$ 的极限分布的方差.

渐进方差一般小于极限方差,我们希望一个统计量在极限状态下表现是UMVUE,有以下定义:

定义 4(渐进有效). 一个估计量序列 $W_n$ 关于一个参数 $\tau(\theta)$ 是渐近有效的,如果 $\sqrt{n}\left[W_n-\tau(\theta)\right] \xrightarrow{L} \mathrm{n}(0, v(\theta))$ 而且 $$ v(\theta)=\frac{\left[\tau^{\prime}(\theta)\right]^2}{\mathrm{E}_\theta\left(\left(\frac{\partial}{\partial \theta} \log f(X \mid \theta)\right)^2\right)}, $$ 就是说,$W_n$ 的渐近方差达到了 Cramér-Rao 下界.

定理 5(MLE 的渐近有效性) 设 $X_1, X_2, \cdots$ 是 iid $f(x \mid \theta), \theta$ 的 MLE 记作 $\hat{\theta}$ ,设 $\tau(\theta)$ 是 $\theta$ 的一个连续函数.那么在对 $L(\theta \mid \boldsymbol{x})$ 的正则性条件之下, $$ \sqrt{n}[\tau(\hat{\theta})-\tau(\theta)] \stackrel{L}{\rightarrow} \mathrm{n}(0, v(\theta)), $$ 其中 $u(\theta)$ 是 Cramér-Rao下界.就是说,$\tau(\hat{\theta})$ 是 $\tau(\theta)$ 的一个相合且渐近有效的估计量.

监禁有效本身就意味着一致性。

稳健性

应该是很有趣的一个主题,Casella的书并没有着重介绍,以后补上。

假设检验

主要关心LRT的渐进分布。

定理 10.3.1(LRT 的渐近分布简单 $H_0$,Wilks定理) 关于检验 $H_0: \theta=\theta_0$ 对 $H_1: \theta \neq \theta_0$ ,设 $X_1, \cdots, X_n$ 是 iid $f(x \mid \theta), \hat{\theta}$ 是 $\theta$ 的 MLE,并且 $f(x \mid \theta)$ 满足在杂录 10.6 .2 中的正则性条件.则在 $H_0$ 之下,当 $n \rightarrow \infty$ , $$ -2 \log \lambda(\boldsymbol{X}) \xrightarrow{L} \chi_1^2, $$ 其中 $\chi_1^2$ 是一个具有自由度 1 的 $\chi^2$ 分布随机变量.
注. wilks定理有高维形式,$\Chi^2$ 分布自由度$=\dim \Theta-\Theta_0$

在LRT中,决策域的具体数值非常难求,利用渐进理论可以简单的求解。

Bootstrap

自助法,一种Resampling的方法,构造插件分布(Plag-in estimator)。

自举的基本思想是,可以通过对样本数据(样本→总体)进行重新采样并对重新采样数据中的样本(→样本)执行有关样本的推理来建模。由于总体未知,因此样本统计量中相对于其总体值的真正误差是未知的。在 bootstrap-resamples 中,“总体”实际上是样本,这是已知的;因此,从重新采样的数据(重新采样→样本)中推断“真实”样本的质量是可衡量的。

使用全部数据构建估计量的计算量是非常大的,使用自助法可以减少计算量,并且达到对统计量的良好估计,一般可以具有一致性。

Resampling领域有很强的理论性和应用性,日后可以读论文去学习。

回归

关于回归的知识可以参考另一篇博客.

采样

采样并不是频率派数理统计的核心,而是贝叶斯学派的核心,频率派处理数据重点在优化,贝叶斯派重点在处理后验分布,由于后验分布非常复杂,所以会使用变分推断和采样的方式进行估计。

逆分布生成

非常朴素的生成方式,参考另一篇博客.

大部分我们熟知的分布都可以由逆分布采样。

但是后验分布的复杂程度是超出我们想象的,逆分布生成法局限性太大,需要构造其他的采样方法,接下来介绍几种Monte Carlo采样方法。

接受拒绝采样

令 $c$ 是一个常数,使得

$$ \frac{f(y)}{g(y)} \leqslant c \quad \forall y $$ 那么,我们用下述技术模拟出具有密度 $f$ 的连续分布函数的随机变量:
步骤 1: 模拟具有密度 $g$ 的随机变量 $Y$ ,并且模拟一个随机数 $U$ .
步骤 2: 如果 $U \leqslant \frac{f(Y)}{c g(Y)}$ ,那么置 $X=Y$ .否则返回步骤 1 .

命题. 由拒绝法生成的随机变量 $X$ 具有密度 $f$ .

证明. 令 $X$ 是得到的值,而以 $N$ 记必需重复的次数.那么

$$ \begin{aligned} \mathrm{P}\{X \leqslant x\} & =\mathrm{P}\left\{Y_N \leqslant x\right\} \\ & =\mathrm{P}\{Y \leqslant x \mid U \leqslant f(Y) / c g(Y)\} \\ & =\frac{\mathrm{P}\{Y \leqslant x, U \leqslant f(Y) / c g(Y)\}}{K} \\ & =\frac{\int \mathrm{P}\{Y \leqslant x, U \leqslant f(Y) / c g(Y) \mid Y=y\} g(y) \mathrm{d} y}{K} \\ & =\frac{\int_{-\infty}^x(f(y) / c g(y)) g(y) \mathrm{d} y}{K} \\ & =\frac{\int_{-\infty}^x f(y) \mathrm{d} y}{K c} \end{aligned} $$

其中 $K=\mathrm{P}\{U \leqslant f(Y) / c g(Y)\}$ ,令 $x \rightarrow \infty$ ,表明 $K=1 / c$ ,故而完成了证明,

可以参考张颢老师在现代数字信号处理(二)中贝叶斯信号处理部分,先对有限支撑集进行采样,再拓展到无限支撑集,非常直观。

但是接受拒绝采样无法处理值域无限的情况,我们先介绍Monte-Carlo积分,再引入重要性采样和MCMC。

Monte-Carlo积分

Monte-Carlo积分是一种数值积分算法,我们不关心其数值积分性质,重点关心其概率方面的性质。

考虑函数 $$ I=\int_{\Omega} f(\overline{\mathbf{x}}) d \overline{\mathbf{x}} $$ 其中 $\Omega$ 是 $\mathbb{R}^m$ 的一个子集,它的体积是: $$ V=\int_{\Omega} d \overline{\mathbf{x}} $$ 朴素的蒙特卡罗方法是在 $\Omega$ 上均匀采样点:给定 $N$ 个均匀样本, $$ \overline{\mathbf{x}}_1, \cdots, \overline{\mathbf{x}}_N \in \Omega $$ 可以被近似为: $$ I \approx Q_N \equiv V \frac{1}{N} \sum_{i=1}^N f\left(\overline{\mathbf{x}}_i\right)=V\langle f\rangle $$ 这是因为大数定律确保以下结论成立: $$ \lim _{N \rightarrow \infty} Q_N=I $$

可以对原始的函数进行拆分,得到如下形式 $$ I=\int_{\Omega} f(\overline{\mathbf{x}})p(\overline{\mathbf{x}}) d \overline{\mathbf{x}} $$ 根据WLLN,若生成 $$ \overline{\mathbf{x}}_1, \cdots, \overline{\mathbf{x}}_N \sim p(x) $$ 则可以利用 $$ I \approx \frac{1}{N} \sum_{i=1}^N f\left(\overline{\mathbf{x}}_i\right) $$ 来近似积分。

值得注意的是,Monte-Carlo积分的收敛速率和精度与梯形公式相似。

由于Monte-Carlo积分可以对后验分布进行处理,所以使用该积分进行采样的后续处理也是合适的,我们将在下一节引入重要性采样,该采样方法一开始是为了减少Monte-Carlo积分的方差,其本身也可以用来做采样。

重要性采样


这里摘取知乎上一位答主的回答:一些常见的方差缩减方法:

假设我们要估计 $\theta=\mathbb{E}_f[h(X)]$ ,其中 $X$ 的PDF是 $f$ 。令 $g$ 为另一PDF,其满足当 $f(x) \neq 0, ~ g(x) \neq 0$ ,即 $f$ 和 $g$ 的支集相等。那么, $$ \theta=\mathbb{E}_f[h(X)]=\int h(x) \frac{f(x)}{g(x)} g(x) d x=\mathbb{E}_g\left[\frac{h(X) f(X)}{g(X)}\right] $$ 初始的模拟是从密度 $f$ 中生成 $n$ 个 $X$ 的值,并取 $\hat{\theta}_n=\sum h\left(X_j\right) / n$ 。另一种模拟是从密度 $g$ 中生成 $n$ 个 $X$ 的值,并设 $$ \hat{\theta}_{n, i s}=\sum_{j=1}^n \frac{h\left(X_j\right) f\left(X_j\right)}{g\left(X_j\right)} $$ 这里 $\hat{\theta}_{n, i s}$ 叫作 $\theta$ 的重要性采样估计。我们通常定义 $h^{*}(x)=h(x) f(x) / g(x)$ ,因此 $\theta=\mathbb{E}_g\left[h^{*}(X)\right]$ 。

令 X 为联合PDF为 $f$ 的随机向量,并不失一般性地假设 $h(\mathrm{X}) \geq 0$ 。重要性采样估计的方差是 $$ \begin{gathered} \operatorname{Var}_g\left(h^{*}(\mathrm{X})\right)=\int h(\mathrm{x})^2 g(\mathrm{x}) d \mathrm{x}-\theta^2 \\ \quad=\int h(\mathrm{x})^2 \frac{f(\mathrm{x})}{g(\mathrm{x})} f(\mathrm{x}) d \mathrm{x}-\theta^2 \end{gathered} $$ 而初始估计的方差是 $\operatorname{Var}_f(h(\mathrm{X}))=\int h(\mathrm{x})^2 f(\mathrm{x}) d \mathrm{x}-\theta^2$ ,所以方差的缩减是 $$ \operatorname{Var}_f(h(\mathrm{X}))-\operatorname{Var}_g\left(h^{*}(\mathrm{X})\right)=\int h(\mathrm{x})^2\left(1-\frac{f(\mathrm{x})}{g(\mathrm{x})}\right) f(\mathrm{x}) d \mathrm{x} $$ 为了实现方差缩减,上面的积分应该是正的,这要求

  1. 当 $h^2(\mathrm{x}) f(\mathrm{x})$ 小时,$f(\mathrm{x}) / g(\mathrm{x})<1$ ,
  2. 当 $h^2(\mathrm{x}) f(\mathrm{x})$ 大时,$f(\mathrm{x}) / g(\mathrm{x})>1$ 。 若密度 $f$ 的重要的部分定义在区域 $A$ ,则上面的要求指出当 x 在 $A$ 中时,$f(\mathrm{x}) / g(\mathrm{x})$ 很小,即密度 $g$ 给 $A$ 更多的权重,所以这个方法叫作重要性采样。由于 $h$ 包含了一个稀有事件,所以在状态空间中 $h(\mathrm{x})=0$ 。我们选择 $g$ ,使得我们经常从状态空间中 $h(\mathrm{x}) \neq 0$ 的部分抽样。因此,重要性采样对于稀有事件模拟很有用。

若我们选择 $g(\mathrm{x})=h(\mathrm{x}) f(\mathrm{x}) / \theta$ ,我们会得到一个方差为 0 的估计,但这显然是不实际的。然而,这告诉我们,我们可以选择 $g(\cdot)$ ,使得它和 $h(\cdot) f(\cdot)$ 相似。特别地,我们可以使得 $g(\mathrm{x})$ 和 $h(\mathrm{x}) f(\mathrm{x})$ 在同一个值 $\mathrm{x}^{*}$ 上取到最大值,这叫作最大值原则。有时最大值原则难以应用,比如 $\mathrm{x}^{*}=\arg \max _{\mathrm{x}} h(\mathrm{x}) f(\mathrm{x})$ 有无数多解。这时我们可以通过缩放 $f$ 来选取 $g$ 。


"重要性"一词的直观解释

"重要性采样"之名正是来源于对高贡献区域的聚焦:在估计 $\mathbb{E}_p[f(X)]$ 时,只有那些使 $f(x) p(x)$ 较大的 $x$ 值对积分贡献显著,被称为"重要点"。理想的 $q(x)$ 会对这些重要点给予更高采样概率,而对贡献小的点减少采样,以降低方差。直观地,当 $f(x)$ 很大时,样本点在该处的权重 $p(x) / q(x)$ 较小,多个样本共同估计可获得稳定的贡献;而在 $f(x)$ 很小的区域,即使权重较大,其对期望的贡献也微不足道。统计理论指出,若 $q$ 对重要区域加大权重,那么估计方差会显著减小 。总之,"重要性"一词反映了:理想的采样策略应更多采样那些对期望值贡献最大的区域,从而实现高效无偏估计。

实际上可以利用Jensen不等式求出重要性采样方差的一个lower bound,在$g\propto|h|f$ 时可以取到等号,即方差最小。值得注意的是,使用重要性采样算出的方差要小于接受拒绝采样。

gpt的一份回答.

重要性采样重采样(SIR)

将重要性采样和重采样进行结合可以得到SIR:

算法步骤
设:

  • 目标分布(如后验):$\pi(x) \propto p(x) \cdot L(x)$ ,其中 $p(x)$ 是先验,$L(x)$ 是似然;
  • 提议分布(容易采样):$q(x)$ ,需满足 $\pi(x)>0 \Rightarrow q(x)>0$ ;
  • 我们的目标是得到近似服从 $\pi(x)$ 的样本集。

步骤如下:

  1. 采样阶段(Sampling)

从提议分布 $q(x)$ 中采样 $N$ 个样本: $$ x_1, x_2, \ldots, x_N \sim q(x) $$

  1. 重要性加权(Importance Weighting)

为每个样本计算未归一化权重: $$ w_i=\frac{\pi\left(x_i\right)}{q\left(x_i\right)} \propto \frac{p\left(x_i\right) L\left(x_i\right)}{q\left(x_i\right)} $$ 归一化权重: $$ \tilde{w}_i=\frac{w_i}{\sum_{j=1}^N w_j} $$

  1. 重采样阶段(Resampling)

根据归一化权重 $\tilde{w}_i$ 从 $\left\{x_1, \ldots, x_N\right\}$ 中进行有放回采样,采样 $M \leq N$ 个样本,得到: $$ x^{*(1)}, x^{*(2)}, \ldots, x^{*(M)} \sim \operatorname{Categorical}\left(\tilde{w}_1, \ldots, \tilde{w}_N\right) $$ 最终样本集 $x^{*(1)}, \ldots, x^{*(M)}$ 可视为来自目标后验 $\pi(x)$ 的近似样本。

SIR难以处理高维数据,需要进一步引入序贯蒙特卡洛采样,转化为低纬度数据后再进行SIR,我们将在后续章节介绍该技术。

方差缩减技术

对偶(立)变量法

定理 1. 设 $g$ 为单调函数,$U \sim U(0,1)$ ,则 $\operatorname{Cov}(g(U), g(1-U)) \leq 0$ 。

定理1可以推广到如下情形。

定理 2. 设 $h\left(x_1, x_2, \ldots, x_n\right)$ 是关于每个自变量单调的函数,$U_1, U_2, \ldots, U_n$ 相互独立同 $\mathrm{U}(0,1)$ 分布,则 $\operatorname{Cov}\left(h\left(U_1, U_2, \ldots, U_n\right), h\left(1-U_1, 1-U_2, \ldots, 1-U_n\right)\right) \leq 0$ 。

对均匀随机数 $U$ 最常见的变换是逆变换 $X=F^{-1}(U)$ 。下面的定理给出了提高 $I=E X$ 估计精度的方法。

定理 3(对立变量法) 设 $F(x)$ 为连续分布函数,$U \sim \mathrm{U}(0,1), X=F^{-1}(U), Y=F^{-1}(1-U)$ , $Z=\frac{X+Y}{2}$ ,则 $X$ 与 $Y$ 同分布 $F(x)$ 且 $\operatorname{Cov}(X, Y) \leq 0$ , $$ \operatorname{Var}(Z) \leq \frac{1}{2} \operatorname{Var}(X) $$

条件期望法

进行统计估计时,如果有额外的相关信息,利用这样的信息可以提高估计精度。比如,对随机变量 $Y$ ,如果 $Y$ 服从某种模型,在估计 $I=E Y$ 时应当尽量利用模型信息。

设变量 $X$ 与 $Y$ 不独立,根据Rao-Blackwell不等式: $$ \operatorname{Var}{\mathrm{E}(Y \mid X)} \leq \operatorname{Var}(Y) $$ 又 $$ \mathrm{E}{\mathrm{E}(Y \mid X)}=\mathrm{E} Y=I $$ 所以,对 $Z=\mathrm{E}(Y \mid X)$ 抽样,用 $Z$ 的样本平均值来估计 $I=\mathrm{E} Y$ 比直接用 $Y$ 的样本平均值的精度更高。这种改善随机模拟估计精度的方法叫做条件期望法,或Rao-Blackwell方法。

例14.5 设 $X \sim p(x), \varepsilon \sim \mathrm{N}\left(0, \sigma^2\right)$ 且与 $X$ 独立, $$ Y=\psi(X)+\varepsilon $$ 估计 $I=\mathrm{E} Y$ 。 可以用条件分布抽样法对二元随机向量 $\boldsymbol{Z}=(X, Y)$ 抽样产生 $Y$ 的样本。从 $p(\cdot)$ 抽样得 $X_1, X_2, \ldots, X_N$ ,独立地从 $\mathrm{N}\left(0, \sigma^2\right)$ 抽样得 $\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_N$ ,令 $$ Y_i=\psi\left(X_i\right)+\varepsilon_i, i=1,2, \ldots, N $$ 然后用 $Y_1, Y_2, \ldots, Y_N$ 的样本平均值估计 $E Y$ : $$ \hat{I}_1=\frac{1}{N} \sum_{i=1}^N Y_i $$ 估计方差为 $$ \operatorname{Var}\left(\hat{I}_1\right)=\frac{\operatorname{Var}\left(Y_1\right)}{N}=\frac{\operatorname{Var}\left(\psi\left(X_1\right)\right)}{N}+\frac{\sigma^2}{N} . $$ 另一方面,注意 $E(Y \mid X)=\psi(X)$ ,也可以只对 $X$ 抽样然后用条件期望法估计 $E Y$ : $$ \hat{I}_2=\frac{1}{N} \sum_i \psi\left(X_i\right) $$ 估计方差为 $$ \operatorname{Var}\left(\hat{I}_2\right)=\frac{\operatorname{Var}\left(\psi\left(X_1\right)\right)}{N}<\operatorname{Var}\left(\hat{I}_1\right) $$ 这个例子演示了条件期望法可以缩减误差方差的原因:对 $Y$ 抽样分成两步进行:第一步对 $X$ 抽样,第二步利用 $Y \mid X$ 的条件分布对 $Y$ 抽样。所以使用 $Y$ 的样本估计 $E Y$ 包含了第二步对 $Y$ 抽样的随机误差,而使用 $X$ 的函数 $E(Y \mid X)=\psi(X)$ 的抽样来估计 $E Y$ 则避免了第二步对 $Y$ 抽样引起的随机误差。

控制变量法

同样假设我们要用模拟来估计 $\mathrm{E}[g(\boldsymbol{X})]$ ,其中 $\boldsymbol{X}=\left(X_1, \cdots, X_n\right)$ .但是现在假设对于某个 $f, f(\boldsymbol{X})$ 的期望值已知(例如 $\mathrm{E}[f(\boldsymbol{X})]=\mu$ ).那么,对于任意常数 $a$ ,我们也可以用 $$ W=g(\boldsymbol{X})+a(f(\boldsymbol{X})-\mu) $$ 作为 $\mathrm{E}[g(\boldsymbol{X})]$ 的估计.现在, $$ \operatorname{Var}(W)=\operatorname{Var}(g(\boldsymbol{X}))+a^2 \operatorname{Var}(f(\boldsymbol{X}))+2 a \operatorname{Cov}(g(\boldsymbol{X}), f(\boldsymbol{X})) $$ 简单的计算说明了上式当 $$ a=\frac{-\operatorname{Cov}(f(\boldsymbol{X}), g(\boldsymbol{X}))}{\operatorname{Var}(f(\boldsymbol{X}))} $$ 时达到最小,而对于这个值 $a$ , $$ \operatorname{Var}(W)=\operatorname{Var}(g(\boldsymbol{X}))-\frac{[\operatorname{Cov}(f(\boldsymbol{X}), g(\boldsymbol{X}))]^2}{\operatorname{Var}(f(\boldsymbol{X}))} $$ 因为 $\operatorname{Var}(f(\boldsymbol{X}))$ 和 $\operatorname{Cov}(f(\boldsymbol{X}), g(\boldsymbol{X}))$ 通常是不知道的,所以应该用模拟的数据来估计这些量。

将上式除以 $\operatorname{Var}(g(\boldsymbol{X}))$ 得出 $$ \frac{\operatorname{Var}(W)}{\operatorname{Var}(g(\boldsymbol{X}))}=1-\operatorname{Corr}^2(f(\boldsymbol{X}), g(\boldsymbol{X})) $$ 其中 $\operatorname{Corr}(X, Y)$ 是 $X$ 和 $Y$ 之间的相关系数.因此,当 $f(\boldsymbol{X})$ 和 $g(\boldsymbol{X})$ 强相关时,使用控制变量将大大地降低模拟估计的方差.

MCMC

之前我们提到过贝叶斯学派的核心,MCMC和变分推断,MCMC速度慢但是精度高,变分推断速度快但是精度低,我们将详细介绍两种MCMC算法,Metropolis-Hasting算法和Gibbs算法。

Metropolis-Hasting算法

个人认为MCMC是一种很逆天的算法,已知Markov过程的稳态分布去求转移概率,在离散情况下的表现就是由一个 $n$ 维向量确定 $n^2$ 维矩阵,是很难处理的问题。

如何得到有平稳分布的转移概率矩阵?一个充分条件是转移概率满足细致平衡条件。如果存在 $\left\{\pi_j, j \in S\right\}$ , $\pi_j \geq 0, \sum_{j \in S} \pi_j=1$ ,使得 $$ \pi_i p_{i j}=\pi_j p_{j i}, \forall i \neq j $$ 称这样的马氏链为细致平衡的(detailed balanced),这时 $\left\{\pi_j\right\}$ 是 $\left\{X_t\right\}$ 的平稳分布。

设随机变量 $X$ 分布为 $\pi(x), x \in \mathcal{X}$ 。 为论述简单起见仍假设 $\mathcal{X}$ 是离散集合。算法需要一个试转移概率函数 $T(y \mid x), x, y \in \mathcal{X}$ ,满足 $0 \leq T(y \mid x) \leq 1, \sum_y T(y \mid x)=1$ ,并且 $$ T(y \mid x)>0 \Leftrightarrow T(x \mid y)>0 . $$ 算法首先从 $\mathcal{X}$ 中任意取初值 $X^{(0)}$ 。设经过 $t$ 步后算法的当前状态为 $X^{(t)}$ ,则下一步由试转移分布 $T\left(y \mid X^{(t)}\right)$ 抽取 $Y$ ,并生成 $U \sim U(0,1)$ ,然后按如下规则转移: $$ X^{(t+1)}= \begin{cases}Y & \text { if } U \leq r\left(X^{(t)}, Y\right) \\ X^{(t)} & \text { other. }\end{cases} $$ 其中(接收$y$的概率为) $$ r(x, y)=\min \left\{1, \frac{\pi(y) T(x \mid y)}{\pi(x) T(y \mid x)}\right\} $$ 在MH算法中如果取 $T(y \mid x)=T(x \mid y)$ ,则 $r(x, y)=\min \left(1, \frac{\pi(y)}{\pi(x)}\right)$ ,相应的算法称为Metropolis抽样法。如果取 $T(y \mid x)=g(y)$(不依赖于 $x$ ),则 $r(x, y)=\min \left(1, \frac{\pi(y) / g(y)}{\pi(x) / g(x)}\right)$ ,相应的算法称为Metropolis独立抽样法,和重要抽样有相似之处,试抽样分布 $g(\cdot)$ 经常取为相对重尾的分布。

在MH算法中,目标分布 $\pi(x)$ 可以用差一个常数倍的 $\tilde{\pi}(x)=C \pi(x)$ 代替,这样关于目标分布仅知道差一个常数倍的 $\tilde{\pi}(x)$ 的情形,也可以使用此算法。

下面说明MH抽样方法的合理性。 我们来验证MH抽样的转移概率 $A(x, y)=P\left(X^{(t+1)}=y \mid X^{(t)}=x\right)$ 满足细致平衡条件。 易见 $$ A(x, y)= \begin{cases}T(y \mid x) r(x, y), & y \neq x \\ T(x \mid x)+\sum_{z \neq x} T(z \mid x)[1-r(x, z)], & y=x\end{cases} $$


说明该“易见”:

情况1:$y \neq x$

  • 只有当抽样 $Y=y$ 且接受它时,才会从 $x$ 转移到 $y$ 。
  • 抽样 $Y=y$ 的概率是 $T(y \mid x)$ ,
  • 接受 $y$ 的概率是 $r(x, y)$ ,

因此, $$ A(x, y)=T(y \mid x) \cdot r(x, y), \quad y \neq x $$ 情况2:$y=x$(即"留在原地") 有两种情况会导致留在原地: 1.从试转移分布中采样到 $Y=x$ ,然后接受(总是接受); 2.从试转移分布中采样到 $Y=z \neq x$ ,但拒绝了(即 $U>r(x, z)$ ) 所以,转移回 $x$ 的概率为:

$$ A(x, x)=T(x \mid x) \cdot 1+\sum_{z \neq x} T(z \mid x) \cdot(1-r(x, z)) $$ 即: $$ A(x, x)=T(x \mid x)+\sum_{z \neq x} T(z \mid x)[1-r(x, z)] $$


于是当 $x \neq y$ 时 $$ \pi(x) A(x, y)=\pi(x) T(y \mid x) \min \left\{1, \frac{\pi(y) T(x \mid y)}{\pi(x) T(y \mid x)}\right\}=\min \{\pi(x) T(y \mid x), \pi(y) T(x \mid y)\} $$ 等式右侧关于 $x, y$ 是对称的,所以等式左侧把 $x, y$ 交换后仍相等。所以,MH构造的马氏链以 $\{\pi(x)\}$ 为平稳分布。多数情况下MH构造的马氏链也以 $\{\pi(x)\}$ 为极限分布。

从构造的 $r$ 可以看出,MH算法也可以由接收拒绝的思想去理解,并且从构造的接受拒绝策略来看,已经有了Langevin采样的思想,即向概率密度高的方向移动。

Gibbs算法

在M-H算法中是需要接收拒绝的,当采样分布高维度时,拒绝率将会无限趋近于1,这将导致我们无法进行正常的采样,为了解决该问题,我们需要使得其接受率尽可能的大,Gibbs采样算法就是可以适合高维采样的算法,可以将其视作M-H算法的特殊情况,即Proposal设置为满条件分布。

定义(满条件分布). 对于一个变量集合 $x_I=\left\{x_i: i \in I\right\}(I \subseteq\{1,2, \cdots, k\})$ 和补集变量 $x_{-I}=\left\{x_i: i \notin I\right\}$ ,满条件分布定义为: $$ p\left(x_I \mid x_{-I}\right)=\frac{p(x)}{\int p(x) d x_I} \propto p(x) $$

算法原理:将 Gibbs 看作 Metropolis-Hastings,其中旧状态是 $x=\left(x_i, x_{-i}\right)$ ,新状态是 $x^{\prime}=\left(x_i^{\prime}, x_{-i}\right)$ 。 由于仅更新了 $x_i$ ,且 proposal 是 $\pi\left(x_i^{\prime} \mid x_{-i}\right)$ ,我们写接受概率为:

$$ \alpha\left(x, x^{\prime}\right)=\min \left(1, \frac{\pi\left(x^{\prime}\right) \pi\left(x_i \mid x_{-i}\right)}{\pi(x) \pi\left(x_i^{\prime} \mid x_{-i}\right)}\right) $$ 注意:

  • $\pi(x)=\pi\left(x_i, x_{-i}\right)$
  • $\pi\left(x^{\prime}\right)=\pi\left(x_i^{\prime}, x_{-i}\right)$

因此上式变为: $$ \alpha\left(x, x^{\prime}\right)=\min \left(1, \frac{\pi\left(x_i^{\prime}, x_{-i}\right) \pi\left(x_i \mid x_{-i}\right)}{\pi\left(x_i, x_{-i}\right) \pi\left(x_i^{\prime} \mid x_{-i}\right)}\right) $$ 利用条件概率定义: $$ \begin{aligned} & \pi\left(x_i, x_{-i}\right)=\pi\left(x_i \mid x_{-i}\right) \pi\left(x_{-i}\right) \\ & \pi\left(x_i^{\prime}, x_{-i}\right)=\pi\left(x_i^{\prime} \mid x_{-i}\right) \pi\left(x_{-i}\right) \end{aligned} $$ 代入后可以发现: $$ \alpha\left(x, x^{\prime}\right)=\min \left(1, \frac{\pi\left(x_i^{\prime} \mid x_{-i}\right) \pi\left(x_{-i}\right) \pi\left(x_i \mid x_{-i}\right)}{\pi\left(x_i \mid x_{-i}\right) \pi\left(x_{-i}\right) \pi\left(x_i^{\prime} \mid x_{-i}\right)}\right)=\min (1,1)=1 $$

则Gibbs采样每步只会更新一个变量,每个循环更新所有的多维变量,接受率始终为1,但是对单个变量的采样中,可能还会需要M-H算法进行采样。注意,Gibbs算法每步只更新一个参数是重要的,只有这样才能保证细致平衡条件,也就是说,Gibbs采样是沿着坐标轴采样的。

MCMC的理论保证

MCMC算法非常的丰富,我们仅仅阐述了两种最基础的MCMC算法,有一些基于物理学思想的MCMC算法也非常有意思,例如哈密顿蒙特卡洛、朗之万蒙特卡洛,以后有机会我们再介绍,本章节我们从Markov链的理论出发,再次回顾MCMC算法,很大部分参考视频:【贝叶斯统计】MCMC理论2 Metropolis-Hastings与two-stage Gibbs的理论.

另外,本人之前阅读过关于Langevin MCMC的一些论文,可见该网页.

定义 (转移核). Let $(X, \mathcal{A})$ and $(Y, \mathcal{B})$ be measurable spaces. A Markov kernel with source $(X, \mathcal{A})$ and target $(Y, \mathcal{B})$, sometimes written as $\kappa:(X, \mathcal{A}) \rightarrow(Y, \mathcal{B})$, is a function $\kappa: \mathcal{B} \times X \rightarrow[0,1]$ with the following properties:

  1. For every (fixed) $B_0 \in \mathcal{B}$, the map $x \mapsto \kappa\left(B_0, x\right)$ is $\mathcal{A}$-measurable
  2. For every (fixed) $x_0 \in X$, the map $B \mapsto \kappa\left(B, x_0\right)$ is a probability measure on $(Y, \mathcal{B})$

In other words it associates to each point $x \in X$ a probability measure $\kappa(d y \mid x): B \mapsto \kappa(B, x)$ on $(Y, \mathcal{B})$ such that, for every measurable set $B \in \mathcal{B}$, the map $x \mapsto \kappa(B, x)$ is measurable with respect to the $\sigma$-algebra $\mathcal{A} .$

在基础的Markov链中,我们应该都接触过不可约的概念,MCMC理论中使用的不可约比常见的不可约条件稍弱:

定义 ( $\varphi$-irreducibility). $\phi$-不可约性:存在一个非零测度 $\phi$ ,使得从任意 $x$ ,对任何测度为正的集合 $A$ ,都有某个时间 $n$ 使 $k^n(x, A)>0$ ;

我们提到过细致平衡条件,以下的运算说明满足细致平衡条件的分布是稳态分布: $$ \begin{aligned} p(z) & =\sum_{z^{\prime}} p^{(\text {prev })}\left(z^{\prime}\right) T\left(z^{\prime}, z\right) \\ & =\sum_{z^{\prime}} p^{(\text {prev })}(z) T\left(z, z^{\prime}\right) \\ & =p^{(\text {prev })}(z) \sum_{z^{\prime}} T\left(z, z^{\prime}\right) \\ & =p^{(\text {prev })}(z) \sum_{z^{\prime}} p\left(z^{\prime} \mid z\right) \\ & =p^{(\text {prev })}(z) \sum_{z^{\prime}} \frac{p\left(z^{\prime}, z\right)}{p(z)} \\ & =p^{(\text {prev })}(z) \frac{p(z)}{p(z)} \\ & =p^{(\text {prev })}(z) \end{aligned} $$

在前面的算法中,我们直接使用细致平衡构造稳态分布,但是对于一个Markov链是否会收敛到稳态分布,以及稳态分布是否可以用来推断缺一无所知,我们接下来回答这两个问题。

对于连续时间Markov链,为了利于问题的分析,我们引入一个“最小单元”:
定义 (Atom). 一个集合 $A \subseteq \mathcal{X}$ 被称为一个 原子(atom),如果存在一个概率测度 $\nu$ 使得对于所有 $x \in A$ 和所有 $B \in$ $\mathcal{B}$ 有: $$ P(x, B)=\nu(B) $$ 也就是说:

  • 从集合 A 中的任何状态出发,下一步转移的分布是完全一样的,等于 $\nu$ ;
  • 因此,集合 $A$"忘记了"进入这个集合前的所有历史,是一个概率上的再生点(regeneration point)。

于是我们可以进而讨论Atom的暂态和常返态,对于MCMC,我们自然希望所有的状态都是常返态,我们将使用更强力的定义来作为以后的条件。

定义(Harris recurrent). 链 $\left\{X_n\right\}$ 被称为是 Harris recurrent,如果它是 $\phi$-不可约的,并且满足: 对于所有 $A \in \mathcal{B}$ 满足 $\phi(A)>0$ ,对所有初始状态 $x \in \mathcal{X}$ ,几乎必然无限次访问 $A$ : $$ \mathbb{P}_x\left(X_n \in A \text { infinitely often }\right)=1 $$ 换句话说: 所有测度为正的集合 $A$ ,链几乎总会返回无数次。

在 Harris recurrent 基础上,若存在不变概率测度 𝜋,则成Markov链为 Positive Harris recurrent.

定理(收敛定理). 任意Positive Harris recurrent且非周期的Markovb链,对于所有初始状态 $x \in \mathcal{X}$ ,转移概率 $P^n(x, \cdot)$ 弱收敛到不变测度 $\pi$ : $$ \lim _{n \rightarrow \infty}\left|P^n(x, \cdot)-\pi(\cdot)\right|{T V}=0 $$ 这里 $|\cdot|{T V}$ 表示 全变差距离(total variation distance)。

收敛定理表明了满足一定条件下的Markov链可以收敛到稳态分布,我们之前推导出来过M-H算法的转移核,可以验证其复合收敛定理的假设。

定理(遍历性定理). 设 $\left\{X_n\right\}$ 是 positive Harris recurrent 且 aperiodic 的 Markov 链,不变分布为 $\pi$ 。那么, 对任何 $f \in L^1(\pi)$ ,有: $$ \lim _{n \rightarrow \infty} \frac{1}{n} \sum_{k=0}^{n-1} f\left(X_k\right)=\int f d \pi \quad \text { a.s. } \quad \forall X_0=x $$

遍历性定理对标的是SLLN,但是这并不是trivial的,因为SLLN要求随机样本,是iid的,但是对于Markov链,样本是强相关的,正因为有了Markov链的遍历性定理,才让我们的采样得以发挥作用,使用生成的后验分布样本去计算一些统计量,这些统计量是一致的。

M-H算法和Gibbs算法是最简单的MCMC算法,其对假设的要求都是最弱的,也容易验证其收敛定理和遍历性定理都是成立的,但是对于一些更复杂的算法,这两个定理未必会成立,这时就需要具体问题具体分析,去检查MCMC带来的误差了。

在实际的数值实验中,MCMC会有各种各样的问题,不同分布的处理方式也多种多样,很可能会碰到MCMC不能良好采样的问题,例如当分布不一致有界时,M-H算法拒绝率可能会很高,甚至最后也不会收敛。

序贯蒙特卡洛采样(粒子滤波器)

Licensed under CC BY-NC-SA 4.0
Last updated on May 01, 2025 22:32 CST
Page views:Loading
Built with Hugo
Theme Stack designed by Jimmy