蒙特卡洛方法及在一些统计模型中的应用

来源 :山东大学 | 被引量 : 0次 | 上传用户:lixuantea
下载到本地 , 更方便阅读
声明 : 本文档内容版权归属内容提供方 , 如果您对本文有版权争议 , 可与客服联系进行内容授权或下架
论文部分内容阅读
随着当代计算机技术的发展,蒙特卡洛方法(又称为统计模拟方法)已经成为人们解决复杂统计模型和高维问题的主要工具.特别是经过最近几十年的发展,蒙特卡洛方法已经极大地拓展了人们的科学视野,不管是从早期计算物理学的发展,还是到现代的计算生物学,都使人们可以研究许多极其复杂的系统.人们第一次系统的应用蒙特卡洛方法解决世纪科学问题是在二十世纪四十年代.由于世界上第一台可编程超级计算机就在这时候诞生,为了能够更好地利用这台电子计算机,一批物理学和统计学家(S. Ulam, J. von Neumann, N. Metropolis, E. Fermi等)发明了一种依赖统计抽样解决复杂问题的方法,并成功将其应用于原子弹设计和薛定谔方程特征根的数值计算.二十世纪五十年代,一批统计物理学家(N. Metropolis. A. Rosenbluth, M. Rosen-bluth, A. Teller和E. Teller等)第一次引入了一种基于马尔科夫链的动态蒙特卡洛方法,随后被用来模拟一些复杂的物理系统,包括Ising模型spin glass模型等.这种马尔科夫链蒙特卡洛方法被称为Metropolis算法,它是也是第一个基于迭代模拟的抽样方法Hastings在1970年将其推广得到Metropolis-Hastings算法,这个算法形式及其简单.这里我们假设目标分布为π,初始值为X0,这个迭代算法从Xu到Xn+1的抽样包含如下两步:Step1.从一个任意给定的建议转移概率密度q(Xn,y)得到“建议抽样”y.Step2.以概率接受“建议抽样”y.并设Xn+1=y:否则,以1-α的概率拒绝“建议抽样”y,并设Xn+1=Xn.二十世纪八十年代,统计学家和计算机学家将蒙特卡洛方法应用于一系列实际应用模型,其中包括组和最优化问题(模拟退火算法[57;70]),非参数统计推断(jackknife和Bootstrap [24;17]),缺失数据下似然函数的计算(EM算法[19;99]和数据扩张算法[101;70;99]),统计遗产算法,贝叶斯建模和贝叶斯计算[100;99]等等,特别是进入二十世纪九十年代,蒙特卡洛方法在计算生物学中扮演了十分重要的角色.并成功解决了基元序列识别(sequence motif identification [70])问题和分析复杂谱系.如今,蒙特卡洛方法已经在生物学,化学,计算机科学,经济与金融,气象学,材料科学,物理学,统计学等等很多领域扮演重要角色,并被广泛应用.作为一种特殊的蒙特卡洛方法,马尔科夫链蒙特卡洛方法由于可以适用于一些更复杂的统计模型和一些更微观的分子模型,最近越来越受到统计学家的关注,并取得了长足发展.马尔科夫链蒙特卡洛方法的最核心的步骤就是抽样,给定目标函数(通常这个目标函数非常复杂多峰甚至是高维目标分布),如何高效的从这个目标分布中抽样成为其它一切问题(例如,参数估计,积分,最优化问题)的关键步骤.随机游走Metropolis算法作为Metropolis-Hastings算法的一种特例,由于其简单明确的形式,被人们广泛应用于实际中.对于同样的目标分布π,随机游走Metropolis算法选择对称的建议分布(即Metropolis-Hastings算法中q(x,y)=q(y,x)),从而Metropolis比率α变为α(Xn,y)=min{1mπ(y)/π(Xn)}、虽然随机游走Metropolis算法从理论上可以对任意给定目标分布π进行抽样,甚至与目标分布维数无关,可是、算法本质上是一个局部更新算法(local update)对于实际应用中的大多数分布,它们都有着复杂的解析形式(多峰或者高维),从而导致了随机游走Metropolis算法极低的抽样效率,甚至在有限的迭代步数内达不到平稳分布.这就是蒙特卡洛抽样算法中所谓的“局部陷阱”(local-trap)问题.从而,设计一个好的马尔科夫链蒙特卡洛抽样方法,并克服“局部陷阱”问题,已经越来越引起人们的关注,并被认为是马尔科夫链蒙特卡洛方法的重要研究方向.最近几十年,为了克服“局部陷阱问题”,许多高级马尔科夫链蒙特卡洛抽样方法先后被人们提出,并在特定的情况下.有它们自己的优势.其中包括模拟回火(simulated tempering)[75;35;70]; Metropolis-Adjusted Langevin算法[91];并行回火算法(parallel tempering)[36;70]自适应拒绝Metropolis算法(adaptive rejection metropolis sampling)[41;40];多测试点Metropolis算法(Multiple-Try Metropolis algorithm)[71;70]进化Monte Carlo (evolutionary Monte Carlo)[66;67;64]动态加权Monte Carlo (dynamic weighting Monte Carlo)[108;63;70;64]; Equi-Energy抽样(Equi-Energy Sampling)[58;64]等等.本文的第一章,提出了一种新的马尔科夫链蒙特卡洛抽样方法:B样条Metropolis-Hastings算法.这种抽样方法是对传统的随机游走Metropolis算法的推广.由于随机游走Metropolis算法采用对称的建议分布(通常人们选择正态分布作为它的建议分布),从而算法的抽样效率就完全依赖于建议分布的方差σ2的大小.过大的建议方差,会使得Metropolis算法拒绝多数的建议抽样.从而使得抽样样本重复在原来样本点;而过小的的建议方差虽然接受了多数建议抽样,可是抽样大多集中在局部区域,从而大大减弱了算法的收敛效率.受到重要抽样的启发[70],我们采用B样条在局部范围近似目标分布,从而可以在不降低算法接受概率的情况下,做大步长建议抽样,这就是B样条Metropolis-Hastings算法的思想.B样条取常见的三次样条作为建议分布:其中n是三次B样条基函数的个数,r定义了y的区域,cTx=(c1,x,…,cn,x)∈[0,+∞]n为B样条函数的系数,BxT(y)=(B1,x(y),…,Bn,x(y))为B样条基函数.下面的定理可以用来高效的从B样条建议分布中抽样:定理1.2.1(1)设B[a,b](x)为在区间[a,b]上等间距节点的三次B样条基函数,定义连续随机变量η:其中{Ui}i=14为独立同分布[0,1]上均匀分布(U[0,1]),则η的概率密度分布函数为B[a,b](x).(2)在建议分布q(Xn,·)中,bxT(y)=(B1,x(y),…,Bn,x(y))的归一化常数是Cxn={1/(24),1/2,(23)/(24),1,…,1,(23)/(24),1/2,1/(24)},即合成法抽样中的系数为:ω={1/(24)c1,Xn,1/2c2,Xn,(23)/(24)c3,Xn,c4,Xn,…cm-3,Xn,(23)/(24)cm-2,Xn,1/2cm-1,Xn,1/(24)cm,Xn}.下面的定理给出了B样条Metropolis-Hastings算法的抽样效率,定理1.2.2设目标函数为π(x),B样条函数的范围为r,如果我们假设B样条的节点数m足够多.则B样条Metropolis-Hastings算法的接受概率(即,马尔科夫链在无限长时间内的平均接受概率)ηr可以写为:为了展示B样条Metropolis-Hastings算法的优势,我们在第一章还给出了高维情况下的B样条Metropolis-Hastings算法的推广,即B样条Hit-and-Run算法和B样条Gibbs抽样.模拟数据和真实数据实验中.我们比较了B样条Metropolis-Hastings算法同其他一些抽样算法(随机游走Metropolis算法,自适应拒绝Metropolis算法,并行回火算法等)的效率,发现B样条Metropolis-Hastings算法有很大的优势,随着最近计算机技术的发展,贝叶斯分析方法已经在医学,社会科学,经济金融等领域有了广泛的应用:并逐渐展示了贝叶斯方法在各个领域的优势.贝叶斯方法最大的特点是将未知参数的看做随机变量,从而有了参数的分布,并给出区间估计.随着马尔科夫链蒙特卡洛技术的发展,Tanner和Wong [101]第一次应用贝叶斯方法解决缺失数据问题,并取得很好的结果.这就是现在被广泛应用的数据扩张算法(Data Augmentation),它的思想与EM算法[19]类似,但是在贝叶斯框架下分析问题.在本文的第二章,我们考虑了混合正态模型参数的贝叶斯分析.这里,我们考虑两个正态分布的混合:其中φ(·;μ、σ)为正态分布的密度函数,y1,…,yn样本值,θ=(a,μ1,σ21,μ2,σ22)T为待估参数.从而我们可以写出未知参数的后验分布:其中π0(θ)取无信息先验分布Little et al和Chen在文献[68]和文献[13]中提出应用EM算法和数据扩张算法分析上述模型.可是,当目标分布为高维情况或者非常粗糙时,蒙特卡洛抽样方法容易遇到“局部陷阱”问题,从而准确地计算参数的后验分布.Kou在文献[58]中首次提出了Equi-Energy抽样,在一定程度上可以克服上述蒙特卡洛方法这些缺点.第二章我们就应用Equi-Energy抽样方法计算模型参数的后验分布,在随后的模拟实验中我们发现,Equi-Energy抽样方法,可以很准确的分析模型参数的后验分布.随着现代金融市场的到来,人们可以掌握到大量的高维数据,统计深度[73]作为多元非参数统计和稳健统计的重要工具越来越多地引起人们的关注.尽管如此,数据深度,特别是广泛应用的投影深度(projection depth)[73;119;120;115]的计算问题却一直迟迟落后于统计深度理论的发展,并阻碍投影深度在实践中的广泛应用[116].在高维情况下,还没有一种高效的算法计算投影深度和它的诱导估计(例如Stahel-Donoho估计).Zuo等在文献[116]给出了一种投影深度的精确计算算法,可是这个算法只能计算二维数据.第三章,我们提出应用模拟退火算法[57;70]计算投影深度.投影深度的计算一个关键步骤就是outlyingness的计算,它可以看做空间的一个全局最优化问题:其中O1(·,·)为一维outlyingnessθ=(θ1,…,θd)∈[0,π]d为高维球坐标变换.于是我们定义最小化的目标函数:h(θ)=-O1(u’(θ))x,u(θ)Xn),并将它引入到模拟退火算法中从而算法的最小值点就是我们希望求的outlyingness的最大值点.随后的数值模拟(二维,四维,七维模拟数据)结果显示,模拟退火算法相比较其他近似算法有较高的计算效率和较高的计算精度.Kendall’s τ和Spearman’s ρ是最常见的非参数相关性度量.众所周知,对于不同的联合密度分布,Kendall’s τ和Spearman’s ρ所计算的相关性度量值不相同,因为它们度量了随机变量间的不同的相关性结构Durbin和Stuart在文献[21]中第一次给出了Kendall’s τ和Spearman’s ρ的关系边界,可是,这个边界并不精确,而且没有给出达到精确边界时候的随机变量的分布.在本文第四章,我们首先将Kendall’s τ和Spearman’s ρ的边界问题看做一个约束条件的全局最优化问题:然后,我们应用模拟退火算法找出了Kendall’s τ和Spearman’s ρ的精确的边界,这个精确的边界可以看做是Durbin和Stuart所给出边界的改进,而且给出了达到精确边界的随机变量的分布.下面是τ∈[-1,0]时候,Kendall’s τ和Spearman’s ρ的改进后的下界:其中n=2,3.4,…,并且τ∈[2/n+1-1,2/n-1],更进一步,我们还可以得到τ在每个小区间[击-1.2/n-1]上ρ对τ的边界的导数:其中τ∈[2/n+1-1,2/n-1],n=2,3,4,5,…,并且α∈[0,1].最近几年.降维问题是统计界的热门问题.解决的办法很多,其中充分降维理论就是解决此类问题的有效途径之一.充分降维的方法大都将问题归结为对称阵的主成分分析,从而对称阵的主成分分析就显得异常重要.Li在文献[61]中应用稀疏主成分分析应用到充分降维问题,但是没有对特征值进行稀疏.zhu在文献[113]中将对称阵的谱分解的估计归为优化问题,并进一步利用LASSO技术[102]对稀疏特征值进行估计,可是没有考虑特征向量的稀疏问题.本文第五章,我们将文献[113]的结果进一步推广应用LASSO技术同时对特征值和特征向量惩罚并得到相应的估计:其中λ=(λ1,…,λp)T为对称阵A的特征值,ζ=(ζ1…ζp)和η=(η1…ηp).λi是A的第i个特征根,且满足|λl|≥|λ2|…≥|λp|≥0,同时,ζ和η是A对应着λi的特征向量.πjk=|ηjk|-γ,ηjk为ηjk的(?)速度收敛的相合估计.随后给出了估计的渐进性质,定理5.3.1(Oracle性质)假设l。(?)→0,lnn(r-1)/2→∞, mn/(?)→0,mnn(r-1)/2→∞.那么,稀疏的主成分估计满足以下条件,(1)模型选择的一致性:当n→∞时,有A(?)A,且对于任意的i,有ηi(?)ηi,其中A={i|ηi≠0},ηi={j|ηij≠0},且(?)表示依概率收敛.(2)渐近正态性:其中W0是(?)(An-A)的极限分布,且(?)表示依概率收敛.最后在模拟实验中.我们重复产生多次抽样样本,基于蒙特卡洛的模拟结果显示,在充分降维问题中稀疏主成分分析在稀疏特征值和稀疏特征向量选择上的优势.
其他文献
<正>英国Tamicare公司开发出世界首创的3D打印纺织成品技术,并冠以"Cosyflex"品牌名称。Cosyflex技术可在同一基面上使用各种纤维打印出各种特征(如图案、压花、穿孔、无缝、
期刊
<正>丁光训主教不但在思想理论建设上为中国基督教的健康发展做出了重大贡献,而且在教会的自身建设实践中也做出了重大贡献。习近平总书记在中央统战工作会议上指出,积极引导
2016年中国乘用车销量超过2400万辆,创全球历史新高,连续八年蝉联全球第一。但与国外成熟汽车市场相比,中国汽车市场仍存在巨大的发展前景,外界普遍预测,五年内中国汽车工业
内腔填充法是通过机械搅拌将Fe3O4粉末引入纤维内腔,而原位复合法是磁性纳米粒子通过化学过程形成沉淀沉积在纤维内腔中。通过X射线衍射证实了原位复合法合成出了Fe3O4磁性粒
本文首先论述了生产社会化的产生和发展不是源于资本主义生产方式, 而是源于科技革命。由于科学技术现已发展到知识经济时代, 生产社会化的特征也有了本质性的改变, 具体表现在
本文论述了提高数学课堂效率,实施创新能力培养的重要性和时代性,并着重从四个方面浅谈了本人对怎样培养学生学习兴趣的认识和在教学中总结出的经验做法。
被称为人的“第二肌肤”的服装,人们对其使用价值、服装机能、服装舒适性问题一直是服装领域研究、开发的首要问题。
以双酚A(BPA),碳酸二丁酯(DBC)为原料,制备出双酚A单丁基碳酸酯(I)和双酚A二丁基碳酸酯(II),用核磁共振波谱表征其结构。通过I、II的熔融自缩聚及II与BPA酯交换反应合成了双
<正>根据我国的具体国情及社会主义现代化建设不同阶段的要求,改革开放以来我国民族理论政策先后形成了"邓小平民族理论"、"三个代表"重要思想的民族理论以及十六大以来我们
当前中国正处于传统乡村景观向现代乡村景观过渡的关键时期,飞速的城市化进程为乡村的景观环境带来前所未有的挑战。为了改善和恢复乡村景观环境、营造美好的乡村生活环境,文