统计推断 - Stellarium

贝叶斯统计分析思维

$f(X|\theta)$ is the conditional density of $X=(X_1,X_2,...,X_n)$ given $\theta$ (likelihood);
$\tau(\theta)$ is the prior density on $\theta$
$g(X)$ is the marginal density of $X$ , i.e. $g(X)=\int_{\Theta} f(X|\theta)\tau(\theta)d\theta$
$f(X,\theta)$ is the joint density of $X$ and $\theta$

f(X,\theta)=f(X|\theta)\tau(\theta)=h(\theta|X)g(X)

$h(\theta|X)$ is the posterior density of $\theta$ given $X=(X_1,X_2,...,X_n)$

h(\theta|X)=\frac{f(X,\theta)}{g(X)}=\frac{f(X|\theta)\tau(\theta)}{\int_{\Theta}f(X|\theta)\tau(\theta)d\theta}

极大似然估计的思想

小概率原理：（小概率事件的实际不可能发生原理），小概率事件在个别试验中几乎是不可能发生的。
重要结论：如果随机事件的概率非常接近1，则可以说明在个别试验中该事件几乎一定发生
直观解释：一个随机试验如果有若干个结果 $A,B,C...$ ，如果在一次试验中结果 $A$ 出现，则一般认为试验条件对 $A$ 出现有利，也即 $A$ 出现的概率最大。

以离散型为例，若总体 $X$ 属于离散型，其分布律为 $P(X=x)=p(x;\theta_1,\theta_2,...,\theta_k)$ ， $\theta=(\theta_1,\theta_2,...,\theta_k) \in \Theta$ 的形式为已知， $\theta$ 为待估参数，属于参数空间 $\Theta$ 。

设 $(X_1,X_2,...,X_n)$ 是来自总体 $X$ 的样本，求 $(\theta_1,\theta_2,...,\theta_k)$ 的极大似然估计量。

建立似然函数 $L(\theta_1,\theta_2,...,\theta_k)=\Pi_{i=1}^{n}p(x_i;\theta_1,\theta_2,...,\theta_k)$
取对数 $lnL=\sum_{i=1}^{n}lnp(x_i;\theta_1,\theta_2,...,\theta_k)$
令对 $\theta_{i}$ 的各个偏导数等于 $0$ ，解方程组求得 $(\theta_1,\theta_2,...,\theta_k)$ 的极大似然估计值（得分函数）

费歇尔信息量

基于极大似然估计， $X_1,...,X_n$ 是一组独立同分布于 $f$ 的随机样本， $x_1,...,x_n$ 是一组观测。

似然函数（Likelihood Function）：
似然函数描述了已知一组观测数据 $x1,x2,...,xn$ 后，未知参数 $\theta$ 的可能取值。
对于独立同分布的随机样本，似然函数可以表示为：

L(\theta | \textbf{x}) = f(\textbf{x} | \theta) = \prod_{i=1}^n f(x_i | \theta)

其中， $f(xi∣\theta)$ 是给定参数 $\theta$ 下观测值 $x_i$ 的概率密度函数或概率质量函数。

对数似然函数（Log-Likelihood Function）：
对数似然函数是似然函数的自然对数。
通常，我们取对数似然函数的负值，因此实际值越大越好。
对数似然函数的一般形式为：

\ln L(\theta | \textbf{x}) = \sum_{i=1}^n \ln f(x_i | \theta)

得分函数（Score Function）：
得分函数是对数似然函数关于参数 $\theta$ 的一阶偏导数。
对于连续数据，得分函数可以表示为：

S(\theta | \textbf{x}) = \frac{\partial}{\partial \theta} \ln L(\theta | \textbf{x})

对于离散数据，得分函数可以表示为：

S(\theta | \textbf{x}) = \sum_{i=1}^n \frac{\partial}{\partial \theta} \ln f(x_i | \theta)

CRLB

不是所有的参数和分布都能达成CRLB最下限分散度估值，得分函数 $V(X,\theta)$ 一定表现为如下形式：

\frac{\alpha}{\alpha \theta}=V(X,\theta)=k_n(\theta)[W(X)-\tau(\theta)]

$k_n(\theta)$ 因子是一个概率系数，所以不包含数据，所要做的就是构造这个得分函数。

如果满足，则就是UMVUE，否则就说明CRLB unattainable

L-S theorem

所以现在想要求得UMVUE，则需要使用L-S定理，不用知道每一个样本 $X_1...X_n$ ，对于比如说正态分布， $\hat{\mu}_{MLE}=\hat{\mu}_{UMVUE}=\bar{X}$ ，仅要一个 $T(X)=\bar{X}$ 值对估计 $\mu$ 就够了，这样我们就称 $T(X)$ sufficient statistical for $\mu$

证明T是sufficient以及complete

想要证明T既是sufficient又是complete的，首先判断是否是exponential family；
如果是：将原分布拆解成如下格式，其中，函数c需要单调

f(x;\theta)=a(\theta)b(x)c(\theta)d(x)

如果不是：则使用L-S或者Neyman-Fisher先证明sufficient
L-S：列出似然函数，附上Indicator-Function，形如

L(x_i;\theta)=\Pi_{i=1}^{n}[f(x;\theta) \cdot I\{x,\theta\}]

$T=t(x)$ 是 $\theta$ 的充分统计量，当且仅当样本值的联合密度可以表示为如下格式，目标正是需要凑出这种格式

L(x_i;\theta)=g(T,\theta)h(x)

比如对于以下概率密度函数， $w=x(1)$ ， $w$ 的定义是最小样本

f(x;\theta)=\frac{\alpha\theta^{\alpha}}{x^{\alpha+1}} x\ge \theta

即可以写出

L(x_i;\theta)=\Pi_{i=1}^{n}[f(x_i;\theta)\cdot I\{x_i\ge\theta\}]

拆解可得

L(x_i;\theta)=I\{x\ge\theta\}\theta^{n\alpha}\Pi_{i=1}^{n}[\frac{\alpha}{x_i^{\alpha+1}}]\\ =g(x(1),\theta)h(x_1,x_2,...,x_n)\\

想要证明T是minimal sufficient的，如果样本只有一个维度，则可根据Neyman-Fisher来说明 $w=x(1)$ 是minimal sufficient；
如果硬是要证明，由于L-S：Consider a partition A of X by define for any $x\in{X}$ , $A(x)=\{y:\frac{L(y,\theta)}{L(x,\theta)}\}$ , the ratio is a function of the type $h(y,x)$ , the defined sets $A(x),x\in{X}$ form a partition of $X$ and this partition is minimal sufficient.
比如对于这道题，可以列如下式子：

\frac{L(y,\theta)}{L(x,\theta)}=\frac{[\Pi_{i=1}^{n}[\alpha/y_i^{\alpha+1}]I{y(1)\ge\theta}]}{[\Pi_{i=1}^{n}[\alpha/x_i^{\alpha+1}]I{x(1)\ge\theta}]}

根据上述化简可得知仅需 $x(1)=y(1)$ ，不再依赖于 $\theta$ ，最小依赖为 $T(x)=x(1)$ 得证

如何证明 $x(1)$ 是complete的，需要用到completeness的定义，假设对于任意的 $\theta>0$ ， $E[g(t)]=0$
假设单个样本的概率密度函数为 $f_{x(1)}x(1)=\frac{n\alpha\theta^{n\alpha} }{x^{n\alpha+1} }$

E(g(t))=\int_{\theta}^{\infty}g(t)\frac{n\alpha\theta^{n\alpha} }{t^{n\alpha+1} }dt=0 \\

\frac{\alpha}{\alpha\theta}E(g(t))=0\\ -g(\theta)\frac{n\alpha\theta^{n\alpha} }{\theta^{n\alpha+1} }=0

因为 $\frac{n\alpha\theta^{n\alpha} }{\theta^{n\alpha+1} }>0$ ，所以 $g(t)=0$ ，当 $t\ge\theta>0$ ，综上， $P(g(T)=0)=1$

找unbiased estimator

比如，现在有 $x_1,x_2,...,x_n$ 是一组泊松分布的样本，其中有8个样本值为 $1,4,1,0,2,1,3,2$
Find unbiased estimator for $h(\theta)=\frac{1}{2}\theta^2e^{-\theta}$ which is $p(x_1=2)$
此处不难看出所要求的参数表达式 $h(\theta)$ 是某个概率值，所以用indicator function

泊松分布的概率密度函数为 $\frac{e^{-\theta}\theta^x}{x!}$

题目这样出可以直接看出

w=I\{x_1=2\}\\ E(w)=E(I\{x_1=2\})=p(x_1=2)=\frac{1}{2}\theta^2e^{-\theta}

但如果题目问 $h(\theta)=\theta^3e^{-\theta}$

w=6I\{x_1=3\}\\ E(w)=E(6I\{x_1=3\})=6p(x_1=3)=h(\theta)

或者 $h(\theta)=\frac{1}{4}\theta^4e^{-2\theta}$ 这样的

w=4I{x_1=2,x_2=2}\\ E(w)=E(4I{x_1=2,x_2=2})=4p(x_1=2)p(x_2=2)=h(\theta)

再比如，题目问 $h(\theta)=e^{-k\theta}$ ，当 $k=1,2,...,n$ ， $p(x_1=0)=e^{-\theta}$

w=I\{x_1=0,...x_k=0\} \\ E(w)=E(I\{x_1=0,...x_k=0\}) \\ =p(x_1=0,...,x_k=0)=p(x_1=0)^k=h(\theta)

找biased estimator

取单样本（最小样本）的期望

比如，现在有 $x_1,x_2,...,x_n$ 是一组相互独立的随机变量，有概率密度函数 $f(x;\theta)=e^{-(x-\theta)}$ ，当 $x>\theta$ ， $T=min\{x_1,...,x_n\}=x(1)$ ，现在我们已知 $f_T(t)=ne^{n(\theta-t)}$ ，当 $t\ge\theta$

Show that the MLE is a biased estimator. Hint: You might want to consider the transformation $Y=T-\theta$ when performing the integral and then utilize the density of an exponential distribution.

E[x(1)]=\int_{\theta}^{\infty}ne^{-n(t-\theta)}t dt

令 $y=t-\theta$ ， $t \in [\theta,\infty)$

E[x(1)]=n\int_{0}^{\infty}e^{-ny}(y+\theta) dy\\ =n\int_{0}^{\infty}ye^{-ny}dy+n\theta\int_{0}^{\infty}e^{-ny}dy=\theta+\frac{1}{n}

不等于 $\theta$ ，所以这个的MLE是一个biased estimator

注意，这里 $\int_{0}^{\infty}ye^{-ny}dy$ 不好求解，我们可以利用Gamma函数性质，其pdf为 $f(x|\alpha,\beta)=\frac{1}{\Gamma(\alpha)\beta^{\alpha} }x^{\alpha-1}e^{-\frac{x}{\beta} }$ ，我们不难得出 $\beta=\frac{1}{n}$ ， $\alpha=2$ ，代入解得

\Gamma(2)(\frac{1}{n})^2\int_{0}^{\infty}\frac{1}{\Gamma(2)(\frac{1}{n})^2}ye^{-ny}dy \\ =(\frac{1}{n})^2

求UMVUE

首先需要T是充分且必要对于 $\theta$ ，并且已经找到了关于h的无偏估计

比如已知 $w=I\{x_1=0,...,x_k=0\}$ ，满足泊松分布， $h(\theta)=e^{-k\theta}$

E(w|T=t)=E(I\{x_1=0,...,x_k=0\}|\sum_{i=1}^nx_i=t)\\ =p(x_1=0,...,x_k=0|\sum_{i=1}^nx_i=t)\\ =\frac{p(x_1=0,...,x_k=0,\sum_{i=k+1}^nx_i=t)}{p(\sum_{i=1}^nx_i=t)}\\ =\frac{p(x_1=0,...,x_k=0)\cdot p(\sum_{i=k+1}^nx_i=t)}{p(\sum_{i=1}^nx_i=t)}\\ =\frac{e^{-k\theta}\cdot e^{-(n-k)\theta}((n-k)\theta)^x/x!}{e^{-n\theta}(n\theta)^x/x!}\\ =(1-\frac{k}{n})^{n\bar{x}}

这就是 $h(\theta)=e^{-k\theta}$ 的UMVUE

求MLE

要使用原样本的likelihood

先求得 $\theta_{MLE}$ ，再求出 $h(\theta)_{MLE}$

求MLE你首先得知道Score Function，使用基于全样本求得的一阶导为Score Function，求解过程在求CRLB中， $V(x,\theta)=-n+\frac{\sum_{i=1}^nx_i}{\theta}=0$ ，解得 $\theta=\bar{x}$ ， $\hat{h(\theta)_{MLE} }=h(\hat{\theta}_{MLE})=e^{-k\bar{x} }$

再对Score Function求一下导， $\frac{\alpha}{\alpha\theta}V(x,\theta)=-\frac{\sum_{i=1}^nx_i}{\theta^2}<0$

还有一种情况，单独只问你求 $\theta$ 的MLE，比如：对于 $L(x_i;\theta)=\Pi_{i=1}^n\frac{\theta}{ {x_i}^2} \cdot I\{x_i\ge\theta\}$

\Pi_{i=1}^n[\frac{\theta}{ {x_i}^2} \cdot I\{x_i\ge\theta\}]=\frac{\theta^n}{\Pi_{i=1}^n{x_i}^2} \cdot I\{x_1\ge\theta\}

当 $x_1<\theta$ ， $L(x_i;\theta)=0$ ， $\theta^n$ 达到最大当 $\theta=x(1)$ ， $\theta_{MLE}=x(1)$

求CRLB

比如说对于 $f(x;\theta)=\frac{e^{-\theta}\theta^x}{x!}$ ，calculate the CRLB for the minimal variance of an unbiased estimator of $h(\theta)=e^{-k\theta}$

直接面向全样本

L(x;\theta)=\Pi_{i=1}^n\frac{e^{-\theta}\theta^{x_i}}{x_i}=\frac{e^{-ne}\theta^{\sum_{i=1}^{n}x_i}}{\Pi_{i=1}^nx_i!}

lnL(x;\theta)=-n\theta+ln\theta\sum_{i=1}^n{x_i}+常数\\ \frac{\alpha}{\alpha\theta}lnL(x;\theta)=-n+\frac{\sum_{i=1}^n{x_i}}{\theta}\\ \frac{\alpha^2}{\alpha\theta^2}lnL(x;\theta)=-\frac{\sum_{i=1}^n{x_i}}{\theta^2}\\

一阶导等于 $0$ ，解得 $\theta=\bar{x}$

I_x(\theta)=-E[\frac{d^2}{d\theta^2}lnL]=\frac{E(\sum_{i=1}^n{x_i})}{\theta^2}=\frac{E(n\theta)}{\theta^2}=\frac{n}{\theta}

代入CRLB公式解得

\frac{(\frac{\alpha}{\alpha\theta}h(\theta))^2}{I_x(\theta)}=\frac{\theta}{n}k^2e^{-2k\theta}

UMVUE的方差

比如问：Show that there does not exist an integer k for which the variance of the UMVUE of $h(\theta)$ attains this bound.

使用基于全样本求得的一阶导为Score Function， $V(x,\theta)=-n+\frac{\sum_{i=1}^n{x_i}}{\theta}=-n+\frac{n\bar{x}}{\theta}$

V(x,\theta)=-n+\frac{n\bar{x}}{\theta}\\ =n(\frac{\bar{x}}{\theta}-1)

如何构造出 $V(x,\theta)=k_n(\theta)[w(x)-\tau(\theta)]$ 这样的格式

首先 $\tau(\theta)=h(\theta)=e^{-k\theta}$

V(x,\theta)=k_n(\theta)[w(x)-\tau(\theta)]\\ =\frac{n}{e^{-k\theta}}[\frac{\bar{x}e^{-k\theta}}{\theta}-e^{-k\theta}]

所以，无论选择哪个整数， $\frac{\bar{x}e^{-k\theta} }{\theta}$ 都不能作为一个统计量，因为仍然需要依赖 $\theta$

求贝叶斯估值

题目问：The prior on $\theta \in (0,1)$ is $Beta(30,3)$ ，determine the Bayesiam estimate of $\theta$ ，w.r.t. quadratic loss $Beta(\alpha,\beta)$

f(x|\alpha,\beta)=\frac{1}{Beta(\alpha,\beta)} \cdot x^{\alpha-1}(1-x)^{\beta-1}

其中 $Beta(\alpha,\beta)$ 是常数，已经求得prior $\tau(x;\theta)$ 为 $\theta^{\sum_{i=1}^nx_i}(1-\theta)^n$ ，也就是似然函数，结果是 $\theta^{136}(1-\theta)^6$

posterior则正比于 $\tau(x;\theta)\cdot f(x|\alpha,\beta)=\theta^{165}(1-\theta)^8~Beta(166,9)$

\hat{\theta}=E[\theta|x]=\frac{\alpha}{\alpha+\beta}=\frac{166}{175}=0.9486

以上求期望参考Beta分布的期望公式

决策论

Consider testing $H_0:\theta \ge 0.95$ versus $H_1:\theta < 0.95$ with a 0-1 loss in a Bayesiam setting with the above $Beta(30,3)$ prior, what is your decision? You may use: $\int_{0.95}^1 \cdot x^{165}(1-x)^8dx \cdot 174 \cdot 173 ... \cdot 166=20428.19$