中国科学技术大学：《实用统计软件》课程课件讲义（统计计算与软件）第十一讲 R中的数值优化方法

Lecture11:数值优化方法 张伟平 Monday 23d November,2009

Lecture 11: Íä`zê{ ‹ï² Monday 23rd November, 2009

Contents 1 Numerical optimization methods in R y 1.1 Root-finding in one dimension............ 1.1.1 Bisection method............... 2 1.1.2 Brent's method .... 6 1.1.3 Newton's method..················ 11.4 Fisher scoring..·······.·······.·17 l.2 multivariate optimization··········.········ 18 1.2.1 Newton's method and Fisher scoring.........19 l.3 Numerical Integration.·........·.···.·.···20 1.4 Maximum Likelihood Problems . 27 1.5 Optimization Problems........ 29 1.5.1One-dimension Optimization...。。.....·.· 29 1.5.2 multi-dimensional Optimization.....·.·..·· 34 1.6 Linear Programming..·.···..·.·····.·..··38 Previous Next First Last Back Forward 1

Contents 1 Numerical optimization methods in R 1 1.1 Root-finding in one dimension . . . . . . . . . . . . . . . . . 1 1.1.1 Bisection method . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Brent’s method . . . . . . . . . . . . . . . . . . . . . 6 1.1.3 Newton’s method . . . . . . . . . . . . . . . . . . . . 8 1.1.4 Fisher scoring . . . . . . . . . . . . . . . . . . . . . . 17 1.2 multivariate optimization . . . . . . . . . . . . . . . . . . . 18 1.2.1 Newton’s method and Fisher scoring . . . . . . . . . 19 1.3 Numerical Integration . . . . . . . . . . . . . . . . . . . . . 20 1.4 Maximum Likelihood Problems . . . . . . . . . . . . . . . . 27 1.5 Optimization Problems . . . . . . . . . . . . . . . . . . . . . 29 1.5.1 One-dimension Optimization . . . . . . . . . . . . . 29 1.5.2 multi-dimensional Optimization . . . . . . . . . . . . 34 1.6 Linear Programming . . . . . . . . . . . . . . . . . . . . . . 38 Previous Next First Last Back Forward 1

Chapter 1 Numerical optimization methods in R 1.1 Root-finding in one dimension 假设f:R→R为一连续函数，则方程f(x)=c的根x,满足g(x)=f(x)一c= 0.因此我们只考虑f(x)=0形式的方程求根问题.使用数值方法求此 方程的根，可以选择是使用f的一阶导数还是不使用导数的方法.Nw- ton方法或者Newton-Raphson方法是使用一阶导数的方法，而Brent的最 小化算法l是不使用导数的一种求根方法.在R中，函数uniroot就是基 于Brent的求根算法.该算法的Fortran程序源代码可以在下面网址上找 到http://www.gnu.org/software/gsl/月 IR.Brent.Algorithms for minimization without derivatives.Prentice-Hall,New Jersey,1973 Previous Next First Last Back Forward 1

Chapter 1 Numerical optimization methods in R 1.1 Root-finding in one dimension bf : R → RèòÎYºÍ, Kêßf(x) = cäx, ˜vg(x) = f(x) − c = 0. œd·Ç êƒf(x) = 0/™ê߶äØK. ¶^Íäê{¶d êßä, 屿J¥¶^fòÍÑ¥ ÿ¶^Íê{. New￾tonê{½ˆNewton-Raphsonê{¥¶^òÍê{, BrentÅ z é{1¥ÿ¶^Íò´¶äê{. 3R•, ºÍuniroot“¥ƒ uBrent¶äé{. Té{Fortran ßS ìËå±3e°å˛È http://www.gnu.org/software/gsl/. 1R. Brent. Algorithms for minimization without derivatives. Prentice-Hall, New Jersey, 1973 Previous Next First Last Back Forward 1

1.1.1 Bisection method 如果f(x)在区间[a,b上连续，以及f(a)和f(b)有相反的符号，则由中值定理知 道存在a0 以及x+1)=(at+1+b+1)/2 3.狮环2直至达到收敛准则。 常用的收敛准则有 Previous Next First Last Back Forward 2

1.1.1 Bisection method XJf(x)3´m[a, b]˛ÎY, ±9f(a)⁄f(b)kÉጓ, Kd•ä½n  3a 0 ±9 x (t+1) = (at+1 + bt+1)/2. 3. ÃÇ2Üñà¬ÒOK. ~^¬ÒOKk Previous Next First Last Back Forward 2

绝对收敛 |x(t+1)-x(1lo92{(b0-a0)/6}-1,将达到容忍 精度x()-x<6. 可以看出，二分法不会失效，达到指定精度所需要的迭代次数也是事先可以 得到的.如果在区间[α，里方程有多个根，则二分法会找到一个根.二分法的 收敛速度是线性的， 相对收敛 相对收敛准则要求 z(t+1)-z(t) <E 2() Previous Next First Last Back Forward 3

˝È¬Ò |x (t+1) − x (t) | log2{(b0 − a0)/δ} − 1, ÚàN= °›|x (t) − x ∗| < δ. å±w—, ©{ÿ¨î, àç½°›§IáSìgÍè¥Økå± . XJ3´m[a, b]p êßkıáä, K©{¨Èòáä. ©{ ¬ÒÑ›¥Ç5. ÉÈ¬Ò ÉȬÒOKᶠ|x (t+1) − x (t) | |x(t) | <  Previous Next First Last Back Forward 3

时停止迭代.此准则可以不考虑x的单位的情况下达到指定的精度. 例1解方程 a2+y2+2ay =n-2 n-1 其中a为常数，n>2为一整数.显然，方程的解为 y=-n2±Vn-2+a2+n 下面我们使用二分法求此方程的一个数值解.我们首先要找一个区间，比 如(0,5n),使得函数f()=a2+y2+兴-n+2在区间两端有着不同的符 号.然后即可以使用二分法 FCode a<-0.5 n<-20 cat("true roots",-a/(n-1)-sqrt(n-2-a*2+(a/(n-1))2), -a/(n-1)+sqrt(n-2-a^2+(a/(n-1))2),"\n") bisec<-function(b0,b1){ Previous Next First Last Back Forward 4

û éSì. dOKå±ÿƒx¸†ú¹eàç½°›. ~1 )êß a 2 + y 2 + 2ay n − 1 = n − 2 Ÿ•aè~Í, n > 2èòÍ. w,, êß)è y = − a n − 1 ± r n − 2 + a 2 + ( a n − 1 ) 2 e°·Ç¶^©{¶dêßòáÍä). ·ÇƒkáÈ òá´m, ' X(0, 5n), ¶ºÍf(y) = a 2 + y 2 + 2ay n−1 − n + 23´m¸‡kXÿ”Œ “. ,￾=屶^©{. ↑Code a <- 0.5 n <- 20 cat("true roots",-a/(n-1)-sqrt(n-2-a^2+(a/(n-1))^2), -a/(n-1)+sqrt(n-2-a^2+(a/(n-1))^2),"\n") bisec<-function(b0,b1){ Previous Next First Last Back Forward 4

f 0) stop("f does not have opposite sign at endpoints") while(it 1000&&abs(y[2])>eps){ 1t<-1t+1 if(y[1]*y[2]<0)[ r[3]<-r[2] y[3]<-y[2] else r[1]<-r[2] y[1]<-y[2] } r[2]<-(r[1]+r[3])/2 y[2]<-f(r[2],a=a,nn) Previous Next First Last Back Forward 5

f 0) stop("f does not have opposite sign at endpoints") while(it eps) { it <- it + 1 if (y[1]*y[2] < 0) { r[3] <- r[2] y[3] <- y[2] } else { r[1] <- r[2] y[1] <- y[2] } r[2] <- (r[1] + r[3]) / 2 y[2] <- f(r[2], a=a, n=n) Previous Next First Last Back Forward 5

pr1nt(c(r[1],y[1],y[3]-y[2]) bisec(0,5*n) Code 精确的解为y=4.186841,-4.239473. 1.1.2 Brent's method 二分法是一种特殊的括入根算法.Bent通过逆二次插值方法将括入 根方法和二分法结合起来.其使用y的二次函数来拟合x.如果三个点 为(a,f(a),(b,f(b),(c,f(c),其中b为当前最好的估计，则通过Lagrange多 项式插值方法(y=0)对方程的根进行估计， 【y-f(alg-f(⑥c ly-f(b)lly-f(c)la x= [f(e)-f(a)llf(e)-f(b)][f(a)-f(6)lIf(a)-f(c)] Previous Next First Last Back Forward 6

print(c(r[1], y[1], y[3]-y[2])) } } bisec(0,5*n) ↓Code °()èy = 4.186841, −4.239473. 1.1.2 Brent’s method ©{¥ò´Aœ)\äé{. Brent œL_gäê{Ú)\ äê{⁄©{(‹Â5. Ÿ¶^ygºÍ5[‹x. XJná: è(a, f(a)), (b, f(b)), (c, f(c)), Ÿ•b ècÅ–O, KœLLagrangeı ë™äê{(y = 0)Èêßä?1O, x = [y − f(a)][y − f(b)]c [f(c) − f(a)][f(c) − f(b)] + [y − f(b)][y − f(c)]a [f(a) − f(b)][f(a) − f(c)] Previous Next First Last Back Forward 6

ly-f(c)]ly-f(a) +[f(b)-f(e)lf(c)-f(a)l 在R中，函数uniroot.就是应用Brent方法求解一元方程的数值根. 例2应用uniroot求例1中方程的根. FCode a<-0.5 n<-20 out <uniroot(function(y){ a2+y^2+2*a*y/(n-1)-(n-2)], lower =0,upper n*5) unlist(out) uniroot(function(y){a2 y2 2*a*y/(n-1)-(n-2)}, interval c(-n*5,0))$root Code 函数polyroot.求一个系数为复数的多项式的根.因此，此例的问题使 用polyroot函数为 Previous Next First Last Back Forward 7

+ [y − f(c)][y − f(a)]b [f(b) − f(c)][f(c) − f(a)] 3R•, ºÍuniroot“¥A^Brentê{¶)òêßÍää. ~2 A^uniroot¶~1•êßä. ↑Code a <- 0.5 n <- 20 out <- uniroot(function(y) { a^2 + y^2 + 2*a*y/(n-1) - (n-2) }, lower = 0, upper = n*5) unlist(out) uniroot(function(y) {a^2 + y^2 + 2*a*y/(n-1) - (n-2)}, interval = c(-n*5, 0))$root ↓Code ºÍpolyroot¶òáXÍèEÍıë™ä. œd, d~ØK¶ ^polyrootºÍè Previous Next First Last Back Forward 7

TCode po1 yroot(c(a2-(n-2),2*a/(n-1),1)) Code 1.1.3 Newton's method Newton方法是一种快速求根的方法，有时也称为Newton-Raphson迭代算法. 假设f(x)是连续且可微的，且f'(x)≠0.假设我们想要求根x,我们已经有了 当前的一个近似值x,则由下图可以看出，通过导数表示斜率这一特点，可以 给出x的一个更好的近似 Previous Next First Last Back Forward 8

↑Code polyroot(c(a^2-(n-2),2*a/(n-1),1)) ↓Code 1.1.3 Newton’s method Newtonê{¥ò´ØѶäê{, kûè°èNewton-RaphsonSìé{. bf(x)¥ÎYÖåá, Öf 0 (x) 6= 0. b·Çéá¶äx, ·ÇƲk còáCqäxn, Kde„å±w—, œLÍL´«˘òA:, å± â—xòáç–Cq. Previous Next First Last Back Forward 8