Questo sito contribuisce alla audience di

Algoritmo di Metropolis-Hastings

L'algoritmo di Metropolis-Hastings serve a generare dei numeri x1, x2, .., xn che presentano una distribuzione p(x) fissata a priori.

Il metodo si basa sulla generazione di numeri di ‘test’ che vengono accettati
o rigettati in modo da ottenere la distribuzione voluta. Il metodo sarà
presentato nel caso di una sola variabile random continua. Il metodo può essere
facilmente esteso al caso di distribuzioni di probabilità P(x1, x2, …, xN) di
un numero qualsiasi di variabili.

target=_BLANK>CONTINUA

Di seguito è riportato un codice fortran per l’implemetazione dell’Algoritmo
di Metropolis-Hastings:

c
********************************

implicit
real*8(a-h,o-z)

c **** Specify
parameters:
c **** Maximum allowed displacement magnitude in one
step
parameter(dxmax=0.1)
c ****
Number of trials
parameter(Ncycle=10000,
Nprint=Ncycle/10)
c **** Temperature

parameter(T=0.1, beta=1./T)
c **** Random number generator’s
seed

parameter(iseed=1234567)

c **** set
counters
fav =
0.0d0
f2av =
0.0d0

c *** set seed
for the random numbers generator (notice, that the name srand
c *** may be
different on different architectures )

call srand(iseed)

c *** Start
from the origin

xold=0.0

fold=f(x)

c *** Start MC
cycle (assuming rand() returns random real*1 from 0 to
1)

do icycle = 1,
Ncycle
c *** get new position, and
energy

xnew=xold+(rand()-0.5)*dxmax

fnew=f(xnew)
w = fnew
- fold
c **** do the Metroplis acceptance
criterium
if( (
w.lt.0.d0 ).or.( exp(-beta*w).gt.rand() )) then
c **** accept the move,
upgrade xold and
fold

xold=xnew

fold=fnew
endif
c
**** increment averages

fav = fav + fold

f2av= f2av + fold*fold

c **** print
out the running
average

if(mod(icycle,Nprint).eq.1)
then

frunav = fav /
icycle

f2runav = f2av /
icycle

f2runav =
sqrt((f2runav-frunav*frunav)/icycle)

write(6,2) icycle, frunav, f2runav
2
format(1x,i8, “cycles Average Energy = “,e12.5,

> ” Std.Dev “,
e12.5)

endif

enddo

c *** Get the
average and the standard deviation

fav =
fav / Ncycle
f2av = f2av /
Ncycle
f2av =
sqrt((f2av-fav*fav)/Ncycle)

c *** Write out
the result

write(6,1)
Ncycle, fav, f2av
1 format(1x,i8, “cycles
Average Energy = “,e12.5,

> ” Std.Dev “, e12.5)

c *** That’s
all folks

end

c ****
Potential
real*8 function
f(x)
real*8 x

f=x*x

return
end