03-参数估计.Rmd

# 参数估计

对于参数估计的模拟,主要分成三步:

1. 生成数据;
2. 估计参数;
3. 评价估计好坏.

对于生成数据,在阅读文献时,文章中会明确给出如何产生.我们不必多费周折.在设计我们自己的模拟实验时,多借鉴其他文献中的例子.一方面方便和其他文献进行比较,另一方面,可以避免落入陷阱[^1].

对于评价估计好坏,有限维离散参数一般采取$\left\|\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\right\|,$无穷维函数一般采取$\int\left(\hat{m}(x)-m_{0}(x)\right)^{2} d x.$

估计量通过解估计方程得到,根据方程的形式,分为M估计量和Z估计量,下面分别说明.

## M估计

定义:极大化(或极小化)目标函数得到参数估计值.

$$
M_{n}(\theta)=\frac{1}{n} \sum_{i=1}^{n} m_{\theta}\left(X_{i}\right)
$$

$$
\hat{\theta}=\arg \max_{\theta \in \Theta} M_{n}(\theta)
$$

其中$m_{\theta}\left(X_{i}\right)$为已知函数.特别的,如果$M_n(\theta)$可导,M估计量和Z估计量有等价形式.

下面以线性模型为例:$\boldsymbol{Y}=\boldsymbol{X}^{T}\boldsymbol{\theta}+\boldsymbol{\varepsilon}$

考虑最小二乘估计,取$m_{\theta}\left(\boldsymbol{X}_{i}\right)=-\left(Y_{i}-\boldsymbol{X}_{i}^{T} \boldsymbol{\theta}\right)^{2}$[^2],则




$$ 
\begin{aligned} 
M_{n}(\theta)&=-\frac{1}{n} \sum_{i=1}^{n}\left(Y_{i}-\boldsymbol{X}_{i}^{T} \boldsymbol{\theta}\right)^{2}\\
&=-\frac{1}{n}\left(\boldsymbol{Y}-\boldsymbol{X}^{T} \boldsymbol{\theta}\right)^{T}\left(\boldsymbol{Y}-\boldsymbol{X}^{T} \boldsymbol{\theta}\right)\\
&=-\frac{1}{n}\left(\boldsymbol{\theta}^{T} \boldsymbol{X} \boldsymbol{X}^{T} \boldsymbol{\theta}-2 \boldsymbol{Y}^{T} \boldsymbol{X}^{T} \boldsymbol{\theta}+\boldsymbol{Y}^{T} \boldsymbol{Y}\right)
\end{aligned} 
$$ 
其中最后一项是与$\boldsymbol{\theta}$无关的常数项,可以不考虑,进而可以整理成如下的二次规划问题:

利用quadprog包中的函数可以求解

```{r}
library(quadprog)
n=100;p=3;beta=c(1,-2,3);
X = matrix(rnorm(n*p),n,p)
Y = X%*%beta + rnorm(n)
lm(Y~X+0)$coef == solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution
```

可以看到结果显示不相等,但是如果我们打印出来显示:
```{r}
lm(Y~X+0)$coef
solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution
```
可以看到结果是一致的.这是由于计算机存储数字精度引起的.比如我们查看
```{r}
sum(abs(lm(Y~X+0)$coef -
solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution))
```


另外,我们可以用all.equal这个函数设置容忍的误差值,判断近似相等:
```{r}
all.equal(unname(lm(Y~X+0)$coef),
solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution
,tolerance=1e-10)
all.equal(unname(lm(Y~X+0)$coef),
solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution
,tolerance=1e-20)
```

下面考虑稍复杂的情况,取
$m_{\theta}\left(\boldsymbol{X}_{i}\right)=\rho_{\tau}\left(y_{i}-\boldsymbol{X}_{i}^{T} \boldsymbol{\theta}\right)$,其中
$\rho_{\tau}(t)=t\left(\tau-I_{\{t<0\}}\right)$.
这称为分位数回归,详细的推导过程可以在[分位数回归总结](https://github.com/dujiangbjut/dujiangbjut.github.io/tree/master/讨论班/分位数回归)查看[^3].


__这里缺少一个M估计迭代求解的例子__

## Z估计

定义:解一个等于0的方程得到参数估计.

$$
\Psi_{n}(\theta)=\frac{1}{n} \sum_{i=1}^{n} \psi_{\theta}\left(X_{i}\right)=0,
$$
其中$\psi_{\theta}\left(X_{i}\right)$为已知函数.

$\hat{\theta}$为
$\Psi_{n}(\theta)=0$的解.

回到线性模型最小二乘的例子,由于目标函数存在导数,可以转化成一个Z估计量.取
$\psi_{\theta}\left(X_{i}\right)=m_{\theta}'\left(X_{i}\right)=2\boldsymbol{X}_{i}^{T}  \left(Y_{i}-\boldsymbol{X}_{i}^{T} \boldsymbol{\theta}\right)$

则[^4]
$$
\Psi_{n}(\theta)=\frac{1}{n} \sum_{i=1}^{n}\boldsymbol{X}_{i}^{T}  \left(Y_{i}-\boldsymbol{X}_{i}^{T} \boldsymbol{\theta}\right) =\boldsymbol{X} \boldsymbol{X}^{T} \boldsymbol{\theta}-\boldsymbol{X} \boldsymbol{Y}=0,
$$

进而得到参数估计值为:$\hat{\boldsymbol{\theta}}=\left(\boldsymbol{X} \boldsymbol{X}^{T}\right)^{-1} \boldsymbol{X} \boldsymbol{Y}.$

通过下面的程序验证:

```{r}
library(quadprog)
all.equal(solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution, as.numeric(solve(t(X)%*%X)%*%t(X)%*%Y))
```


下面是一个迭代求解Z估计量的例子[^5].

**补充哪一篇参考文献**

```{r}
source("code/glm.R")
n=500
m1=2
m2=2
m3=2
m=m1+m2+m3
beta=c(rep(-1,m1),rep(0,m2),rep(1,m3))
X=runif(n*m,-1,1)
X=matrix(X,n,m)
eta=X%*%beta
mu=1/(1+exp(-eta))

# Y~B(1,p)
Y=runif(n)
Y[Y>=mu]=0
Y[Y>0]=1
glm(Y~X+0,family=binomial(link="logit"))
myglm(Y,X,distribution  = "binom")
```
可以看到和真值相差不大,但是很快收敛.



[^1]:给定模型和估计方法,不同数据的估计效果是不一样的.我们要清楚自己的方法对于什么样的数据应该有较好的结果,什么样的数据可能估计效果不好.对此设计不同代表性的实验.这也给阅读文献模拟部分添加了任务,思考作者为什么用这种方式产生数据、设置参数.
[^2]:添加负号与求极大值对应.
[^3]:后续重新整理成Rmd格式.
[^4]:省略了常数系数.
[^5]:目前只有单次二项分布的程序正确.完整的程序可以从[这里](code/glm.R)下载.