Bonjour,

Voilà je suis un cours de Monte Carlo suivant les chaînes de Markov et nous avons à répondre à ces questions :

a) Using the Metropolis algorithm, write a code that generates random numbers
according to the Gauss law : p(x) = A exp( -x^2/2). Is the constant A pertinent ?

b) Analyze the influence of , the maximum amplitude of the jumps, on the
acceptance rate and on the equilibration time. For this purpose, you may com-
pute the variance < x^2 >, that should converge to 1

VOICI MON CODE SOUS FORTRAN :

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
program Markov
real(kind=8) ksi, lambda, lambda_i, x0, x_trial,p0, p_trial,t, beta, sum_sqr_x
lambda = 0.001
print *,' Nbre de Pas ?'
read *, n_steps
x0 = 0.
p0 = exp(-x0**2/2)
sum_sqr_x = x0**2
 
do k = 1, n_steps
 
!--- Definition de la fonction p(x) ----------------------
 
  call random_number(ksi)
  lambda_i = lambda*(2*ksi-1)
  x_trial = x0 + lambda_i
  p_trial = exp(-x_trial**2/2)  
  t = min (1.0, p_trial/p0)  
 
  call random_number(beta)
  if (beta < t ) then
     x0 = x_trial
     p0 = p_trial      
  end if
 
  sum_sqr_x = sum_sqr_x + x0**2
 
end do
sum_sqr_x = sum_sqr_x / n_steps
print *, '<x**2 > = ', sum_sqr_x
end program MARKOV
À partir de là, tout se passe bien ; tout marche. Mais quand il s’agit de tracer (grâce à la fonction "plot") ...the ”asymptotic” probability density... je n'y arrive pas :/ Auriez-vous des solutions ?

Merci d'avances,
Mes amitiés à mes confrères les programmeurs