拓端tecdat|R语言用Hessian-free 、Nelder-Mead优化方法对数据进行参数估计
原文链接:http://tecdat.cn/?p=22828
原文出处:拓端数据部落公众号
主要优化方法的快速概述
我们介绍主要的优化方法。我们考虑以下问题 .
无导数优化方法
Nelder-Mead方法是最著名的无导数方法之一,它只使用f的值来搜索最小值。过程:
- 设置初始点x1,...,xn+1
- 对点进行排序,使得f(x1)≤f(x2)≤⋯≤f(xn+1)。
- 计算xo作为x1,...,xn的中心点。
- 反射
- 计算反射点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指标
拓端tecdat|R语言用Hessian-free 、Nelder-Mead优化方法对数据进行参数估计相关推荐
- 拓端tecdat|R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险
最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险 逻辑回归Logistic模型原理和R语言分类预测冠 ...
- 拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解(STL)进行时间序列异常检测
最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...
- 拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...
- 拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例
最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的 房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子 人口众多 独生子女政策: ...
- 拓端tecdat荣获掘金社区入驻新人奖
2021年7月,由掘金发起了"入驻成长礼"颁奖活动.本次活动邀请到知名开发者.服务机构代表等业界人士. 据了解,掘金社区"新入驻创作者礼"主要对已经积累了一定历 ...
- 拓端tecdat荣获2022年度51CTO博主之星
相信技术,传递价值,这是51CTO每一个技术创作者的动力与信念,2022 年度,拓端tecdat 作为新锐的数据分析咨询公司,在51CTO平台上,不断的输出优质的技术文章,分享前沿创新技术,输出最佳生 ...
- 数据分享|R语言关联规则挖掘apriori算法挖掘评估汽车性能数据
全文链接:http://tecdat.cn/?p=32092 我们一般把一件事情发生,对另一件事情也会产生影响的关系叫做关联.而关联分析就是在大量数据中发现项集之间有趣的关联和相关联系(形如" ...
- R语言可视化散点图、气泡图、动态气泡图、数据点重合的散点图、数据点计数图、抖动数据点图、基于lm方法或者loess方法拟合数据点之间的趋势关系曲线、自定义数据点的大小、色彩、添加主标题、副标题、题注
R语言可视化散点图.气泡图.动态气泡图.数据点重合的散点图.数据点计数图.抖动数据点图.基于
- R语言可视化斜率图、扩充图像纵横比为数据标签显示更整齐、ggrepel包来帮忙
R语言可视化斜率图.扩充图像纵横比为数据标签显示更整齐.ggrepel包来帮忙 目录
- R语言使用t.test函数计算两组独立数据的t检验(Independent t-test)
R语言使用t.test函数计算两组独立数据的t检验(Independent t-test) 目录 R语言使用t.test函数计算两组独立数据的t检验(Independent t-test) #仿真数据
最新文章
- Android 打包 aar文件的流程以及aar的引用
- OpenCV(22)SIFT尺度不变特征变换(纯理论)
- Welcome to Swift (苹果官方Swift文档初译与注解三十四)---241~247页(第五章-- 函数)
- linux内核时钟管理,Linux时钟管理透彻分析
- Zabbix学习之路(五)之MySQL监控
- 团队行为心理学读书笔记(8)绩效考核背后的行为心理学
- 真香警告!2021Android高级面试题,挥泪整理面经
- django-celery定时任务以及异步任务and服务器部署并且运行全部过程
- namespace! 报错
- 斗鱼上市首日低开平收 总市值37.3亿美元
- 抓包工具 tcpdump tshark
- 开发中一些常用的css小技巧
- java垃圾回收的方法_java垃圾回收的方法都有哪些
- IOS13以上抓https包,基于win7+Fiddler,操作记录
- 简单的avr c语言程序,avr单片机c语言编程风格介绍 - 全文
- Android软键盘遮挡问题解决
- 蓝凌ekp开发_蓝凌EKP在eclipse中启动报错
- java实验目的_Java实验报告(实验一)
- 今之君子,其责人也详,其待己也廉
- 【STM32利用CuBe MX生成HID设备】1-熟悉软件以及生成一个8键的游戏控制器
热门文章
- 新手android中ListView实现音乐列表
- 设备管理(最近考试有考到,就转一下)
- OCS Inventory NG使用之win平台下的AGENT端安装与信息收集(一)
- htons htonl ntohl ntohs 的区别和作用
- 【TDA4系列】 IPC applications应用举例
- 【翻译】3D Bounding Box Estimation Using Deep Learning and Geometry
- 【Webcam设计】USB摄像头(V4L2接口)的图片采集
- Introduction to Conditional Random Fields
- 数据--第35课 - 创建二叉树
- 数据--第28课 - 进阶星移