拓端tecdat|R語言用Hessian-free 、Nelder-Mead優(yōu)化方法對(duì)數(shù)據(jù)進(jìn)行參數(shù)估計(jì)
原文鏈接：http://tecdat.cn/?p=22828?
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
主要優(yōu)化方法的快速概述
我們介紹主要的優(yōu)化方法。我們考慮以下問題?
.
無導(dǎo)數(shù)優(yōu)化方法
Nelder-Mead方法是最著名的無導(dǎo)數(shù)方法之一，它只使用f的值來搜索最小值。過程：
設(shè)置初始點(diǎn)x1,...,xn+1
對(duì)點(diǎn)進(jìn)行排序，使得f(x1)≤f(x2)≤?≤f(xn+1)。
計(jì)算xo作為x1,...,xn的中心點(diǎn)。
反射
計(jì)算反射點(diǎn)xr=xo+α（xo-xn+1）。
如果f(x1)≤f(xr)<f(xn)，那么用xr替換xn+1，轉(zhuǎn)到步驟2。
否則轉(zhuǎn)到第5步。
擴(kuò)展:
如果f(xr)<f(x1)，那么計(jì)算擴(kuò)展點(diǎn)xe=xo+γ(xo?xn+1).
如果f(xe)<f(xr)，那么用xe替換xn+1，轉(zhuǎn)到步驟2。
否則用xr替換xn+1，轉(zhuǎn)到第2步。
否則轉(zhuǎn)到第6步。
收縮:
計(jì)算收縮點(diǎn)xc=xo+β（xo-xn+1）.
如果f(xc)<f(xn+1)，那么用xc替換xn+1，進(jìn)入第2步。
否則轉(zhuǎn)到第7步.
減少:
對(duì)于i=2,...,n+1，計(jì)算xi=x1+σ（xi-x1）.
Nelder-Mead方法在optim中可用。默認(rèn)情況下，在optim中，α=1，β=1/2，γ=2，σ=1/2。
Hessian-free?優(yōu)化方法
對(duì)于光滑的非線性函數(shù)，一般采用以下方法：局部方法結(jié)合直線搜索工作的方案xk+1=xk+tkdk，其中局部方法將指定方向dk，直線搜索將指定步長(zhǎng)tk∈R。
基準(zhǔn)
為了簡(jiǎn)化優(yōu)化方法的基準(zhǔn)，我們創(chuàng)建一個(gè)函數(shù)，用于計(jì)算所有優(yōu)化方法的理想估計(jì)方法。
benchfit <- function(data, distr, ...)
β分布的數(shù)值說明
β分布的對(duì)數(shù)似然函數(shù)及其梯度
理論值
β分布的密度由以下公式給出
其中β表示β函數(shù)。我們記得β(a,b)=Γ(a)Γ(b)/Γ(a+b)。在這里，一組觀測(cè)值（x1,...,xn）的對(duì)數(shù)似然性為
與a和b有關(guān)的梯度為
R實(shí)現(xiàn)
我們最小化了對(duì)數(shù)似然的相反數(shù)：實(shí)現(xiàn)了梯度的相反數(shù)。對(duì)數(shù)似然和它的梯度都不被輸出。
function(par)loglikelihood(par, fix.arg ,...)
?樣本的隨機(jī)生成?
#(1) beta分布n <- 200x <- rbeta(n, 3, 3/4)lnl(c(3, 4), x) #檢驗(yàn)
?hist(x, prob=TRUE)
擬合Beta分布
定義控制參數(shù)。
list(REPORT=1, maxit=1000)
用默認(rèn)的優(yōu)化函數(shù)調(diào)用，對(duì)于不同的優(yōu)化方法，有梯度和無梯度。
fit(x, "beta", "mle", lower=0,...)
?
在約束優(yōu)化的情況下，我們通過使用對(duì)數(shù)障礙允許線性不平等約束。
使用形狀參數(shù)δ1和δ2的exp/log變換，來確保形狀參數(shù)嚴(yán)格為正。
#取起始值的對(duì)數(shù)lapply(default(x, "beta"), log)#為新的參數(shù)化重新定義梯度exp <- function(par,...) beta(exp(par), obs) * exp(par)fit(x, distr="beta2", method="mle")
?
#返回到原始參數(shù)化expopt <- exp(expopt)
然后，我們提取擬合參數(shù)的值、相應(yīng)的對(duì)數(shù)似然值和要最小化的函數(shù)的計(jì)數(shù)及其梯度（無論是理論上的梯度還是數(shù)值上的近似值）。
數(shù)值調(diào)查的結(jié)果
結(jié)果顯示在以下表格中。1）沒有指定梯度的原始參數(shù)（-B代表有界版本），（2）具有（真實(shí)）梯度的原始參數(shù)（-B代表有界版本，-G代表梯度），（3）沒有指定梯度的對(duì)數(shù)轉(zhuǎn)換參數(shù)，（4）具有（真實(shí)）梯度的對(duì)數(shù)轉(zhuǎn)換參數(shù)（-G代表梯度）。
?
?
?我們繪制了真實(shí)值（綠色）和擬合參數(shù)（紅色）周圍的對(duì)數(shù)似然曲面圖。
llsurface(min.arg=c(0.1, 0.1), max.arg=c(7, 3),plot.arg=c("shape1", "shape2"), nlev=25,plot.np=50, data=x, distr="beta", back.col = FALSE)points(unconstropt[1,"BFGS"], unconstropt[2,"BFGS"], pch="+", col="red")points(3, 3/4, pch="x", col="green")
?
我們可以用bootdist函數(shù)來模擬bootstrap?復(fù)制的情況。
boot(fit(x, "beta", method="mle", optim.method="BFGS"))
plot(b1)abline(v=3, h=3/4, col="red", lwd=1.5)
負(fù)二項(xiàng)分布的演示
負(fù)二項(xiàng)分布的對(duì)數(shù)似然函數(shù)及其梯度
理論值
負(fù)二項(xiàng)分布的p.m.f.由以下公式給出
其中Γ表示β函數(shù)。存在另一種表示方法，即μ=m(1-p)/p或等價(jià)于p=m/(m+μ)。因此，一組觀測(cè)值（x1,...,xn）的對(duì)數(shù)似然性是
相對(duì)于m和p的梯度是
R實(shí)現(xiàn)
我們最小化對(duì)數(shù)似然性的相反數(shù)：實(shí)現(xiàn)梯度的相反數(shù)。
m <- x[1]p <- x[2]c(sum(psigamma(obs+m)) - n*psigamma(m) + n*log(p),m*n/p - sum(obs)/(1-p))
?
樣本的隨機(jī)生成
#(1) β分布trueval <- c("size"=10, "prob"=3/4, "mu"=10/3)x <- rnbinom(n, trueval["size"], trueval["prob"])hist(x, prob=TRUE, ylim=c(0, .3))
擬合負(fù)二項(xiàng)分布
定義控制參數(shù)并做基準(zhǔn)。
list(trace=0, REPORT=1, maxit=1000)fit(x, "nbinom", "mle", lower=0)
在約束優(yōu)化的情況下，我們通過使用對(duì)數(shù)障礙允許線性不平等約束。
使用形狀參數(shù)δ1和δ2的exp/log變換，來確保形狀參數(shù)嚴(yán)格為正。
#對(duì)起始值進(jìn)行變換mu <- size / (size+mu)arg <- list(size=log(start), prob=log(start/(1-start)))#為新的參數(shù)化重新定義梯度function(x)c(exp(x[1]), plogis(x[2]))fit(x, distr="nbinom2", method="mle")
#返回到原始參數(shù)化expo <- apply(expo, 2, Trans)
?
然后，我們提取擬合參數(shù)的值、相應(yīng)的對(duì)數(shù)似然值和要最小化的函數(shù)的計(jì)數(shù)及其梯度（無論是理論上的梯度還是數(shù)值上的近似值）。
數(shù)值調(diào)查的結(jié)果
結(jié)果顯示在以下表格中。1）沒有指定梯度的原始參數(shù)（-B代表有界版本），（2）具有（真實(shí)）梯度的原始參數(shù)（-B代表有界版本，-G代表梯度），（3）沒有指定梯度的對(duì)數(shù)轉(zhuǎn)換參數(shù)，（4）具有（真實(shí)）梯度的對(duì)數(shù)轉(zhuǎn)換參數(shù)（-G代表梯度）。
?
?我們繪制了真實(shí)值（綠色）和擬合參數(shù)（紅色）周圍的對(duì)數(shù)似然曲面圖。
surface(min.arg=c(5, 0.3), max.arg=c(15, 1),)points(trueval , pch="x")
?
我們可以用bootdist函數(shù)來模擬bootstrap?復(fù)制的情況。
boot(fit(x, "nbinom", method="mle")
?
plot(b1)abline(v=trueval)
?
結(jié)論
基于前面的兩個(gè)例子，我們觀察到所有的方法都收斂到了同一個(gè)點(diǎn)。
然而，不同方法的函數(shù)評(píng)價(jià)（和梯度評(píng)價(jià)）的結(jié)果是非常不同的。此外，指定對(duì)數(shù)似然性的真實(shí)梯度對(duì)擬合過程沒有任何幫助，通常會(huì)減慢收斂速度。一般來說，最好的方法是標(biāo)準(zhǔn)BFGS方法或?qū)?shù)進(jìn)行指數(shù)變換的BFGS方法。由于指數(shù)函數(shù)是可微的，所以漸進(jìn)特性仍被保留（通過Delta方法），但對(duì)于有限樣本來說，這可能會(huì)產(chǎn)生一個(gè)小的偏差。
最受歡迎的見解
1.Matlab馬爾可夫鏈蒙特卡羅法（MCMC）估計(jì)隨機(jī)波動(dòng)率（SV，Stochastic Volatility） 模型
2.基于R語言的疾病制圖中自適應(yīng)核密度估計(jì)的閾值選擇方法
3.WinBUGS對(duì)多元隨機(jī)波動(dòng)率模型：貝葉斯估計(jì)與模型比較
4.R語言回歸中的hosmer-lemeshow擬合優(yōu)度檢驗(yàn)
5.matlab實(shí)現(xiàn)MCMC的馬爾可夫切換ARMA – GARCH模型估計(jì)
6.R語言區(qū)間數(shù)據(jù)回歸分析
7.R語言WALD檢驗(yàn) VS 似然比檢驗(yàn)
8.python用線性回歸預(yù)測(cè)股票價(jià)格
9.R語言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)
