R语言包_stats::optim

stats中的optim函数是解决优化问题的一个简易的方法。

Univariate Optimization

f = function(x,a) (x-a)^2
xmin = optimize(f,interval = c(0,1),a=1/3)
xmin

General Optimization

optim函数包含了几种不同的算法。 算法的选择依赖于求解导数的难易程度,通常最好提供原函数的导数。

在求解之前,一般需要scale。 可以尝试用不同的方法求解同样的问题。

Nelder-Mead method

optim默认的方法。

又称下山单纯形法,可做非线性函数的极值以及曲线拟合。

其主要思想是: 在n维空间构建(n+1)顶点的多面体,通过reflection,expansion,contraction,来逐步逼近最佳点x∗x^*。

特点是: 1. 不适用函数的导数信息 2. 对不可导函数适用 3. 可能很慢

BFGS method

属于quasi-Newton方法。

首先,简单介绍牛顿法: 牛顿法基于目标函数的二阶导数(海森矩阵),收敛速度快,迭代次数少,尤其在最优值附近,收敛速度是二次的。 缺点是:海森矩阵稠密时,每次迭代计算量交大,且每次都会重新计算目标函数的海森矩阵的逆。这样以来,问题规模大时,其计算量以及存储空间都很大。

拟牛顿法是在牛顿法基础上的改进,其引入了海森矩阵的近似矩阵,避免了每次迭代都需要计算海森矩阵的逆,其收敛速度介于梯度下降和牛顿法之间,属于超线性。

同时,牛顿法在每次迭代时不能保证海森矩阵总是正定的,一旦其不是正定,优化方向就会跑偏,从而使牛顿法失效,也证明了牛顿法的鲁棒性较差。 拟牛顿法利用海森矩阵的逆矩阵代替海森矩阵,虽然每次迭代不一定保证最优化的方向,但是近似矩阵始终正定,因此算法总是朝着最优值搜索。

注意: 1. 使用函数导数信息,通过人工提供或者有限微分 2. 高维的数据存储会很大

CG method

一种共轭梯度法(conjugate gradient),选择连续的、与椭圆轴线相仿的路径。

特点: 1. 不存储海森矩阵 2. 三种不同的路径搜索方法 3. 与BFGS相比,较差的鲁棒性 4. 使用函数导数信息

L-BFGS-B method

A limited memory version of BFGS

特点: 1. 不存储海森矩阵,只有一个对海森矩阵大小受限的更新步骤。 2. 使用导数信息 3. 可以把解决方法限制到box里,是optim中仅有的方法。

SANN method

模拟退火法(simulated annealing)的变种。

特点: 1. 随机算法 2. 接受以正概率提升目标的改变 3. 不使用导数信息 4. 收敛很慢,但是找到一个good solution很快

Brent method

An interface to optimize

特点: 1. 仅适用于一维问题 2. 可以在其他函数中包含

How to use

optim

control options

components of returned value

constrained optimization

Exampels

One Dimensional Ex1

假定

f(x)=e−(x−2)2

f(x)=e^{-(x-2)^2} 其导数为

f′(x)=−2(x−2)f(x)

f'(x)=-2(x-2)f(x)

# we supply negative f, since we want to maximize.
f <- function(x) -exp(-( (x-2)^2 ))
######### without derivative
# I am using 1 at the initial value
# $par extracts only the argmax and nothing else
optim(1, f)$par
######### with derivative
df <- function(x) -2*(x-2)*f(x)
optim(1, f, df, method="CG")$par
######### with "Brent" method
optim(1,f,method="Brent",lower=-0,upper=3)$par

# figure
x = seq(0,3,length.out = 100)
y = f(x)
plot(x,y)

One Dimensional Ex2

假定

f(x)=sin(xcos(x))

f(x)=sin(xcos(x))

#算法可以很快地发现与初值相近的局部最优值。
f <- function(x) sin(x*cos(x))
optim(2, f)$par
optim(4, f)$par
optim(6, f)$par
optim(8, f)$par

Two Dimensional Ex3

Rosenbrock function:

f(x,y)=(1−x)2+100(y−x2)2

f(x,y)=(1-x)^2+100(y-x^2)^2

This function is strictly positive, but is 0 when y = x^2, and x = 1, so (1, 1) is a minimum. Let’s see if optim can figure this out. When using optim for multidimensional optimization, the input in your function definition must be a single vector.

# 绘图
f <- function(x1,y1) (1-x1)^2 + 100*(y1 - x1^2)^2
x <- seq(-2,2,by=.15)
y <- seq(-1,3,by=.15)
z <- outer(x,y,f)
persp(x,y,z,phi=45,theta=-45,col="yellow",shade=.00000001,ticktype="detailed")
# 求解
f <- function(x) (1-x[1])^2 + 100*(x[2]-x[1]^2)^2
# starting values must be a vector now
optim( c(0,0), f )$par
[1] 0.9999564 0.9999085

Two Dimensional Ex4

Himmelblau’s function:

f(x,y)=(x2+y−11)2+(x+y2−7)2

f(x,y)=(x^2+y-11)^2+(x+y^2-7)^2

There appear to be four “bumps” that look like minimums in the realm of (-4,-4), (2,-2),(2,2) and (-4,4). Again this function is strictly positive so the function is minimized when x^2 + y − 11 = 0 and x + y^2 − 7 = 0.

#画图
f <- function(x1,y1) (x1^2 + y1 - 11)^2 + (x1 + y1^2 - 7)^2
x <- seq(-4.5,4.5,by=.2)
y <- seq(-4.5,4.5,by=.2)
z <- outer(x,y,f)
persp(x,y,z,phi=-45,theta=45,col="yellow",shade=.65 ,ticktype="detailed")
#求解局部最优值
f <- function(x) (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
optim(c(-4,-4),f)$par
optim(c(2,-2), f)$par
optim(c(2,2), f)$par
optim(c(-4,4),f)$par
#which are indeed the true minimums. This can be checked by seeing that these inputs
correspond to function values that are about 0.

Fit a probit regression model

pass

Minimise residual sum of squares

# 初始化数据
d = data_frame(x=1:6,y=c(1,3,5,6,8,12))
d
ggplot(d,aes(x,y)) + geom_point() +  stat_smooth(method = "lm")
# 最小问题的优化函数
min.RSS = function(data,par) {
    with(data,sum((par[1]+par[2]*x-y)^2))
}
#optim函数调用的格式
result = optim(par=c(0,0),min.RSS,data=d)
#optim调用的结果参数
result$par
result$value
result$counts
result$convergence
result$message
#可视化分析结果
ggplot(d,aes(x,y)) + geom_point() +geom_abline(intercept=result$par[1],
    slope=result$par[2],color="red")
#optim与lm的结果对比分析
lm(y~x,data=d)
result$par

Maximum likelihood

To fit a Poisson distribution to x I don’t minimise the residual sum of squares, instead I maximise the likelihood for the chosen parameter lambda.

# 观测数据
obs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 17, 42, 43)
freq = c(1392, 1711, 914, 468, 306, 192, 96, 56, 35, 17, 15, 6, 2, 2, 1, 1)
x <- rep(obs, freq)
plot(table(x), main="Count data")
qplot(x,stat_bin=1)
# 优化函数,注意“-”号
lklh.poisson <- function(x, lambda) lambda^x/factorial(x) * exp(-lambda)
log.lklh.poisson <- function(x, lambda){ 
    -sum(x * log(lambda) - log(factorial(x)) - lambda) 
}
# 调用optim
optim(par=2,log.lklh.poisson,x=x)
optim(par=2,log.lklh.poisson,x=x,method="Brent",lower=2,upper=3)
# 比较结果
library(MASS)
fitdistr(x, "Poisson")
mean(x)
# 系统信息
sessionInfo()

参考资料

  1. magesBlog
  2. ime.unicamp

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券