# 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函数包含了几种不同的算法。 算法的选择依赖于求解导数的难易程度，通常最好提供原函数的导数。

optim默认的方法。

## L-BFGS-B method

A limited memory version of BFGS

## Brent method

An interface to optimize

# 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)
# 求解
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)
#求解局部最优值
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")
# 最小问题的优化函数
with(data,sum((par[1]+par[2]*x-y)^2))
}
#optim函数调用的格式
#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()```

145 篇文章46 人订阅

0 条评论