蒙特卡羅最優(yōu)化課件

MonteCarloOptimization主要內(nèi)容一、數(shù)值優(yōu)化方法〔Numericaloptimizationmethods〕二、應(yīng)用于求解隨機(jī)優(yōu)化問(wèn)題的蒙特卡羅方法〔1〕模擬退火算法〔SimulatedAnnealing)〔2〕EM算法〔TheEMalgorithm)1.NumericaloptimizationmethodsinR1.1Root-findinginonedimension假設(shè)f:R→R為一連續(xù)函數(shù)，那么方程f(x)=c的根x,滿足g(x)=f(x)-c=0.為此我們只考慮f(x)=0形式的方程求根問(wèn)題。使用數(shù)值方法求此方程的根，可以選擇是使用f的一階導(dǎo)數(shù)還是不使用導(dǎo)數(shù)的方法。Newton方法或者Newton-Raphson方法是使用一階導(dǎo)數(shù)的方法，而B(niǎo)rent的最小化算法是不使用導(dǎo)數(shù)的一種求根方法。1.1.1Bisectionmethod(二分法〕如果f(x)在區(qū)間[a,b]上連續(xù)，以及f(a)和f(b)有相反的符號(hào)，那么由中值定理知道存在a<c<b，使得f(c)=0。二分法通過(guò)在每次迭代中簡(jiǎn)單的判斷f(x)在中點(diǎn)x=(a+b)/2處的符號(hào)來(lái)尋求方程的根。如果f(a)和f(x)有相反的符號(hào)那么區(qū)間就被[a,x]代替，否那么就被[x,b]代替。在每次迭代中，包含根的區(qū)間長(zhǎng)度減少一半。即下面我們使用二分法求此方程的一個(gè)數(shù)值解。我們首先要找到一個(gè)區(qū)間，比方(0,5n)，使得函數(shù)在區(qū)間兩端有著不同的符號(hào)。然后即可使用二分法。例1解方程其中a為常數(shù)，n>2為一整數(shù)。顯然，方程的解為程序：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)運(yùn)行結(jié)果：trueroots-4.2394734.1868411.1.2Brent’smethod二分法是一種特殊的括入根算法。Brent通過(guò)逆二次插值方法將括入根方法和二分法結(jié)合起來(lái)。其使用y的二次函數(shù)來(lái)擬合x(chóng)。如果三個(gè)點(diǎn)為(a,f(a)),(b,f(b)),(c,f(c)〕,其中b為當(dāng)前最好的估計(jì)，那么通過(guò)Lagrange多項(xiàng)式插值方法〔y=0)對(duì)方程的根進(jìn)行估計(jì)，在R中，函數(shù)uniroot就是應(yīng)用Brent方法求解一元方程的數(shù)值根。例2應(yīng)用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)輸出結(jié)果：

1

54.2526180.747382224.2526184.1873470.0652709534.1873474.1868410.000505533844.1868414.1868413.032932e-08

Newton方法依賴于f的形狀和初值。該方法從初值開(kāi)始就發(fā)散。運(yùn)行結(jié)果：運(yùn)行結(jié)果：

2.應(yīng)用于求解隨機(jī)優(yōu)化問(wèn)題的蒙特卡羅方法2.1模擬退火算法模擬退火算法來(lái)源于固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時(shí)，固體內(nèi)部粒子隨溫升變?yōu)闊o(wú)序狀，內(nèi)能增大，而徐徐冷卻時(shí)粒子漸趨有序，在每個(gè)溫度都到達(dá)平衡態(tài)，最后在常溫時(shí)到達(dá)基態(tài)，內(nèi)能減為最小。根據(jù)Metropolis準(zhǔn)那么，粒子在溫度T時(shí)趨于平衡的概率為e-ΔE/(kT)，其中E為溫度T時(shí)的內(nèi)能，ΔE為其改變量，k為Boltzmann常數(shù)。用固體退火模擬組合優(yōu)化問(wèn)題，將內(nèi)能E模擬為目標(biāo)函數(shù)值f，溫度T演化成控制參數(shù)t，即得到解組合優(yōu)化問(wèn)題的模擬退火算法：由初始解i和控制參數(shù)初值t開(kāi)始，對(duì)當(dāng)前解重復(fù)“產(chǎn)生新解→計(jì)算目標(biāo)函數(shù)差→接受或舍棄〞的迭代，并逐步衰減t值，算法終止時(shí)的當(dāng)前解即為所得近似最優(yōu)解，這是基于蒙特卡羅迭代求解法的一種啟發(fā)式隨機(jī)搜索過(guò)程。退火過(guò)程由冷卻進(jìn)度表(CoolingSchedule)控制，包括控制參數(shù)的初值t及其衰減因子Δt、每個(gè)t值時(shí)的迭代次數(shù)L和停止條件S給定一些觀察數(shù)據(jù)x，假設(shè)x符合如下高斯分布：求混合高斯分布的三組參數(shù)2.2EM算法問(wèn)題來(lái)源EM算法是個(gè)聚類算法，即根據(jù)給定觀察數(shù)據(jù)自動(dòng)對(duì)數(shù)據(jù)進(jìn)行分類。該混合高斯分布一共有K個(gè)分布，并且對(duì)于每個(gè)觀察到的x，如果我們同時(shí)還知道它屬于K中的哪一個(gè)分布，那么我們可以根據(jù)最大似然估計(jì)求出每個(gè)參數(shù)。結(jié)論：簡(jiǎn)單問(wèn)題特別注意是個(gè)向量，而是個(gè)數(shù)值。表示屬于第k個(gè)高斯分布的觀察數(shù)據(jù)x。實(shí)際問(wèn)題觀察數(shù)據(jù)x屬于哪個(gè)高斯分布是未知的所以要用EM算法來(lái)解決這種實(shí)際問(wèn)題。EM算法過(guò)程：1、用隨機(jī)函數(shù)初始化K個(gè)高斯分布的參數(shù)，同時(shí)保證Expectation2、依次取觀察數(shù)據(jù)x，比較x在K個(gè)高斯函數(shù)中概率的大小，把x歸類到這K個(gè)高斯中概率最大的一個(gè)。