Box-Cox变换

1. Box-Cox变换

如果需要对一个变量进行变换,使之成为正态成布,可以使用Box-Cox变换。

Box-Cox变换对因变量 Y 作如下变换: math Y^{(\lambda)}=\begin{cases} \frac{Y^\lambda-1}{\lambda},\lambda \neq0 \\ ln Y, \lambda=0 \end{cases} 这里$\lambda$是一个待定变换参数。对不同的$\lambda$,所做的变换自然就不同,所以是一个变换族。它包括了对数变换($\lambda=0$),平方根变换($\lambda=\frac{1}{2}$)和倒数变换($\lambda=-1$)等常用变换。

Image Alt
变换前变量的分布
Image Alt
变换后变量分布

对因变量的n个观测值$y1,y2,...yn$,应用上述变换,得到变换后的向量$y^{(\lambda)}=(y1^{\lambda},y2^{\lambda},...yn^{\lambda})$。 即要确定变换参数$\lambda$,使得$y^{(\lambda)}$满足 math Y^{(\lambda)}=X\beta+e,e\sim N(0,\sigma^2I) 也就是说,通过对因变量的变换,使得变换过的向量$y^{(\lambda)}$与回归自变量具有线性相依关系,误差也服从正态分布,误差各分量是等方差且相互独立。 以极大似然法来确定$\lambda$。因为$Y^{(\lambda)}\sim N(X\beta,\sigma^2I)$,所以对固定的$\lambda,\beta,\sigma^2$的似然函数为 math L(\beta,\sigma^2)=\frac{1}{(\sqrt {2\pi})^2}exp\{ -\frac{1}{2\sigma ^2}(y^{(\lambda)}-X\beta)'(y^{(\lambda)}-X\beta)\} J 这里 J 为变换 Jacobi 的行列式 math J=\prod_{i=1}^n|\frac{dy_i^{(\lambda)}}{dy_i}|=\prod_{i=1}^ny_i^{\lambda-1} 当$\lambda$固定时,$J$是不依赖于参数$\beta$和$\sigma^2$的常数因子。$L(\beta,\sigma^2)$的其余部分关于$\beta$和$\sigma^2$求导数,令其等于0,可以求得$\beta$和$\sigma^2$的极大似然估计 math \hat{\beta}(\lambda)=(X'X)^{-1}X'y^{(\lambda)} \\ \hat{\sigma}^2(\lambda)=\frac{1}{n}y^{(\lambda)}(I-X(X'X)^{-1}X')y^{(\lambda)}=\frac{1}{n}RSS(\lambda,y^{(\lambda)}) 其中 math RSS(\lambda,z^{(\lambda)})=z^{(\lambda)}(I-X(X'X)^{-1}X')z^{(\lambda)} ...(*) 这里,$z^{(\lambda)}=(z1^{(\lambda)},...,zn^{(\lambda)})'=\frac{y^{(\lambda)}}{J^{\frac{1}{n}}}$ math z_i^{(\lambda)}=\begin{cases} \frac{y_i^\lambda}{(\prod_{i=1}^ny_i)^{\frac{\lambda-1}{n}}}, \lambda \neq 0 \\ (lny_i)(\prod_{i=1}^ny_i)^{\frac{1}{n}},\lambda=0 \end{cases} (*)式对Box-Cox变换带来很大方便,因为为了求$lnL_{max}(\lambda)$的最大值,只需求残差平方和的$RSS(\lambda,y^{(\lambda)})$最小值。

2.单变量的Box-Cox变换

设变量$y$经变换后,

$y^{(\lambda)}\sim N(\mu,\sigma^2)$

对固定的$\lambda,\beta,\sigma^2$的似然函数为 math L(\beta,\sigma^2)=\frac{1}{(\sqrt{2\pi}\sigma^2)^n} exp\{-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i^{(\lambda)}-\mu)^2 \} J $J$同为变换 Jacobi 的行列式 math J=\prod_{i=1}^n|\frac{dy_i^{(\lambda)}}{dy_i}|=\lambda^n \prod_{i=1}^ny_i^{\lambda-1} 求得$\beta$和$\sigma^2$的极大似然估计为

$\hat\mu=\overline y^{(\lambda)}=\frac{1}{n}\sum{i=1}^nyi^{(\lambda)}$

$\hat\sigma^2=\frac{1}{n}\sum{i=1}^{n}(yi^{(\lambda)}-\hat\mu)^2$

对极大似然函数作对数变换 math lnL_{max}(\lambda)=-\frac{n}{2}-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\hat\sigma^2)+lnJ 化简得 math lnL_{max}(\lambda)=-\frac{n}{2}-\frac{n}{2}ln(2\pi)- \frac{n}{2}ln(\frac{1}{n}\sum_{i=1}^n(z_i^{(\lambda)}-\overline z^{(\lambda)})^2) 其中

$z^{(\lambda)}=gm(y^{(\lambda)})^{1-\lambda}\cdot y^{(\lambda)}$

$gm(y^{(\lambda)})=exp(\sum{i=1}^n \frac{1}{n}ln(yi^{(\lambda)}))$

为了简单起见,重新将 Box-Cox 变换定义为 math Y^{(\lambda)}=\begin{cases} gm(Y)^{1-\lambda} \cdot \frac{Y^\lambda-1}{\lambda},\lambda \neq0 \\ gm(Y), \lambda=0 \end{cases} 为了最大化$lnL_{max}(\lambda)$,只须最小化$\frac{n}{2}ln(\frac{n-1}{n}var(y^{(\lambda)}))$。

3. 黄金分割搜索法

黄金分割法(Golden Section Method),是用于在单峰函数区间上求极小值的一种方法。其基本思想是通过取试探点和函数值比较,使包含极小点的搜索区间不断减少,当区间长度缩短到一定程度时,就得到函数极小点的近似值。

设$\tau$是一元二次方程$\tau^2+\tau-1=0$的正根,即$\tau=\frac{\sqrt5-1}{2}$。

对于函数$\phi(x)$,先在搜索区间[a,b]上确定两个试探点,其中左试探点为$al=a+(1-\tau)(b-a)$,右试探点为$ar=a+\tau(b-a)$。

再分别计算这两个试探点的函数值$\phil=\phi(al)$,$\phir=\phi(ar)$。由单峰函数的性质,若$\phil<\phir$,则区间$[ar,b]$内不可能有极小点,因此去掉区间$[ar,b]$,令$a'=a,b'=ar$,得到一个新的搜索区间。若$\phil>\phir$,则区间$[a,al]$内不可能有极小点,去掉区间$[a,al]$,令$a'=al,b'=b$,得到一个新的搜索区间。

类似上面的步骤,在区间$[a',b']$内再计算两个新的试探点

$a_l'=a'+(1-\tau)(b'-a')$

$a_r'=a'+\tau(b'-a')$

比较函数值,得到新的区间。

在上述方中,事实上每次迭代并不需要计算两个试探点及函数值。下面对新的试探点进行分析。 (1) 若$\phil<\phir$,则去掉区间$[a_r,b]$,那么新的右试探点为

$ar'=a'+\tau(b'-a')=a+\tau(ar-a)=a+\tau^2(b-a)$

注意到$\tau$是方程$\tau^2+\tau-1=0$的根,因此有

$ar'=a+\tau^2(b-a)=a+(1-\tau)(b-a)=al$

即原区间的左试探点。

(2) 若$\phil>\phir$,则去掉区间$[a,a_l]$,那么新的左试探点为

$al'=a'+(1-\tau)(b'-a')=al+(1-\tau)(b-a_l) $

$=a+(1-\tau)(b-a)+\tau(1-\tau)(b-a) $

$=a+(1-\tau^2)(b-a)=a+\tau(b-a) $

$=a_r$

即原区间的右试探点。

  因此在上述计算过程中,只需要计算一个新试探点和一个点的函数值。

算法:

(1) 置初始搜索区间$[a,b]$,并置精度要求$\epsilon$,并计算左右试探点

$al=a+(1-\tau)(b-a),ar=a+\tau(b-a)$,其中$\tau=\frac{\sqrt5-1}{2}$

及相应的函数值$\phil=\phi(al)$,$\phir=\phi(ar)$。 (2) 如果$\phil<\phir$,则置

$b=ar, ar=al, \phir=\phi_l$

并计算$al=a+(1-\tau)(b-a), \phil=\phi(a_l)$。

否则

$a=al,al=ar,\phil=\phi_r$

并计算$ar=a+\tau(b-a), \phir=\phi(a_r)$。

(3)若$|b-a| \leq \epsilon$。如果$\phil<\phir$,则置问题的解$\mu=al$;否则置$\mu=ar$,停止计算。否则转到(2)继续计算。

4. 正态分布检验

对变量进行变换后,需要检验变量是否服丛正态分布。这里用到正态分布检验。

1. W检验

W检验是S.S.Shapiro和M.B.Wilk1965年提出来的,这种方法在样本容量$3\leq50$时适用。

W检验即检验假设$H_0$:总体服从正态分布

利用W检验的方法检验原假设$H_0$的步骤如下

(1)把n个样本观测值按由小到大的次序排列成 math x_{(1)}\leq x_{(2)}\leq ...x_{(i)}\leq ...\leq x_{(n)}

(2)W检验的统计量为 math W=\frac{(\sum_{i=1}^{[\frac{n}{2}]}a_i(w)(x_{(n+1-i)}-x_{(i)}))^2}{\sum_{i=1}^n(x_{(i)}-\overline x)^2} 其中$\overline x$表示样本均值,$a_i(w)$的值可查表得。$[\frac{n}{2}]$表示数$\frac{n}{2}$的整数部分。

$x_{(1)},x_{(2)}, ...x_{(n)}$的值代入上式计算统计量W的值。

(3) 根据给定的检验水平$\alpha$和样本容量n查表得统计量W的$\alpha$的分位数$W_{\alpha}$。

(4) 作出判断:

若$W<W{\alpha}$,则拒绝$H0$,认为总体不服从正态分布;

若$W \ge W{\alpha}$,则不拒绝$H0$。

2. D检验

W检验是一种有效的正态性检验方法,可惜它只适用于容量为3至50的样本。1971年D’Agostino提出了D’Agostino检验(简称D检验)。这种检验不需要附系数表,它所适用的样本容量n的范围为$50 \le n \ \le1000$。

  进行D检验的步骤如下:

(1) 把n个样本观测值按由小到大的次序排列成

$$x_{(1)}\leq x_{(2)}\leq ...x_{(i)}\leq ...\leq x_{(n)} $$

(2) D检验的统计量为 math Y=\frac{\sqrt n(D-0.28209479)}{0.02998598} 其中 math D=\frac{\sum_{i=1}^n(i-\frac{n+1}{2})x_{(i)}}{(\sqrt n)^3\sqrt {\sum_{i=1}^n(x_{(i)}-\overline x)^2}}

按上两式计算统计量Y的值。

(3) 根据给定的检验水平$\alpha$和样本容量n查表,得统计量$Y$$\frac{\alpha}{2}$分位数$Y_{\frac{\alpha}{2}}$$1-\frac{\alpha}{2}$分位数$Y_{1-\frac{\alpha}{2}}$

(4) 作出判断:若$Y \lt Y_{\frac{\alpha}{2}}$$Y \gt Y_{1-\frac{\alpha}{2}}$,则拒绝$H_0$,否则不拒绝$H_0$

附:Box-Cox变换的R代码 R BoxCox_Trans<-function(x,interval,loop=1000,epsilon= .Machine$double.eps){ Likelihood_Log<-function(x,lambda){ y_lambda<-function(x,lambda){ gm<-exp(mean(log(x))) if (lambda == 0) log(x) * gm else (gm^(1 - lambda)) * ((x^lambda) - 1)/lambda } y <- y_lambda(x, lambda) (length(y)/2) * log(((length(y) - 1)/length(y)) * var(y)) } GoldenSecSearch<-function(f,x,interval,loop, epsilon) { t<-(sqrt(5) - 1)/2 a<-min(interval) b<-max(interval) a_l<-a+(1-t)*(b-a) a_r<-a+t*(b-a) f_l<-f(x,a_l) f_r<-f(x,a_r) i<-1 while(abs(b-a)>epsilon){ i<-i+1 if(f_l<f_r){ b<-a_r a_r<-a_l f_r<-f_l a_l<-a+(1-t)*(b-a) f_l<-f(x,a_l) } else { a<-a_l a_l<-a_r f_l<-f_r a_r<-a+t*(b-a) f_r<-f(x,a_r) } if(i>loop) break } Result<-list() if(f_l<f_r) { Result$minimum<-a_l Result$Objective<-f_l } else { Result$minimum<-a_r Result$Objective<-f_r } Result } Output<-list() Output<- GoldenSecSearch(Likelihood_Log,x,interval,loop,epsilon) Output } 进行变换的R代码 ```R

attach(Prestige) //car package hist(income) BoxCox_Trans(income,seq(-3,3)) $minimum [1] 0.1792894 $Objective [1] 827.9459 hist(income^0.1792894) ```