Algoritmi

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.

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 
        

• » Invia tramite EMAIL
• » Versione per la STAMPA
• » Le vostre opinioni