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