黄金分割法python 黄金分割法求极小点

1 Box－Cox变换

在回归模型号中，Box－Cox变换是对因变量Y作如下变换：

　　　　　　　　　　　　（1.1）

这里是一个待定变换参数。对不同的，所做的变换自然就不同，所以是一个变换族。它包括了对数变换（＝0），平方根变换（

）和倒数变换（＝-１）等常用变换。

图１.　变换前变量的分布

图２.变换后变量分布

对因变量的n个观测值

，应用上述变换，得到变换后的向量

　　　　　　　　　　（1.2）

即要确定变换参数，使得

满足

　　　　　　　　　（1.3）

也就是说，通过对因变量的变换，使得变换过的向量

与回归自变量具有线性相依关系，误差也服从正态分布，误差各分量是等方差且相互独立。

　　以极大似然法来确定。因为

，所以对固定的，，的似然函数为

　（1.4）

这里为变换Jacobi的行列式

　　　　　　　　　　　（1.5）

当固定时，是不依赖于参数和的常数因子。

的其余部分关于和求导数，令其等于0，可以求得和的极大似然估计

　　　　　　　　　　（1.6）

　（1.7）

为了求

的最大值，考虑到lnx是x的单调函数，对

求对数。略去与无关的常数项，得到

　（1.8）

其中

(1.9)

(1.10)

(1.1１)

(1.9)式对Box－Cox变换带来很大方便，因为为了求

的最大值，只需求残差平方和的

最小值。

2 单变量的Box-Cox变换

设变量经变换后，

(2.1)

对固定的，，的似然函数为

(2.2)

同为变换Jacobi的行列式

(2.3)

　　求得和的极大似然估计为

　　　　　　　　　　　（2.4）

　　　　　　　　　　　（2.5）

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

(2.6)

化简得

(2.7)

其中

(2.8)

(2.9)

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

为了简单起见，重新将Box－Cox变换定义为

(2.10)

为了最大化

，只须最小化

。

3 黄金分割搜索法

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

　　设是一元二次方程

(3.1)

的正根，即

。

　　对于函数

，先在搜索区间[a,b]上确定两个试探点，其中左试探点为

(3.2)

右试探点为

(3.3)

再分别计算这两个试探点的函数值

，

。由单峰函数的性质，若

，则区间

内不可能有极小点，因此去掉区间

，令a’=a,b’=,得到一个新的搜索区间。若

，则区间

内不可能有极小点，去掉区间

，令a’=,b’=b,得到一个新的搜索区间。

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

(3.4)

(3.5)

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

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

（１） 若

，则去掉区间

,那么新的右试探点为

(3.6)

注意到是方程(3.1)的根，因此有

(3.7)

即原区间的左试探点。

（２） 若

，则去掉区间

，那么新的左试探点为

(3.8)

即原区间的右试探点。

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

算法：

（１） 置初始搜索区间[a,b],并置精度要求,并计算左右试探点

，

，其中

，

及相应的函数值

，

。

（２） 如果

，则置

b=,=,

,

并计算

，

否则

a=,

,

并计算

，

（３） 若|b-a|,如果

，则置问题的解

；否则置

，停止计算。否解转到（２）继续计算。

4 正态分布检验

I. Ｗ检验

Ｗ检验是S.S.Shapiro和M.B.Wilk1965年提出来的，这种方法在样本容量3n50时适用。

　　Ｗ检验即检验假设

：总体服从正态分布

　　利用Ｗ检验的方法检验原假设的步骤如下

（１） 把n个样本观测值按由小到大的次序排列成

（２） Ｗ检验的统计量为

(4.1)

其中表示样本均值，

的值可查表得。表示数的整数部分。

将

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

（３） 根据给定的检验水平和样本容量n查表得统计量Ｗ的的分位数。

（４） 作出间判断：若Ｗ<，则拒绝,认为总体不服从正态分布；若W，则不拒绝。

II. Ｄ检验

Ｗ检验是一种有效的正态性检验方法，可惜它只适用于容量为3至50的样本。1971年D’Agostino提出了D’Agostino检验（简称Ｄ检验）。这种检验不需要附系数表，它所适用的样本容量n的范围为50n1000。

　　进行Ｄ检验的步骤如下：

（１） 把n个样本观测值按由小到大的次序排列成

（２） Ｄ检验的统计量为

(4.2)

其中

(4.3)

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

（３） 根据给定的检验水平和样本容量n查表，得统计量Ｙ的分位数和１-分位数；

（４） 作出判断：若Ｙ<或Ｙ>，则拒绝，否则不拒绝。

附：Ｂox-Cox变换的Ｒ代码

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 
} 
进行变换的Ｒ代码 
> attach(Prestige) //car package 
> hist(income) 
> BoxCox_Trans(income,seq(-3,3)) 
$minimum 
[1] 0.1792894 
$Objective 
[1] 827.9459 
> hist(income^0.1792894)