metropolis算法

❶ 求用fortran语言解释 伊辛模型 程序，最好在程序中注释出用蒙特卡洛方法的五个步骤在哪里

//我是按C/C++的注释格式写的，对此Fortran90程序直接改成!就可以了

programmain1
	parameter(n=50)
	integeri,j
	commona(n,n)
	realt,m,z
//设置输出文件
	open(1,file='out.dat')
//初始化Ising格点，均赋值为1

	do10i=1,n
	do20j=1,n
	a(i,j)=1
20continue
10continue

//体系按温度改变(0.1~4.0)进行MonteCarlo模拟
do30t=0.1,4.0,0.1
//体系的预热阶段，在此期间不采样
	do40,i=1,5000
	callabc(t)//调用metropolis算法子例程
40continue
//预热结束，可以对体系进行采样了
m=0.0
	do50,k=1,1000
	callabc(t)
	if(mod(k,5).eq.0)then
	z=0.0
	do60,i=1,n
	do70,j=1,n
		z=z+a(i,j)
70continue
60continue
m=m+z/(n*n)//计算体系的磁化强度M
	endif
50continue
m=m/200
	write(1,200)t,m//记录M随着T的变化曲线，向文件输出
	write(*,200)t,m//向屏幕输出
30continue

200format(2x,f6.4,8x,f10.6)
close(1)
	end

//实现Metropolis算法

	subroutineabc(t)
	parameter(n=50)
	commona(n,n)
	do10,i=1,n
	do20,j=1,n

//x1,x2,x3,x4是当前格点(i,j)的四个近邻
	x1=a(i+1,j)
	x2=a(i,j+1)
	x3=a(i-1,j)
	x4=a(i,j-1)
	if(i.eq.1)x3=a(n,j)
	if(i.eq.n)x1=a(1,j)
	if(j.eq.1)x4=a(i,n)
	if(j.eq.n)x2=a(i,1)


//h表示体系尝试翻转前后的Hamiltonian(能量)变化
	h=2*a(i,j)*(x1+x2+x3+x4)
	w=exp(-h/t)
	b=rand()
//按Metropolis规则，判断是否接受当前格点的尝试翻转
	if(w.gt.b)a(i,j)=-a(i,j)
20continue
10continue
end

//产生0~1均匀分布的随机数
	functionrand()
	realrand
	integerseed,c1,c2,c3
	parameter(c1=29,c2=217,c3=1024)
	dataseed/0/
	seed=mod(seed*c1+c2,c3)
	rand=real(seed)/c3
	end

❷ 如何理解贝叶斯估计

根据贝叶斯公式，进行统计推断，
在垃圾邮件分类方面应用很广，方法简单，具有很好的稳定性和健壮性

❸ Matlab下Metropolis-Hastings算法问题的求助

这个是M文件中的函数啊,只有运行在主界面并且这样运行: floyd(a)(把a输入括号这里才行) 单独运行当然没有定义a啊 %floyd.m %采用floyd算法计算图a中每对顶点最短路 %d是矩离矩阵 %r是路由矩阵 function [d,r]=floyd(a) n=size(a,1);

❹ 什么是mcmc算法

1.MCMC方法主要是为了解决有些baysian推断中参数期望E(f(v)|D)不能直接计算得到的问题的。
其中v是要估计的参数，D是数据观察值
2. Markov chain monte carlo概念包含了两部分：markov chain 和monte carlo integration。
2.1. 首先是monte carlo integration： monte carlo integration是利用采样的方法解决参数期望不能直接计算解决的问题的：即根据v的后验概率密度函数对v进行n次随机采样，计算n个f(v)，然后将这n个值求平均。根据大数定理当n足够大并且采样服从独立原则的时候，该值趋向于期望的真实值。但是当v的后验概率函数很难得到的时候该方法并不适用。
2.2 而在此基础上产生的 Markov chain monte carlo，虽然也是通过采样的方法进行的，却将马尔科夫链的概念引进来。它的想法是这样的：如果某条马尔科夫链具有irrecible 和aperiodic的特性的时候，该马尔科夫链具有一个唯一的静态点，即Pt=Pt-1;因此当马尔科夫链足够长后（设为N），产生的值会收敛到一个恒定的值（m）。这样对f（v）产生马尔科夫链，在N次之后f（v）的值收敛于恒定的值m，一般假设n>N后，f(v)服从N(m,scale)的正态分布。
即当n足够大的时候，用马尔科夫链产生的f(v)相当于在N(m,scale)独立抽样产生的值。
3. 如何产生具有这样特性的马尔科夫链：主要的方法是M-H算法
M-H算法有两部分组成：1. 根据条件概率密度函数，抽样得到下一个时间点的参数值Vt+1；2计算产生的这个值的接受概率a。如果a有显着性，就接受抽样得到的值，否则下一时间点的值保持不变。
在1中引入了proposal distribution的概念。在参数取值连续的情况下，后一个时间点的值服从一个分布，而这个分布函数只和前一个时间点的值有关q（.|Vt）
a的计算这里不贴了，一般都一样。重要的计算完a之后，如果决定接受采样获得的值还是保持原来值不变。在这个问题上，一般的处理方法是假设a服从0~1的均匀分布，每次采样计算a后都从U(0,1)中随机抽样一个值a'，如果a>=a'则接受抽样的值，否则保持原来的值不变。
根据q（.|Vt）的不同，M-H算法又有不同的分类：
3.1 Metroplis algorithm:在这里假设q(Vt+1|Vt)=q(|Vt+1 - Vt|)因此a被化解。该方法叫做random walk metropolis。
3.2 independence sampler:在这个算法里，假设q(Vt+1|Vt)=q(Vt+1)
对于多参数的情况，既可以同时产生多向量的马尔科夫链，又可以对每个参数分别进行更新：即，如果有h个参数需要估计，那么，在每次迭代的时候，分h次每次更新一个参数。
4. Gibbs抽样，H-M算法的一个变体

❺ 贴代码】R语言用metropolis-hastings算法产生4个自由度的t分布的均值,采用如下建议分布:1,N(0,1);2,t(2)

没太懂你说的建议分布是什么意思。比如N(0,1)，是说X^{n+1}=X^{n}+X（X~N(0,1)）还是直接X^{n+1}~N(0,1)？我当是后一种了哈。

python">#1.N(0,1)
set.seed(2015)#setrandomseed
N<-1e4#samplesize
x1<-rep(NA,N)
x1[1]<-rnorm(1)#oranotherinitialguess
for(nin2:N){
x1[n]<-rnorm(1)
log.rho<-dt(x1[n],df=4,log=TRUE)-dt(x1[n-1],df=4,log=TRUE)+dnorm(x1[n-1],log=TRUE)-dnorm(x1[n],log=TRUE)
if(runif(1)>exp(log.rho))x1[n]<-x1[n-1]
}
mean(x1)

#2.t(2)
N<-1e4
x2<-rep(NA,N)
x2[1]<-rt(1,df=2)#oranotherinitialguess
for(nin2:N){
x2[n]<-rt(1,df=2)
log.rho<-dt(x2[n],df=4,log=TRUE)-dt(x2[n-1],df=4,log=TRUE)+dt(x2[n-1],df=2,log=TRUE)-dt(x2[n],df=2,log=TRUE)
if(runif(1)>exp(log.rho))x2[n]<-x2[n-1]
}
mean(x2)

比较的话你就画histogram之类的吧，比如：

hist(x1,freq=FALSE,breaks=40);curve(dt(x,4),-4,4,add=TRUE)
hist(x2,freq=FALSE,breaks=40);curve(dt(x,4),-4,4,add=TRUE)

感觉t(2)要好一些，毕竟proposal不像N(0,1)那样保守。

❻ 什么是metropolis算法

算法是在有限步骤内求解某一问题所使用的一组定义明确的规则。通俗点说，就是计算机解题的过程。在这个过程中，无论是形成解题思路还是编写程序，都是在实施某种算法。前者是推理实现的算法，后者是操作实现的算法。

❼ 蒙特卡罗模拟

蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程，反复生成时间序列，计算参数估计量和统计量，进而研究其分布特征的方法。具体的，当系统中各个单元的可靠性特征量已知，但系统的可靠性过于复杂，难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时，可用随机模拟法近似计算出系统可靠性的预计值；随着模拟次数的增多，其预计精度也逐渐增高。由于涉及到时间序列的反复生成，蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的，因此只是在近些年才得到广泛推广。
蒙特卡洛(Monte Carlo)模拟这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。
蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时，可以用计算机模拟的方法产生抽样结果，根据抽样计算统计量或者参数的值；随着模拟次数的增多，可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。
蒙特卡洛模拟法求解步骤
应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。
解题步骤如下：
1.根据提出的问题构造一个简单、适用的概率模型或随机模型，使问题的解对应于该模型中随机变量的某些特征（如概率、均值和方差等），所构造的模型在主要特征参量方面要与实际问题或系统相一致
2 .根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，方可进行随机模拟试验。
3. 根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。
4.按照所建立的模型进行仿真试验、计算，求出问题的随机解。
5. 统计分析模拟试验结果，给出问题的概率解以及解的精度估计。

蒙特卡洛模拟法的应用领域
蒙特卡洛模拟法的应用领域主要有：
1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统，得到某些参数或重要指标。
2.蒙特卡洛积分:利用随机数列计算积分，维数越高，积分效率越高。
3.MCMC:这是直接应用蒙特卡洛模拟方法的推广，该方法中随机数的产生是采用的马尔科夫链形式。

❽ 到底什么是蒙特卡罗仿真方法

蒙特卡罗仿真原理
蒙特卡罗(MonteCarlo)方法，又称随机抽样或统计模拟方法，泛指所有基于统计采样进行数值计算的方法。在第二次世界大战期间，美国参与“曼哈顿计划’’的几位科学家Stanislaw Ulam，John Von Neumann 和 N．Metropolis等首先将这种方法用于解决原子弹研制中的一个关键问题。后来N．Metropolis用驰名世界的赌城－－－摩纳哥的MonteCarlo一来命名这种方法，为它蒙上了一层神秘色彩。随着现代计算机技术的飞速发展，蒙特卡罗方法已经在统计物理、经济学、社会学甚至气象学等方面的科学研究中发挥了极其重要的作用，将蒙特卡罗方法用于仿真即为蒙特卡罗仿真。蒙特卡罗方法适用于两类问题，第一类是本身就具有随机性的问题，第二类是能够转化为概率模型进行求解的确定性问题。
※蒙特卡罗方法求解问题的一般步骤
用蒙特卡罗方法求解问题一般包括构造或描述概率过程、从已知概率分布抽样和建立估计量三个步骤。
构造或描述概率过程实际上就是建立随机试验模型，构造概率过程是对确定性问题而言的，描述概率过程是对随机性问题而言的，不同的问题所需要建立的随机试验模型各不相同。
所谓的从已知概率分布抽样指的是随机试验过程，随机模型中必要包含某些已知概率分布的随机变量或随机过程作为输入，进行随机试验的过程就是对这些随机变量的样本或随机过程的样本函数作为输入产生相应输出的过程，因此通常被称为对已知概率分布的抽样。如何产生已知分布的随机变量或随机过程是蒙特卡罗方法中的一个关键问题。
最后一个步骤是获得估计量，蒙特卡罗方法所得到的问题的解总是对真实解的一个估计，本身也是一个随机变量，这个随机变量是由随机试验模型输出通过统计处理得到的。

❾ 如何在matlab中使用metropolis-hasting算法

MH算法也是一种基于模拟的MCMC技术，一个很重要的应用是从给定的概率分布中抽样。主要原理是构造了一个精妙的Markov链，使得该链的稳态 是你给定的概率密度。它的好处，不用多说，自然是可以对付数学形式复杂的概率密度。有人说，单维的MH算法配上Gibbs Sampler几乎是“无敌”了。今天试验的过程中发现，MH算法想用好也还不简单，里面的转移参数设定就不是很好弄。即使用最简单的高斯漂移项，方差的确定也是个头疼的问题，需要不同问题不同对待，多试验几次。当然你也可以始终选择“理想”参数。还是拿上次的混合高斯分布来做模拟，模拟次数为500000次的时候，概率分布逼近的程度如下图。虽然几个明显的"峰"已经出来了，但是数值上还是 有很大差异的。估计是我的漂移方差没有选好。感觉还是inverse sampling好用，迭代次数不用很多，就可以达到相当的逼近程度。

试了一下MH算法，R Code：

p=function(x,u1,sig1,u2,sig2){
(1/3)*(1/(sqrt(2*pi)*15)*exp(-0.5*(x-70)^2/15^2)+1/(sqrt(2*pi)*11)*exp(-0.5*(x+80)^2/11^2)+1/(sqrt(2*pi)*sig1)*exp(-0.5*(x-u1)^2/sig1^2)+1/(sqrt(2*pi)*sig2)*exp(-0.5*(x-u2)^2/sig2^2))
}

MH=function(x0,n){
x=NULL
x[1] = x0
for (i in 1:n){
x_can= x[i]+rnorm(1,0,3.25)
d= p(x_can,10,30,-10,10)/p(x[i],10,30,-10,10)
alpha= min(1,d)
u=runif(1,0,1)
if (u<alpha){
x[i+1]=x_can}
else{
x[i+1]=x[i]
}
if (round(i/100)==i/100) print(i)
}
x
}
z=MH(10,99999)
z=z[-10000]
a=seq(-100,100,0.2)

plot(density(z),col=1,main='Estimated Density',ylim=c(0,0.02),lty=1)
points(a, p(a,10,30,-10,10),pch='.',col=2,lty=2)
legend(60,0.02,c("True","Sim (MH)"),col=c(1,2),lty=c(1,2))