注册 登录  
 加关注

网易博客网站关停、迁移的公告:

将从2018年11月30日00:00起正式停止网易博客运营
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

数据挖掘

学习数据挖掘

 
 
 

日志

 
 

Box-Cox变换  

2013-07-05 15:25:19|  分类: SAS |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

BoxCox变换

在回归模型号中,Box-Cox变换是对因变量Y作如下变换:

clip_image002            (1.1)

这里clip_image004是一个待定变换参数。对不同的clip_image004[1],所做的变换自然就不同,所以是一个变换族。它包括了对数变换(clip_image004[2]=0),平方根变换(clip_image008)和倒数变换(clip_image004[3]=-1)等常用变换。

clip_image011

图1. 变换前变量的分布

clip_image013

图2.变换后变量分布

对因变量的n个观测值clip_image015,应用上述变换,得到变换后的向量

clip_image017          (1.2)

即要确定变换参数clip_image004[4],使得clip_image019满足

clip_image021         (1.3)

也就是说,通过对因变量的变换,使得变换过的向量clip_image019[1]与回归自变量具有线性相依关系,误差也服从正态分布,误差各分量是等方差且相互独立。

  以极大似然法来确定clip_image004[5]。因为clip_image023,所以对固定的clip_image004[6]clip_image025clip_image027的似然函数为

clip_image029 (1.4)

这里clip_image031为变换Jacobi的行列式

clip_image033           (1.5)

clip_image004[7]固定时,clip_image031[1]是不依赖于参数clip_image025[1]clip_image027[1]的常数因子。clip_image035的其余部分关于clip_image025[2]clip_image027[2]求导数,令其等于0,可以求得clip_image025[3]clip_image027[3]的极大似然估计

clip_image037          (1.6)

clip_image039 (1.7)

为了求clip_image041的最大值,考虑到lnx是x的单调函数,对clip_image041[1]求对数。略去与clip_image004[8]无关的常数项,得到

clip_image044 (1.8)

其中

clip_image046 (1.9)

clip_image048 (1.10)

clip_image050 (1.11)

(1.9)式对Box-Cox变换带来很大方便,因为为了求clip_image052的最大值,只需求残差平方和的clip_image054最小值。

单变量的Box-Cox变换

设变量clip_image056经变换后,

clip_image058 (2.1)

对固定的clip_image004[9]clip_image025[4]clip_image027[4]的似然函数为

clip_image060 (2.2)

clip_image031[2]同为变换Jacobi的行列式

clip_image062 (2.3)

  求得clip_image025[5]clip_image027[5]的极大似然估计为

clip_image064           (2.4)

clip_image066           (2.5)

对极大似然函数作对数变换

clip_image068 (2.6)

化简得

clip_image070 (2.7)

其中

clip_image072 (2.8)

clip_image074 (2.9)

(2.9)亦即为几何平均值。

为了简单起见,重新将Box-Cox变换定义为

clip_image076 (2.10)

为了最大化clip_image052[1],只须最小化clip_image079

黄金分割搜索法

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

  设clip_image081是一元二次方程

clip_image083 (3.1)

的正根,即clip_image085

  对于函数clip_image087,先在搜索区间[a,b]上确定两个试探点,其中左试探点为

clip_image089 (3.2)

右试探点为

clip_image091 (3.3)

再分别计算这两个试探点的函数值clip_image093clip_image095。由单峰函数的性质,若clip_image097,则区间clip_image099内不可能有极小点,因此去掉区间clip_image099[1],令a’=a,b’=clip_image102,得到一个新的搜索区间。若clip_image104,则区间clip_image106内不可能有极小点,去掉区间clip_image106[1],令a’=clip_image108,b’=b,得到一个新的搜索区间。

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

clip_image110 (3.4)

clip_image112 (3.5)

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

  在上述方中,事实上每次迭代并不需要计算两个试探点及函数值。下面对新的试探点进行分析。

(1) 若clip_image097[1],则去掉区间clip_image099[2],那么新的右试探点为

clip_image114 (3.6)

注意到clip_image081[1]是方程(3.1)的根,因此有

clip_image116 (3.7)

即原区间的左试探点。

(2) 若clip_image104[1],则去掉区间clip_image106[2],那么新的左试探点为

clip_image118 (3.8)

即原区间的右试探点。

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

算法:

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

clip_image089[1]clip_image091[1],其中clip_image085[1]

及相应的函数值clip_image093[1]clip_image095[1]

(2) 如果clip_image097[2],则置

b=clip_image102[1],clip_image102[2]=clip_image108[1],clip_image125,

并计算

clip_image089[2]clip_image093[2]

否则

a=clip_image108[2],clip_image129,clip_image131

并计算

clip_image091[2]clip_image095[2]

(3) 若|b-a|clip_image134,如果clip_image097[3],则置问题的解clip_image136;否则置clip_image138,停止计算。否解转到(2)继续计算。

正态分布检验

I. W检验

W检验是S.S.Shapiro和M.B.Wilk1965年提出来的,这种方法在样本容量3clip_image140nclip_image140[1]50时适用。

  W检验即检验假设

clip_image143:总体服从正态分布

  利用W检验的方法检验原假设clip_image143[1]的步骤如下

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

clip_image146

(2) W检验的统计量为

clip_image148 (4.1)

其中clip_image150表示样本均值,clip_image152的值可查表得。clip_image154表示数clip_image156的整数部分。

clip_image158的值代入(3.1)式计算统计量W的值。

(3) 根据给定的检验水平clip_image160和样本容量n查表得统计量W的clip_image160[1]的分位数clip_image163

(4) 作出间判断:若W<clip_image163[1],则拒绝clip_image143[2],认为总体不服从正态分布;若Wclip_image166clip_image163[2],则不拒绝clip_image143[3]

II. D检验

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

  进行D检验的步骤如下:

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

clip_image146[1]

(2) D检验的统计量为

clip_image170 (4.2)

其中

clip_image172 (4.3)

按(4.2)和(4.3)式计算统计量Y的值。

(3) 根据给定的检验水平clip_image160[2]和样本容量n查表,得统计量Y的clip_image175分位数clip_image177和1-clip_image175[1]分位数clip_image179

(4) 作出判断:若Y<clip_image177[1]或Y>clip_image179[1],则拒绝clip_image143[4],否则不拒绝clip_image143[5]

附:Box-Cox变换的SAS代码

%macro boxcox(phimin=,phimax=, steps=, model=, class=, stmts=, data=, response=);

/***********************************************************/

/* This macro performs thr Box-Cox-Power Transformation */

/* Box, G. E. P. and Cox, D. R. (1964) An analysis of */

/* transformations (with discussion). */

/* J. R. Statist. Soc. B., 26, 211-246. */

/* */

/* phimin = smallest phi */

/* phimax = largest phi */

/* steps = number of steps in grid search */

/* class = list of variables for class statment */

/* model = list of terms in model statement */

/* stmts = %str(other mixed statements, incl. random) */

/* data = dataset with original data */

/* response = name of response variable in dataset */

/* */

/***********************************************************/

data t;

set &data;

phi=&phimin;

size=(&phimax-&phimin)/&steps;

do i=0 to &steps;

if phi ne 0 then _y_=(&response**phi-1)/phi;

if abs(phi)<1e-10 then _y_=log(&response);

jac=(phi-1)*log(&response);

output;

phi=phi+size;

end;

proc sort data=t out=t;

by phi;

run;

ods output fitstatistics=fit;

proc mixed data=t method=ml;

class &class;

model _y_=&model;

&stmts;

by phi;

run;

data fit2;

set fit;

if descr='-2 Log Likelihood';

loglik=-value/2;

proc means data=t;

var jac;

output out=jac sum=;

by phi;

run;

data loglik;

merge jac fit2;

loglik=loglik+jac;

proc print data=loglik;

run;

symbol value=dot i=join;

proc gplot data=loglik;

plot loglik*phi;

run;

proc means data=loglik;

var loglik;

id phi;

output out=max max=max;

run;

data r;

if _n_=1 then set max;

set loglik;

if loglik=max;

title 'Optimal value of phi';

proc print data=r;

var phi loglik;

run;

title ' ';

%mend;


原文链接:http://www.cnblogs.com/zgw21cn/archive/2008/08/29/1279681.html

  评论这张
 
阅读(2308)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018