函数与优化

练习：
编写一个函数计算向量的最大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()-ptm
user  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 11
NA运算的结果还是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=5
optim(c(-5,5),fn=f,gr=grr)

$par
[1] -2.804927  3.131364  #最终的参数值

$value    #极小目标值
[1] 1.306985e-06

$counts
function gradient 
      55       NA 

$convergence
[1] 0

$message
NULL

• 约束条件下的优化求解
  有两种求法：
  一种是利用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>0
0*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   

$message
NULL

$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

$counts
function gradient 
     115       NA 

$convergence
[1] 0

$message
NULL

#将z值转换为x的取值
exp(optran$par)
[1] 3.000094 1.999839