第 2 章 参数估计
对于参数估计的模拟,主要分成三步:
- 生成数据;
- 估计参数;
- 评价估计好坏.
对于生成数据,在阅读文献时,文章中会明确给出如何产生.我们不必多费周折.在设计我们自己的模拟实验时,多借鉴其他文献中的例子.一方面方便和其他文献进行比较,另一方面,可以避免落入陷阱6.
对于评价估计好坏,有限维离散参数一般采取\(\left\|\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_{0}\right\|,\)无穷维函数一般采取\(\int\left(\hat{m}(x)-m_{0}(x)\right)^{2} d x.\)
估计量通过解估计方程得到,根据方程的形式,分为M估计量和Z估计量,下面分别说明.
2.1 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}\)7,则
\[ \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包中的函数可以求解
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## X1 X2 X3
## FALSE FALSE FALSE可以看到结果显示不相等,但是如果我们打印出来显示:
lm(Y~X+0)$coef## X1 X2 X3
## 1.122 -2.070 2.893solve.QP(Dmat = t(X)%*%X/n, dvec = (t(Y)%*%X)/n, Amat = matrix(0,p,p))$solution## [1] 1.122 -2.070 2.893可以看到结果是一致的.这是由于计算机存储数字精度引起的.比如我们查看
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))## [1] 3.109e-15另外,我们可以用all.equal这个函数设置容忍的误差值,判断近似相等:
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)## [1] TRUEall.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)## [1] "Mean relative difference: 5.109e-16"下面考虑稍复杂的情况,取 \(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)\). 这称为分位数回归,详细的推导过程可以在分位数回归总结查看8.
这里缺少一个M估计迭代求解的例子
2.2 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)\)
则9 \[ \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}.\)
通过下面的程序验证:
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))## [1] TRUE下面是一个迭代求解Z估计量的例子10.
补充哪一篇参考文献
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"))##
## Call: glm(formula = Y ~ X + 0, family = binomial(link = "logit"))
##
## Coefficients:
## X1 X2 X3 X4 X5 X6
## -1.030 -0.762 0.453 -0.304 0.823 0.958
##
## Degrees of Freedom: 500 Total (i.e. Null); 494 Residual
## Null Deviance: 693
## Residual Deviance: 583 AIC: 595myglm(Y,X,distribution = "binom")## Loading required package: MASS## $theta
## [,1]
## [1,] -0.9946
## [2,] -0.7397
## [3,] 0.4485
## [4,] -0.2557
## [5,] 0.8269
## [6,] 0.8908
##
## $steps
## [1] 4可以看到和真值相差不大,但是很快收敛.