蒙特卡罗最优化课件

MonteCarloOptimization主要内容一、数值优化方法〔Numericaloptimizationmethods〕二、应用于求解随机优化问题的蒙特卡罗方法〔1〕模拟退火算法〔SimulatedAnnealing)〔2〕EM算法〔TheEMalgorithm)1.NumericaloptimizationmethodsinR1.1Root-findinginonedimension假设f:R→R为一连续函数，那么方程f(x)=c的根x,满足g(x)=f(x)-c=0.为此我们只考虑f(x)=0形式的方程求根问题。使用数值方法求此方程的根，可以选择是使用f的一阶导数还是不使用导数的方法。Newton方法或者Newton-Raphson方法是使用一阶导数的方法，而Brent的最小化算法是不使用导数的一种求根方法。1.1.1Bisectionmethod(二分法〕如果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]代替。在每次迭代中，包含根的区间长度减少一半。即下面我们使用二分法求此方程的一个数值解。我们首先要找到一个区间，比方(0,5n)，使得函数在区间两端有着不同的符号。然后即可使用二分法。例1解方程其中a为常数，n>2为一整数。显然，方程的解为程序：a<-0.5n<-20cat("trueroots",-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<-0eps<-.Machine$double.eps^0.25r<-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("fdoesnothaveoppositesignatendpoints")while(it<1000&&abs(y[2])>eps){it<-it+1if(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])/2y[2]<-f(r[2],a=a,n=n)print(c(r[1],y[1],y[3]-y[2]))}}bisec(0,5*n)运行结果：trueroots-4.2394734.1868411.1.2Brent’smethod二分法是一种特殊的括入根算法。Brent通过逆二次插值方法将括入根方法和二分法结合起来。其使用y的二次函数来拟合x。如果三个点为(a,f(a)),(b,f(b)),(c,f(c)〕,其中b为当前最好的估计，那么通过Lagrange多项式插值方法〔y=0)对方程的根进行估计，在R中，函数uniroot就是应用Brent方法求解一元方程的数值根。例2应用uniroot求例1中的方程的根。程序：a<-0.5n<-20out<-uniroot(function(y){a^2+y^2+2*a*y/(n-1)-(n-2)},lower=0,upper=n*5)unlist(out)rootf.rootiterestim.prec4.186870e+002.381408e-041.400000e+016.103516e-05uniroot(function(y){a^2+y^2+2*a*y/(n-1)-(n-2)},interval=c(-n*5,0))$root[1]-4.239501

1.1.3Newton’smethod例3使用Newton方法求例1方程的根。程序：nt<-function(b0){a<-0.5n<-20f<-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)}b1<-b0b0<-b0-1eps<-.Machine$double.eps^0.25it<-0while(it<1000&&abs(b1-b0)>eps){it<-it+1b0<-b1b1<-b0-f(b0,a,n)/fd(b0,a,n)cat(it,c(b0,b1,abs(b1-b0)),"\n")}}输入：nt(5)输出结果：

1

54.2526180.747382224.2526184.1873470.0652709534.1873474.1868410.000505533844.1868414.1868413.032932e-08

Newton方法依赖于f的形状和初值。该方法从初值开始就发散。运行结果：运行结果：

2.应用于求解随机优化问题的蒙特卡罗方法2.1模拟退火算法模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都到达平衡态，最后在常温时到达基态，内能减为最小。根据Metropolis准那么，粒子在温度T时趋于平衡的概率为e-ΔE/(kT)，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃〞的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(CoolingSchedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S给定一些观察数据x，假设x符合如下高斯分布：求混合高斯分布的三组参数2.2EM算法问题来源EM算法是个聚类算法，即根据给定观察数据自动对数据进行分类。该混合高斯分布一共有K个分布，并且对于每个观察到的x，如果我们同时还知道它属于K中的哪一个分布，那么我们可以根据最大似然估计求出每个参数。结论：简单问题特别注意是个向量，而是个数值。表示属于第k个高斯分布的观察数据x。实际问题观察数据x属于哪个高斯分布是未知的所以要用EM算法来解决这种实际问题。EM算法过程：1、用随机函数初始化K个高斯分布的参数，同时保证Expectation2、依次取观察数据x，比较x在K个高斯函数中概率的大小，把x归类到这K个高斯中概率最大的一个。