蒙特卡罗最优化

假设f:R→R为一连续函数，则方程f(x)=c的 根x,满足g(x)=f(x)-c=0.为此我们只考虑 f(x)=0形式的方程求根问题。使用数值方法 求此方程的根，可以选择是使用f的一阶导 数还是不使用导数的方法。Newton方法或 者Newton-Raphson方法是使用一阶导数的 方法，而Brent的最小化算法是不使用导数 的一种求根方法。
2.2EM算法
EM算法是个聚类算法，即根据给定观察 数据自动对数据进行分类。 问题来源
K
•给定一些观察数据x，假设x符合如下高斯分布：
p ( x) k N ( x k , )
k 1 2 k
•求混合高斯分布的三组参数
k k k2
简单问题
该混合高斯分布一共有K个分布，并且对于每个观察 到的x，如果我们同时还知道它属于K中的哪一个分布， 则我们可以根据最大似然估计求出每个参数。 结论：
Monte Carlo Optimization
主要内容
一、数值优化方法（Numerical optimization methods） 二、应用于求解随机优化问题的蒙特卡罗方 法 （1）模拟退火算法（Simulated Annealing) （2）EM算法（The EM algorithm)
1.Numerical optimization methods in R 1.1 Root-finding in one dimension
EM算法过程：
1、用随机函数初始化K个高斯分布的参数，同时 保证

k 1
K
k
1
Expectation 2、依次取观察数据x，比较x在K个高斯函数 中概率的大小，把x归类到这K个高斯中 概率最大的一个。 ( x )( x )T
k
1

k
k
2 k
e
2 k 2
Maximum 3、 用最大似然估计，使观察数据是x的概率 最大，因为已经在第2步中分好类了，所 以，即简单问题的求法。
1 k Nk
x
xk
k
1 Nk
2 k
( xk k )( xk k )T
xk
Nk k N
ห้องสมุดไป่ตู้4、返回第2步用第3步新得到的参数来对观察数据x 重新分类。直到下式概率（最大似然函数）达 到最大。 ( x )( x )T


x
1
k
k
k
2 k
e
2 k 2
• • • • • • • • • • •
b1<-b0 b0<-b0-1 eps <- .Machine$double.eps^0.25 it<-0 while(it<1000 && abs(b1-b0)>eps){ it<-it+1 b0<-b1 b1<-b0-f(b0,a,n)/fd(b0,a,n) cat(it,c(b0,b1,abs(b1-b0)),"\n") } }
运行结果：
运行结果：
2.应用于求解随机优化问题的蒙特卡罗方法
2.1模拟退火算法
模拟退火算法来源于固体退火原理，将 固体加温至充分高，再让其徐徐冷却，加 温时，固体内部粒子随温升变为无序状， 内能增大，而徐徐冷却时粒子渐趋有序， 在每个温度都达到平衡态，最后在常温时 达到基态，内能减为最小。根据Metropolis 准则，粒子在温度T时趋于平衡的概率为eΔE/(kT)，其中E为温度T时的内能，ΔE为 其改变量，k为Boltzmann常数。用固体退 火模拟组合优化问题，将内能E模拟为目标 函数值f，温度T演化成控制参数t， 即得到
1.1.3 Newton’s method
例3 使用Newton方法求例1方程的根。 程序： • nt<-function(b0){ • a <- 0.5 • n <- 20 • f <- function(y, a, n) { • a^2 + y^2 + 2*a*y/(n-1) - (n-2) • } • fd<-function(y,a,n){ • 2*y+2*a/(n-1) • }
• 运行结果： • true roots -4.239473 4.186841
1.1.2 Brent’s method
二分法是一种特殊的括入根算法。Brent通 过逆二次插值方法将括入根方法和二分法 结合起来。其使用y的二次函数来拟合x。 如果三个点为(a,f(a)),(b,f(b)),(c,f(c)）,其中b 为当前最好的估计，则通过Lagrange多项 式插值方法（y=0)对方程的根进行估计，
while(it < 1000 && abs(y[2]) > 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) print(c(r[1], y[1], y[3]-y[2])) } } bisec(0,5*n)
解组合优化问题的模拟退火算法：由初始 解i和控制参数初值t开始，对当前解重复 “产生新解→计算目标函数差→接受或舍 弃”的迭代，并逐步衰减t值，算法终止时 的当前解即为所得近似最优解，这是基于 蒙特卡罗迭代求解法的一种启发式随机搜 索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其 衰减因子Δt、每个t值时的迭代次数L和停止 条件S
下面我们使用二分法求此方程的一个数值解。我 们首先要找到一个区间，比如(0,5n)，使得函数 在区间两端有着不同的符号。然后即可使用二分 法。
程序： 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){ f <- function(y, a, n) { a^2 + y^2 + 2*a*y/(n-1) - (n-2) } it <- 0 eps <- .Machine$double.eps^0.25 r <- seq(b0, b1, length=3) y <- c(f(r[1], a, n), f(r[2], a, n), f(r[3], a, n)) if (y[1] * y[3] > 0) stop("f does not have opposite sign at endpoints")
1.1.1 Bisection method(二分法）
如果f(x)在区间[a,b]上连续，以及f(a)和f(b) 有相反的符号，则由中值定理知道存在 a<c<b，使得f(c)=0。二分法通过在每次迭 代中简单的判断f(x)在中点x=(a+b)/2处的符 号来寻求方程的根。如果f(a)和f(x)有相反的 符号则区间就被[a,x]代替，否则就被[x,b]代 替。在每次迭代中，包含根的区间长度减 少一半。即
xk
1 k Nk
x
xk
k
表示属于第k个高斯 分布的观察数据x。
1 Nk
2 k
( xk k )( xk k )T
xk
Nk k N
特别注意 k 是个 向量，而 k 2 是个 数值。
实际问题 观察数据x属于哪个高斯分布是未知的 所以要用EM算法来解决这种实际问题。
常用的收敛准则有：
绝对收敛
可以看出，二分法不会失效，达到指定精 度所需要的迭代次数也是事先可以得到的。 如果在区间[a,b]里方程有多个根，则二分
法会找到一个根。二分法的收敛速度是线 性的。 相对收敛
时停止迭代。此准则可以不考虑x的单位情况下达 到指定的精度。
例1 解方程
其中a为常数，n>2为一整数。显然，方程的 解为
输入：nt(5) 输出结果： 1 5 4.252618 0.7473822 2 4.252618 4.187347 0.06527095 3 4.187347 4.186841 0.0005055338 4 4.186841 4.186841 3.032932e-08
Newton方法依赖于f的形状和初值。该方法 从初值开始就发散。
在R中，函数uniroot就是应用Brent方法 求解一元方程的数值根。
例2 应用uniroot求例1中的方程的根。 程序： 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) root f.root iter estim.prec 4.186870e+00 2.381408e-04 1.400000e+01 6.103516e-05 uniroot(function(y) {a^2 + y^2 + 2*a*y/(n-1) - (n-2)}, interval = c(-n*5, 0))$root [1] -4.239501