函数与优化

练习：
编写一个函数计算向量的最大5个数的均值，并返回最大的5个值
思路：
1.对一个向量排序，用sort()函数，此函数默认从小到大排序；再用函数rev()转成由大到小的序列
2.提取前5个值，求均值
3.返回函数值

vms <- function(x){  x1 = rev(sort(x))  x2 = sum(x1[1:5])/5  return(list(xbar=x2,top5=x1[1:5]))}y <- c(5,13,4,20,30,23,27,28,47,35,37)vms(y)$xbar[1] 35.4$top5[1] 47 37 35 30 28

程序运行时间和效率

• R中的proc.time()函数可以返回当前R已经运行的时间(单位是s)

  user system elapsed
  4.10 0.79 1592.78
  其中，user是指R执行用户指令的CPU运行时间，system是系统所需要的时间,elapsed是指R打开到现在总共运行的时间。

• 计算程序运行时间的函数是:system.time(expr,gcFirst=TRUE)

system.time(for(i in 1:100) mad(runif(1000))) ser  system elapsed 0.02    0.00    0.04 该程序等价于：ptm <- proc.time()for(i in 1:100) mad(runif(1000))proc.time()-ptmuser  system elapsed 0.02    0.00    0.01 

topics:R中若知道向量Z的长度，可以首先给出Z的长度，这样可以提高运行效率

---

缺失值与空值的比较

z <- rep(NA,4)for(i in 1:10){  z <- c(z,i+1)} z [1] NA NA NA NA  2  3  4  5  6  7  8  9 10 11NA运算的结果还是NA；z <- c()for(i in 1:10){  z <- c(z,i+1)} z [1]  2  3  4  5  6  7  8  9 10 11

---

1. 优化求解
   一元函数优化求解

optimize(f,interval,...,lower=min(interval),upper=max(interval),maximum=FALSE,tol=.Machine$double.eps^0.25)

其中：f是需要求优化的函数;interval是参数搜索的区间;lower是参数搜索的下限,如果缺失，用interval的最小值代替；
upper是参数搜索的上限,缺失时用interval的最大值代替;maximum 默认为求极小值;tol是精度的容忍度。
练习:
求解f(x)=ln(x)-x^2此函数的极大值
思路：
1.先编写函数f;
2.作f函数在区间(0,2)上的图
3.求f函数的最大值的返回结果

> f <- function(x) log(x)-x^2> curve(f,xlim=c(0,2))> optimize(f,c(0.1,10),tol=0.0001,maximum=T)$maximum[1] 0.7071232$objective[1] -0.8465736

结果显示：当x=0.707时，目标函数的最大值为-0.8465736

• 多元函数优化求解
  用optim()求解

  optim(par,fn,gr=NULL,...,method=c("Nelder-Mead","BFGS&","CG&","L-BFGS-B&","SANN"),      lower=-INF,upper=INF,control=list(),hessian=FALSE)

其中：par是设定初始值,fn是需要优化的目标函数，gr是梯度向量,如果为空，则由optim()计算所得的近似值替代。
control是用来控制optim函数的一些参数,hessian=FALSE表示不需要返回海塞矩阵；
method()是一些迭代算法;
对多元函数：1.写出目标函数；2.写出梯度表达式；3.画出三维视图，从直观上判断其极值(最大值或最小值)；4.作优化计算

##目标函数f <- function(x){  x1=x[1]  x2=x[2]  (x1^2+x2-11)^2+(x1+x2^2-7)^2}#梯度向量grr <- function(x){  x1=x[1]  x2=x[2]  c(2*(x1^2+x2-11)*2*x1+2*(x1+x2^2-7),2*(x1^2+x2-11)+2*(x1+x2^2-7)*2*x2)}#画图判断极值点x1<-x2<-seq(-10,10,length=100)f <- function(x1,x2){  (x1^2+x2-11)^2+(x1+x2^2-7)^2}z<- outer(x1,x2,f)  #求外积contour(x1,x2,z)   #画等高线persp(x1,x2,z) #画三维图

由上图可以看出函数存在极小值点

#求最优化，设定初始值x1=-5,x2=5optim(c(-5,5),fn=f,gr=grr)$par[1] -2.804927  3.131364  #最终的参数值$value    #极小目标值[1] 1.306985e-06$countsfunction gradient       55       NA $convergence[1] 0$messageNULL

• 约束条件下的优化求解
  有两种求法：
  一种是利用constrOptim()函数直接求解：

constrOptim(theta,f,grad,ui,ci,mu=1e-04,control=list(),method=if(is.null(grad))"Nelder-Mead" else "BFGS",outer.iterations = 100,outer.eps = 1e-05,...,hessian=FALSE)   ```         其中，theta是初始值向量;f是目标函数;grad是梯度向量，可以是空值NULL；ui是约束矩阵的左边系数矩阵;ci是约束矩阵的右边的值例如：对一个二元函数，其约束条件是x1>0,x2>0，等价于1*x1+0*x2>00*x1+1*x2>0,  ui即为x1,x2的系数矩阵,ci=t(c(0,0)).

#编写目标函数fr1 <- function(x){  x1=x[1]  x2=x[2]  (x1^2+x2-11)^2+(x1+x2^2-7)^2}#一阶梯度表达式grr <- function(x){  x1=x[1]  x2=x[2]  c(2*(x1^2+x2-11)*2*x1+2*(x1+x2^2-7),2*(x1^2+x2-11)+2*(x1+x2^2-7)*2*x2)}uimat <- matrix(c(1,0,0,1),nr=2) #约束矩阵的左边系数矩阵cimat <- matrix(c(0,0),nr=2)  #约束矩阵的右边的值cop <- constrOptim(c(0.2,0.5),fr1,grr,ui=uimat,ci=cimat)$par  #目标函数达到极值时,x1和x2的取值分别为（3,2）[1] 3 2$value  #目标函数的极值[1] 1.066434e-22#表示迭代过程中用到目标函数function52次，梯度函数gradient16次$counts  function gradient       52       16 $convergence   #convergence取0表示成功收敛[1] 0   $messageNULL$outer.iterations[1] 3$barrier.value[1] 3.178688e-05

另外一种方法是通过参数转换为无约束下的最优化问题
例如：原来的约束条件是x1>0,x2>0，可以用指数变换x1=exp(z1),x2=exp(z2),z1,z2没有任何限制

#编写目标函数fr2 <- function(z){  x1=exp(z[1])  x2=exp(z[2])  (x1^2+x2-11)^2+(x1+x2^2-7)^2}#一阶梯度表达式grr <- function(z){  x1 <- exp(z[1])  x2 <- exp(z[2])  c(2*(x1^2+x2-11)*2*x1*exp(z[1])+2*(x1+x2^2-7)*exp(z[1]),    2*(x1^2+x2-11)*exp(z[2])+2*(x1+x2^2-7)*2*x2*exp(z[2]))}#求解最优化值，此处最大值返回的是z值的取值optran <- optim(c(-1.6,-0.7),fr2,grr)optran$par[1] 1.0986436 0.6930667$value[1] 4.647052e-07$countsfunction gradient      115       NA $convergence[1] 0$messageNULL#将z值转换为x的取值exp(optran$par)[1] 3.000094 1.999839