热门标签 | HotTags
当前位置:  开发笔记 > 编程语言 > 正文

opcr参数ua_R语言非参时间序列(一):MetropolisHastings算法——不同转移核的思考…

MCMC在金融工程中有着重要作用,它可以用来仿真产品价格分布以及波动,可以用来求解模型中的参数。Metropolis-Hastings算法是一种常用的MCMC算法,本文是我在使

MCMC在金融工程中有着重要作用,它可以用来仿真产品价格分布以及波动,可以用来求解模型中的参数。Metropolis-Hastings算法是一种常用的MCMC算法,本文是我在使用MCMC时候的一些思考,和大家一起分享。阅读本文需要对MCMC有一定的认识,本文默认读者熟悉MCMC,并在此基础上对一些技术细节给出了自己的思考。文中算例采用R语言实现,每一个算例都非常简单,但是很具有代表性,如有不正确的地方,希望大家多多批评,多多指正。




0、Metropolis-Hastings算法回顾

算法的作用,是在已知某种概率密度函数(pdf)的情况下,生成一族服从这种概率密度函数的随机数。这里的概率分布函数应该是比较奇葩的,不然我们也不会用到MH算法。或者,是在已知某种函数的情况下,生成一族随机数,这族随机数的直方图的中点连线就是这个函数。说白了,这个函数就是非正则pdf。我这里说的pdf,可以是连续的,也可以是离散的。

转移核Q: 这个转移核的作用可以以某一个概率把i变为j,这个概率我们表示为Q(i,j)。例如某一个转移核能够以3/14的概率把7变成12.8,那么Q(7,12.8)=3/14。如果是离散的转移核,那我们很好理解,它可以是马尔科夫过程,Q就是马氏链的转移矩阵。如果是连续的转移核,那Q就应该是概率密度函数。

仿真目标P:算法的任务是生成一组随机数{




那么




我们这里使用MH算法,并采用一种随机游走的建议转移核。为什么要使用随机游走的建议转移核?事实上我们用什么核都行,这里我们先使用一种典型转移核,如下:




9658582de7a3d2ee76a241ab2146adc7.png

不带吸收壁的9×9对称马尔科夫随机游走转移矩阵

library(ggplot2)
set.seed(1234)#随机数种子
P<-c(9,1,8,2,7,3,6,4,5)/45#我们的目标概率分布
theta1<-2#初始值
iter<-10000#迭代次数
chain<-vector(length = iter)#要生成的一族随机数
for (i in seq_len(iter)){
if (theta1==1){
ifelse(runif(1)<0.5,theta<-theta1,theta<-theta1+1)
}else if(theta1==11){
ifelse(runif(1)<0.5,theta<-theta1,theta<-theta1-1)
}else{
ifelse(runif(1)<0.5,theta<-theta1+1,theta<-theta1-1)
}#这个部分其实是转移核Q,对theta1进行转移,Q服从无吸收壁随机游走
u<-runif(1)
a<-min(1,P[theta]*0.5/P[theta1]/0.5)#注意:Q(i,j)=0.5,Q(j,i)=0.5
if (u theta1<-theta
}
chain[i]<-theta1
}
chain1<-as.data.frame(chain)
ggplot(chain1, aes(x=chain))+
geom_histogram(bins = 9)#直方图



448afd93b37651040bf0e202aa3c2b4a.png

仿真结果

2、对(非正则)离散分布仿真,采用一个很奇怪的转移核(体现出转移核的随意性)

假设我们要随机生成的随机变量X,需要满足如下的离散分布:




812a22410aadbb8d68ec7b40dcd21099.png

我们说过转移核可以随意选取,到底可以多随意呢?我们不妨采用如下的一个转移核:




3eb7c0de71dc538fa1e2753305ea5ebb.png

瞎编的数字代表我们可以随意找一个转移核

数字是我瞎编的,我们来试一试:


library(ggplot2)
set.seed(1234)
P<-c(1,3,2)/6
theta1<-2#初始值
iter<-100000#迭代次数
Q<-matrix(c(0.234,0.145,0.621,0.66,0,0.34,0.1111,0.2222,0.6667),nr=3,byrow = F)
chain<-vector(length = iter)#要生成的一族随机数
for (i in seq_len(iter)){
if (theta1==1){
theta<-sample(c(1:3),prob = c(0.234,0.145,0.621),size = 1)
}else if(theta1==2){
theta<-sample(c(1:3),prob = c(0.66,0,0.34),size = 1)
}else{
theta<-sample(c(1:3),prob = c(0.1111,0.2222,0.6667),size = 1)
}#这个部分其实是转移核Q,对theta1进行转移
u<-runif(1)
a<-min(1,P[theta]*Q[theta1,theta]/P[theta1]/Q[theta,theta1])
if (u theta1<-theta
}
chain[i]<-theta1
}
chain1<-as.data.frame(x1<-c(14,7,22)[chain])
ggplot(chain1, aes(x=x1))+
geom_histogram(bins = 3)



34bb7402c649238ac606b291f8811741.png

仿真结果

3、对(非正则)离散分布仿真,采用退化转移核(为后面的内容做铺垫)

假设我们要随机生成的随机变量X,需要满足如下的离散分布:




f6a386f8f01d1676604b2817851e21b2.png

采用退化转移核:




762c8e24cb33829519569893e89ba2f8.png

9*9退化转移核

退化转移核作用是对于任意的i,j,Q(i,j)=1/9,我们于是发现一个问题Q(i,j)=Q(j,i)=c了。


library(ggplot2)
set.seed(1234)#随机数种子
P<-c(9,1,8,2,7,3,6,4,5)/45#我们的目标概率分布
theta1<-1#初始值
iter<-10000#迭代次数
chain<-vector(length = iter)#要生成的一族随机数
for (i in seq_len(iter)){
theta<-sample(1:9,size=1)#这个部分其实是转移核Q,对theta1进行转移
u<-runif(1)
a<-min(1,P[theta]/P[theta1])#注意:Q(i,j)=Q(j,i)=c,已经约去
if (u theta1<-theta
}
chain[i]<-theta1
}
chain1<-as.data.frame(chain)
ggplot(chain1, aes(x=chain))+
geom_histogram(bins = 9)#直方图



3f422698c5feaad5622c36226101407c.png

4、对(非正则)连续分布仿真,采用连续转移核

这里才是我们讨论的重点,连续转移核和离散转移核其实是一样的,只不过i和j可以去连续数值而已。举个例子,假如我们现在在处在位置是8.18,假如我们的转移核是连续的,并且是正态分布(均值为8.18,假设标准差为3),那么我们可以计算一下从8.18转移到13的概率Q(8.18,13)=0.4459,也就是说它将以44.59%的概率转移到13,那向左移动到2.22呢?Q(8.18,2.21)=0.023,也就是说它将以2.3%的概率移动到2.21。我们从图上可以感知一下,什么是连续的转移核:




108a74d99d565f17e6e078c1e2fb309f.png

蓝线所在位置是8.18,它将以进行正态分布(标准差假设为3)为转移核进行转移



64409b6107893ed3de8c4faa84f839cb.png

现在蓝线转移到绿线(x=3)了,我们计算出来Q(8.18,13)=0.4459

那么我们可以发现转移是对称的,什么叫对称呢,就是Q(i,j)=Q(j,i),也就是Q(8.18,13)=Q(13,8.18),我们可以从图上一眼看出来:




0bf767ccdeeba342429b708cf7307fb5.png

从蓝线转移到绿线和从绿线转移到蓝线有着相同的概率分布

那么我们来试一试:假如我要实现X~





我们假设此时我们在位置1(红线),我们采用正态转移核(标准差为3),可是我们不能转移到小于0的的地方,所以我们必须对转移核进行截断!形成一种截断转移核。



5c00ab28c23bef08827c07db09d688d0.png

我们假设此时我们在位置3.4(红线),我们采用正态转移核(标准差为3),可是我们不能转移到小于0的的地方,所以我们必须对转移核进行截断!形成一种截断转移核。

我们来试一试:

实现X~







d6e3b21f23da8f9e774861f27b7a731d.png

上图是MH结果,下图是真实的pdf

5、对(非正则)连续分布仿真,采用连续退化转移核

5其实就是对3的衍生,退化转移核在连续分布下其实就是均匀分布,但是有一个问题,这个均匀分布的取值范围是什么?我们首先来试一试,取值范围是【x-3,x+3】的情形:

此时,Q(i,j)=Q(j,i),


library(ggplot2)
set.seed(1234)#随机数种子
P<-function(x) 0.5*x^2*exp(-x)
theta1<-rnorm(1,mean = 5)
iter<-20000
chain<-vector(length = iter)
for (i in seq_len(iter)){
theta<-runif(1,min = theta1-3,max=theta1+3)
if (theta<0){
theta<-theta1
}
u<-runif(1)
a<-min(P(theta)/P(theta1),1)#建议分布是概率密度函数
if (u theta1<-theta
}
chain[i]<-theta1
}
chain1<-as.data.frame(chain)
ggplot(chain1, aes(x=chain))+
geom_histogram(bins = 100)



16dcf91be98ca5b4d33c739bfeb31815.png

6、对(非正则)连续分布仿真,转移核已经退化,没有转移核了

什么意思?如果5中,我们的均匀分布无限延长就是没有转移核了,我们会取一个大于0的随机点,我们先假设是0到1000,我们先试一试。


library(ggplot2)
set.seed(1234)#随机数种子
P<-function(x) 0.5*x^2*exp(-x)
theta1<-4
iter<-200000
chain<-vector(length = iter)
for (i in seq_len(iter)){
theta<-runif(1,min =0,max=1000)
if (theta<0){
theta<-theta1
}
u<-runif(1)
a<-min(P(theta)/P(theta1),1)#建议分布是概率密度函数
if (u theta1<-theta
}
chain[i]<-theta1
}
chain1<-as.data.frame(chain)
ggplot(chain1, aes(x=chain))+
geom_histogram(bins = 100)



934b4c28f2cfbb3058c1bef48d0d8a22.png

问题出在20万次的仿真对于0,1000这个随机转移区间来说远远不够,也体现出了这种方法效率的低下


有点困了,先写这么多吧。star19950818@foxmail.com。




推荐阅读
  • 本文详细介绍了在 Vue.js 前端框架中集成 vue-i18n 插件以实现多语言支持的方法。通过具体的配置步骤和示例代码,帮助开发者快速掌握如何在项目中实现国际化功能,提升用户体验。同时,文章还探讨了常见的多语言切换问题及解决方案,为开发人员提供了实用的参考。 ... [详细]
  • 探索偶数次幂二项式系数的求和方法及其数学意义 ... [详细]
  • 本文详细探讨了Laravel框架中路由和控制器的操作实例,旨在为开发者提供深入的理解和实用的参考。通过具体的示例,读者可以更好地掌握如何在Laravel中配置和使用路由及控制器。此外,文章还提供了丰富的PHP编程技巧和最佳实践,帮助读者提升开发效率和代码质量。 ... [详细]
  • PHP 数组逆序排列方法及常用排序函数详解 ... [详细]
  • 本文深入探讨了二叉树路径和问题的算法优化方法。具体而言,给定一棵二叉树,需要找出所有从根节点到叶节点的路径,其中各节点值的总和等于指定的目标值。通过详细分析和优化,提出了一种高效的解决方案,并通过多个样例验证了其有效性和性能。 ... [详细]
  • 本文详细解析了九度编程平台上的斐波那契数列高效算法挑战(题目编号:1387)。该挑战要求在1秒的时间限制和32兆的内存限制下,设计出高效的斐波那契数列计算方法。通过多种算法的对比和性能分析,本文提供了优化方案,帮助参赛者在限定资源条件下实现高效计算。 ... [详细]
  • Python与R语言在功能和应用场景上各有优势。尽管R语言在统计分析和数据可视化方面具有更强的专业性,但Python作为一种通用编程语言,适用于更广泛的领域,包括Web开发、自动化脚本和机器学习等。对于初学者而言,Python的学习曲线更为平缓,上手更加容易。此外,Python拥有庞大的社区支持和丰富的第三方库,使其在实际应用中更具灵活性和扩展性。 ... [详细]
  • 本文作为“实现简易版Spring系列”的第五篇,继前文深入探讨了Spring框架的核心技术之一——控制反转(IoC)之后,将重点转向另一个关键技术——面向切面编程(AOP)。对于使用Spring框架进行开发的开发者来说,AOP是一个不可或缺的概念。了解AOP的背景及其基本原理,对于掌握这一技术至关重要。本文将通过具体示例,详细解析AOP的实现机制,帮助读者更好地理解和应用这一技术。 ... [详细]
  • JavaScript最初并非设计为纯粹的面向对象编程(OOP)语言,因为直到ES5标准中仍未引入类的概念。然而,随着ES6的发布,JavaScript正式引入了类的语法,使得开发者能够更加直观地实现继承机制。本文将深入探讨JavaScript中多样的继承实现方法,包括原型链、寄生组合式继承等技术,并分析它们的优缺点及适用场景。 ... [详细]
  • 宏基因组学经典文献重现(一):利用ggplot2进行散点图可视化分析
    宏基因组学经典文献重现(一):利用ggplot2进行散点图可视化分析 ... [详细]
  • 本文探讨了通过JavaScript实现密码强度检测功能的优化方法。密码强度分为三个等级:使用单一类型字符(如纯数字或纯字母)且长度在1-6位的密码为“弱”;包含两种类型字符(例如字母与数字组合)的密码为“中等”;而长度超过12位并结合三种或四种类型字符(如大小写字母、数字及特殊符号)的密码则被认定为“强”。此外,文章还介绍了如何利用正则表达式和条件判断语句来高效地评估密码强度,确保用户账户的安全性。 ... [详细]
  • Python正则表达式详解:掌握数量词用法轻松上手
    Python正则表达式详解:掌握数量词用法轻松上手 ... [详细]
  • PHP编程中的命名规则与最佳实践 ... [详细]
  • Kafka核心理论问题汇编【持续更新中】
    本文汇总了Kafka的核心理论问题,涵盖了常见的技术难点和解决方案。内容将持续更新,旨在为开发者提供全面的参考。文章源自博客园,作者呱嗒呱嗒,转载时请注明出处。 ... [详细]
  • 在一系列的学习与实践后,Jsoup学习笔记系列即将进入尾声。本文详细介绍了如何使用Jsoup实现从Saz文件到Csv格式的数据解析功能。未来,计划将此功能进一步封装,开发成具有用户界面的独立应用程序,以增强其实用性和便捷性。对于希望深入掌握Jsoup技术的开发者,本文提供了宝贵的参考和实践案例。 ... [详细]
author-avatar
PHP1.CN | 中国最专业的PHP中文社区 | DevBox开发工具箱 | json解析格式化 |PHP资讯 | PHP教程 | 数据库技术 | 服务器技术 | 前端开发技术 | PHP框架 | 开发工具 | 在线工具
Copyright © 1998 - 2020 PHP1.CN. All Rights Reserved | 京公网安备 11010802041100号 | 京ICP备19059560号-4 | PHP1.CN 第一PHP社区 版权所有