关于算法:R语言用Hessianfree-NelderMead优化方法对数据进行参数估计

35次阅读

共计 3566 个字符,预计需要花费 9 分钟才能阅读完成。

原文链接:http://tecdat.cn/?p=22828 

次要优化办法的疾速概述

咱们介绍次要的优化办法。咱们思考以下问题 .

无导数优化办法

Nelder-Mead 办法是最驰名的无导数办法之一,它只应用 f 的值来搜寻最小值。过程:

  1. 设置初始点 x1,…,xn+1
  2. 对点进行排序,使得 f(x1)≤f(x2)≤⋯≤f(xn+1)。
  3. 计算 xo 作为 x1,…,xn 的中心点。
  4. 反射
  • 计算反射点 xr=xo+α(xo-xn+1)。
  • 如果 f(x1)≤f(xr)<f(xn),那么用 xr 替换 xn+1,转到步骤 2。
  • 否则转到第 5 步。

扩大:

  • 如果 f(xr)<f(x1),那么计算扩大点 xe=xo+γ(xo−xn+1).
  • 如果 f(xe)<f(xr),那么用 xe 替换 xn+1,转到步骤 2。
  • 否则用 xr 替换 xn+1,转到第 2 步。
  • 否则转到第 6 步。

膨胀:

  • 计算膨胀点 xc=xo+β(xo-xn+1).
  • 如果 f(xc)<f(xn+1),那么用 xc 替换 xn+1,进入第 2 步。
  • 否则转到第 7 步.

缩小:

  • 对于 i =2,…,n+1,计算 xi=x1+σ(xi-x1).

Nelder-Mead 办法在 optim 中可用。默认状况下,在 optim 中,α=1,β=1/2,γ=2,σ=1/2。

Hessian-free 优化办法

对于润滑的非线性函数,个别采纳以下办法:部分办法联合直线搜寻工作的计划 xk+1=xk+tkdk,其中部分办法将指定方向 dk,直线搜寻将指定步长 tk∈R。

基准

为了简化优化办法的基准,咱们创立一个函数,用于计算所有优化办法的现实预计办法。

benchfit <- function(data, distr, ...)

β 散布的数值阐明

β 散布的对数似然函数及其梯度

理论值

β 散布的密度由以下公式给出

其中 β 示意 β 函数。咱们记得 β(a,b)=Γ(a)Γ(b)/Γ(a+b)。在这里,一组观测值(x1,…,xn)的对数似然性为

与 a 和 b 无关的梯度为

R 实现

咱们最小化了对数似然的 相同 _数_:实现了梯度的 相同_数_。对数似然和它的梯度都不被输入。

function(par) 
loglikelihood(par, fix.arg ,...)

 样本的随机生成 

#(1) beta 散布
n <- 200
x <- rbeta(n, 3, 3/4)
lnl(c(3, 4), x) #测验
   

``````
hist(x, prob=TRUE)

拟合 Beta 散布

定义控制参数。

list(REPORT=1, maxit=1000)

用默认的优化函数调用,对于不同的优化办法,有梯度和无梯度。

fit(x, "beta", "mle", lower=0,...)

 

在束缚优化的状况下,咱们通过应用对数阻碍容许线性不平等束缚。

应用形态参数 δ1 和 δ2 的 exp/log 变换,来确保形态参数严格为正。

# 取起始值的对数
lapply(default(x, "beta"), log)
#为新的参数化从新定义梯度
exp <- function(par,...) beta(exp(par), obs) * exp(par)
fit(x, distr="beta2", method="mle")

 

# 返回到原始参数化
expopt <- exp(expopt)

而后,咱们提取拟合参数的值、相应的对数似然值和要最小化的函数的计数及其梯度(无论是实践上的梯度还是数值上的近似值)。

数值考察的后果

结果显示在以下表格中。1)没有指定梯度的原始参数(- B 代表有界版本),(2)具备(实在)梯度的原始参数(- B 代表有界版本,- G 代表梯度),(3)没有指定梯度的对数转换参数,(4)具备(实在)梯度的对数转换参数(- G 代表梯度)。

 

 

 咱们绘制了实在值(绿色)和拟合参数(红色)四周的对数似然曲面图。

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 函数来模仿 bootstrap 复制的状况。

boot(fit(x, "beta", method="mle", optim.method="BFGS"))

plot(b1)
abline(v=3, h=3/4, col="red", lwd=1.5)

负二项分布的演示

负二项分布的对数似然函数及其梯度

理论值

负二项分布的 p.m.f. 由以下公式给出

其中 Γ 示意 β 函数。存在另一种示意办法,即 μ =m(1-p)/ p 或等价于 p =m/(m+μ)。因而,一组观测值(x1,…,xn)的对数似然性是

绝对于 m 和 p 的梯度是

R 实现

咱们最小化对数似然性的 相同 _数_:实现梯度的 相同_数_。

  m <- x\[1\]
  p <- x\[2\]
  c(sum(psigamma(obs+m)) - n\*psigamma(m) + n\*log(p),
    m*n/p - sum(obs)/(1-p))
   

样本的随机生成

#(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))

拟合负二项分布

定义控制参数并做基准。

list(trace=0, REPORT=1, maxit=1000)
fit(x, "nbinom", "mle", lower=0)

在束缚优化的状况下,咱们通过应用对数阻碍容许线性不平等束缚。

应用形态参数 δ1 和 δ2 的 exp/log 变换,来确保形态参数严格为正。

# 对起始值进行变换
mu <- size / (size+mu)
arg <- list(size=log(start), prob=log(start/(1-start)))

#为新的参数化从新定义梯度
function(x)
  c(exp(x\[1\]), plogis(x\[2\]))

fit(x, distr="nbinom2", method="mle")

# 返回到原始参数化
expo <- apply(expo, 2, Trans)
   

而后,咱们提取拟合参数的值、相应的对数似然值和要最小化的函数的计数及其梯度(无论是实践上的梯度还是数值上的近似值)。

数值考察的后果

结果显示在以下表格中。1)没有指定梯度的原始参数(- B 代表有界版本),(2)具备(实在)梯度的原始参数(- B 代表有界版本,- G 代表梯度),(3)没有指定梯度的对数转换参数,(4)具备(实在)梯度的对数转换参数(- G 代表梯度)。

 

 咱们绘制了实在值(绿色)和拟合参数(红色)四周的对数似然曲面图。

surface(min.arg=c(5, 0.3), max.arg=c(15, 1), 
         )
points(trueval , pch="x")
   

咱们能够用 bootdist 函数来模仿 bootstrap 复制的状况。

boot(fit(x, "nbinom", method="mle")

 

plot(b1)
abline(v=trueval)
   

论断

基于后面的两个例子,咱们察看到所有的办法都收敛到了同一个点。

然而,不同办法的函数评估(和梯度评估)的后果是十分不同的。此外,指定对数似然性的实在梯度对拟合过程没有任何帮忙,通常会减慢收敛速度。一般来说,最好的办法是规范 BFGS 办法或对参数进行指数变换的 BFGS 办法。因为指数函数是可微的,所以渐进个性仍被保留(通过 Delta 办法),但对于无限样本来说,这可能会产生一个小的偏差。


最受欢迎的见解

1.Matlab 马尔可夫链蒙特卡罗法(MCMC)预计随机稳定率(SV,Stochastic Volatility)模型

2.基于 R 语言的疾病制图中自适应核密度估计的阈值抉择办法

3.WinBUGS 对多元随机稳定率模型:贝叶斯预计与模型比拟

4. R 语言回归中的 hosmer-lemeshow 拟合优度测验

5.matlab 实现 MCMC 的马尔可夫切换 ARMA – GARCH 模型预计

6. R 语言区间数据回归剖析

7. R 语言 WALD 测验 VS 似然比测验

8.python 用线性回归预测股票价格

9. R 语言如何在生存剖析与 Cox 回归中计算 IDI,NRI 指标

正文完
 0